Multifrontal Method Research Articles

Computational cost estimates for parallel shared memory isogeometric multi-frontal solvers

In this paper we present computational cost estimates for parallel shared memory isogeometric multi-frontal solvers. The estimates show that the ideal isogeometric shared memory parallel direct solver scales as O(p2log(N/p)) for one dimensional problems, O(Np2) for two dimensional problems, and O(N4/3p2) for three dimensional problems, where N is the number of degrees of freedom, and p is the polynomial order of approximation. The computational costs of the shared memory parallel isogeometric direct solver are compared with those corresponding to the sequential isogeometric direct solver, being the latest equal to O(Np2) for the one dimensional case, O(N1.5p3) for the two dimensional case, and O(N2p3) for the three dimensional case. The shared memory version significantly reduces both the scalability in terms of N and p. Theoretical estimates are compared with numerical experiments performed with linear, quadratic, cubic, quartic, and quintic B-splines, in one and two spatial dimensions.

Computational Complexity Analysis for Multifrontal Method With Limited-Depth of Assembly Tree

The computational complexity and memory requirement of multifrontal method is analyzed for solving finite element system of equations defined on 2- and 3-D regular meshes with respect to the number of fronts and the depth of assembly tree. Accurate estimation of operation counts obtained along with the relevant parameters for both local and multilevel global condensation phases is more suitable to practical applications than previous estimation. It, moreover, enables us to illustrate how the number of fronts, the depth of assembly tree, and the performance of the method are related, and how the multilevel multifrontal scheme reduces the total numerical operations intrinsically. The operation count, memory usage, and execution time of the method are experimentally measured by running a sequence of 2- and 3-D examples, the obtained results consistently follow the theoretical estimations.

Performance models and workload distribution algorithms for optimizing a hybrid CPU–GPU multifrontal solver

Problems that involve large and sparse linear systems are ubiquitous in scientific computing, and there are strong needs to accelerate the solution processes. Hybrid CPU–GPU systems have recently become a new platform trend with powerful computing capabilities. However, it is not clear how such systems can accelerate the solvers. We study how to make the best use of the CPU and the GPU to minimize the total time required to solve symmetric positive definite systems using the multifrontal method. We analyze the computation and communication costs of the multifrontal method on such hybrid systems to build up timing performance models. Workload distribution algorithms are proposed to determine if a frontal matrix should be factored on the CPU or on the GPU to minimize the total execution time of the overall computation. We provide theoretical analyses and numerical results to illustrate the characteristics and efficiency of the proposed algorithms. Because the performance models and workload distribution algorithms can accommodate different CPUs and GPUs adaptively, we expect the applicability and significance of these techniques to continue to grow as heterogeneous hardware and software evolve.

A study of shared-memory parallelism in a multifrontal solver

We introduce shared-memory parallelism in a parallel distributed-memory solver, targeting multi-core architectures. Our concern in this paper is pure shared-memory parallelism, although the work will also impact distributed-memory parallelism. Our approach avoids a deep redesign and fully benefits from the numerical kernels and features of the original code. We use performance models to exploit coarse-grain parallelism in an OpenMP environment while, at the same time, also relying on third-party optimized multithreaded libraries. In this context, we propose simple approaches to take advantage of NUMA architectures, and original optimizations to limit thread synchronization costs. The performance gains are analyzed in detail on test problems from various application areas. Although the studied code is a direct solver for sparse systems of linear equations, the contributions of this paper are more general and could be useful in a wider range of situations.

Parallel Computation of Entries of ${A}^{-1}$

In this paper, we are concerned about computing in parallel several entries of the inverse of a large sparse matrix. We assume that the matrix has already been factorized by a direct method and that the factors are distributed. Entries are efficiently computed by exploiting sparsity of the right-hand sides and the solution vectors in the triangular solution phase. We demonstrate that in this setting, parallelism and computational efficiency are two contrasting objectives. We develop an efficient approach and show its efficacy by runs using the MUMPS code that implements a parallel multifrontal method.

Hypergraph Grammar based Linear Computational Cost Solver for Three Dimensional Grids with Point Singularities

In this paper we present a hypergraph grammar based multi-frontal solver for three dimensional grids with point singularities. We show experimentally that the computational cost of the resulting solver algorithm is linear with respect to the number of degrees of freedom. We also propose a reutilization algorithm that enables to reuse LU factorizations over unrefined parts of the mesh when new local refinements are executed by the hypergraph grammar productions.

Dynamic Programming Algorithm for Generation of Optimal Elimination Trees for Multi-frontal Direct Solver Over H-refined Grids

In this paper we present a dynamic programming algorithm for finding optimal elimination trees for computational grids refined towards point or edge singularities. The elimination tree is utilized to guide the multi-frontal direct solver algorithm. Thus, the criterion for the optimization of the elimination tree is the computational cost associated with the multi-frontal solver algorithm executed over such tree. We illustrate the paper with several examples of optimal trees found for grids with point, isotropic edge and anisotropic edge mixed with point singularity. We show the comparison of the execution time of the multi-frontal solver algorithm with results of MUMPS solver with METIS library, implementing the nested dissection algorithm.

Graph Grammar based Multi-thread Multi-frontal Direct Solver with Galois Scheduler

In this paper, we present a multi-frontal solver algorithm for the adaptive finite element method expressed by graph grammar productions. The graph grammar productions construct first the binary elimination tree, and then process frontal matrices stored in distributed manner in nodes of the elimination tree. The solver is specialized for a class of one, two and three dimensional h refined meshes whose elimination tree has a regular structure. In particular, this class contains all one dimensional grids, two and three dimensional grids refined towards point singularities, two dimensional grids refined in an anisotropic way towards edge singularity as well as three dimensional grids refined in an anisotropic way towards edge or face singularities. In all these cases, the structure of the elimination tree and the structure of the frontal matrices are similar. The solver is implemented within the Galois environment, which allows parallel execution of graph grammar productions. We also compare the performance of the Galois implementation of our graph grammar based solver with the MUMPS solver.

확장 유한 요소 법을 적용한 RTM 공정 해석

ABSTRACT: Numerical simulation for Resin Transfer Molding manufacturing process is attempted by using theeXtended Finite Element Method (XFEM) combined with the level set method. XFEM allows to obtaining a goodnumerical precision of the pressure near the resin flow front, where its gradient is discontinuous. The enriched shapefunctions of XFEM are derived by using the level set values so as to correctly describe the interpolation with the resinflow front. In addition, the level set method is used to transport the resin flow front at each time step during themold filling. The level set values are calculated by an implicit characteristic Galerkin FEM. The multi-frontal solver ofIPSAP is adopted to solve the system. This work is validated by comparing the obtained results with analyticsolutions. Moreover, a localization method of XFEM and level set method is proposed to increase the computingefficiency. The computation domain is reduced to the small region near the resin flow front. Therefore, the totalcomputing time is strongly reduced by it. The efficiency test is made with a simple channel flow model. Severalapplication examples are analyzed to demonstrate ability of this method.초록

An O(N) algorithm for constructing the solution operator to 2D elliptic boundary value problems in the absence of body loads

The large sparse linear systems arising from the finite element or finite difference discretization of elliptic PDEs can be solved directly via, e.g., nested dissection or multifrontal methods. Such techniques reorder the nodes in the grid to reduce the asymptotic complexity of Gaussian elimination from O(N 2) to O(N 1.5) for typical problems in two dimensions. It has recently been demonstrated that the complexity can be further reduced to O(N) by exploiting structure in the dense matrices that arise in such computations (using, e.g., $\mathcal {H}$ -matrix arithmetic). This paper demonstrates that such accelerated nested dissection techniques become particularly effective for boundary value problems without body loads when the solution is sought for several different sets of boundary data, and the solution is required only near the boundary (as happens, e.g., in the computational modeling of scattering problems, or in engineering design of linearly elastic solids). In this case, a modified version of the accelerated nested dissection scheme can execute any solve beyond the first in O(N boundary) operations, where N boundary denotes the number of points on the boundary. Typically, N boundary ∼ N 0.5. Numerical examples demonstrate the effectiveness of the procedure for a broad range of elliptic PDEs that includes both the Laplace and Helmholtz equations.

A rapidly converging domain decomposition method for the Helmholtz equation

A new domain decomposition method is introduced for the heterogeneous 2-D and 3-D Helmholtz equations. Transmission conditions based on the perfectly matched layer (PML) are derived that avoid artificial reflections and match incoming and outgoing waves at the subdomain interfaces. We focus on a subdivision of the rectangular domain into many thin subdomains along one of the axes, in combination with a certain ordering for solving the subdomain problems and a GMRES outer iteration. When combined with multifrontal methods, the solver has near-linear cost in examples, due to very small iteration numbers that are essentially independent of problem size and number of subdomains. It is to our knowledge only the second method with this property next to the moving PML sweeping method.

Improving Multifrontal Methods by Means of Block Low-Rank Representations

Matrices coming from elliptic Partial Differential Equations (PDEs) have been shown to have a low-rank property: well defined off-diagonal blocks of their Schur complements can be approximated by low-rank products. Given a suitable ordering of the matrix which gives to the blocks a geometrical meaning, such approximations can be computed using an SVD or a rank-revealing QR factorization. The resulting representation offers a substantial reduction of the memory requirement and gives efficient ways to perform many of the basic dense algebra operations. Several strategies have been proposed to exploit this property. We propose a low-rank format called Block Low-Rank (BLR), and explain how it can be used to reduce the memory footprint and the complexity of direct solvers for sparse matrices based on the multifrontal method. We present experimental results that show how the BLR format delivers gains that are comparable to those obtained with hierarchical formats such as Hierarchical matrices (H matrices) and Hierarchically Semi-Separable (HSS matrices) but provides much greater flexibility and ease of use which are essential in the context of a general purpose, algebraic solver.

Efficient Structured Multifrontal Factorization for General Large Sparse Matrices

Rank structures provide an opportunity to develop new efficient numerical methods for practical problems, when the off-diagonal blocks of certain dense intermediate matrices have small (numerical) ranks. In this work, we present a framework of structured direct factorizations for general sparse matrices, including discretized PDEs on general meshes, based on the multifrontal method and hierarchically semiseparable (HSS) matrices. We prove the idea of replacing certain complex structured operations by fast simple ones performed on compact reduced matrix forms. Such forms result from the hierarchical factorization of a tree-structured HSS matrix in a ULV-type scheme, so that the tree structure is reduced into a single node, the root of the original tree. This idea is shown to be very useful in the partial ULV factorization of an HSS matrix (for quickly computing Schur complements) as well as the solution stage. These techniques are then built into the multifrontal method for sparse factorizations after nested dissection, so as to convert the intermediate dense factorizations into fast structured ones. This method keeps certain Schur complements dense so as to avoid complicated data assembly, and is much simpler and more general than some existing methods. In particular, if the matrix arises from the discretization of certain PDEs, the factorization costs roughly $O(n)$ flops in two dimensions, and roughly $O(n^{4/3})$ flops or less in three dimensions. The solution cost and memory are nearly $O(n)$ in both cases. These counts are obtained with an idea of rank relaxation, so that this method is more generally applicable to problems where the intermediate off-diagonal ranks are not small. We demonstrate the performance of the method with two- and three-dimensional discretized equations, as well as various examples from a sparse matrix collection. The ideas here are also useful in future developments of fast structured solvers.

A Parallel Sweeping Preconditioner for Heterogeneous 3D Helmholtz Equations

A parallelization of a sweeping preconditioner for 3D Helmholtz equations without large cavities is introduced and benchmarked for several challenging velocity models. The setup and application costs of the sequential preconditioner are shown to be O({\gamma}^2 N^{4/3}) and O({\gamma} N log N), where {\gamma}({\omega}) denotes the modestly frequency-dependent number of grid points per Perfectly Matched Layer. Several computational and memory improvements are introduced relative to using black-box sparse-direct solvers for the auxiliary problems, and competitive runtimes and iteration counts are reported for high-frequency problems distributed over thousands of cores. Two open-source packages are released along with this paper: "Parallel Sweeping Preconditioner (PSP)" and the underlying distributed multifrontal solver, "Clique".

Graph Grammar Based Multi-Frontal Solver for Isogeometric FEM Simulations On GPU

We present a graph grammar model of a multi-frontal direct solver for one and two dimensional isogeometric finite element method simulations with NVIDIA CUDA and perform numerical experiments for linear, quadratic, cubic, quartic and quintic B-splines. Finally, we estimate the computational cost of the for a parallel shared memory implementation: O(p^2log(N/p)) for 1D, O(Np^2) for 2D, as well as O(N^1.33 p^2) for 3D problems. We compare the parallel costs with the corresponding estimates for a standard sequential implementation: O((N/p)p^2) for 1D, O(N^1.5p^3) for 2D and O(N^2p^3) for 3D. We conclude the presentation with observation that computational cost of the shared memory direct solver scales like p^2 when we increase the global continuity of the isogeometric solution, which is an adventage with respect to sequential isogeometric solver scalability of the order of p^3. We conclude with a discussion on practical limitations of such GPU implementations, related to memory usage especially for 3D.

Grammar-Based Multi-Frontal Solver for One Dimensional Isogeometric Analysis with Multiple Right-Hand-Sides

This paper introduces a grammar-based model for developing a multi-thread multi-frontal parallel direct solver for one- dimensional isogeometric finite element method. The model includes the integration of B-splines for construction of the element local matrices and the multi-frontal solver algorithm. The integration and the solver algorithm are partitioned into basic indivisible tasks, namely the grammar productions, that can be executed squentially. The partial order of execution of the basic tasks is analyzed to provide the scheduling for the execution of the concurrent integration and multi-frontal solver algo- rithm. This graph grammar analysis allows for optimal concurrent execution of all tasks. The model has been implemented and tested on NVIDIA CUDA GPU, delivering logarithmic execution time for linear, quadratic, cubic and higher order B-splines. Thus, the CUDA implementation delivers the optimal performance predicted by our graph grammar analysis. We utilize the solver for multiple right hand sides related to the solution of non-stationary or inverse problems.

Fine-Grained Multithreading for the Multifrontal $QR$ Factorization of Sparse Matrices

The advent of multicore processors represents a disruptive event in the history of computer science as conventional parallel programming paradigms are proving incapable of fully exploiting their potential for concurrent computations. The need for different or new programming models clearly arises from recent studies which identify fine-granularity and dynamic execution as the keys to achieving high efficiency on multicore systems. This work presents an approach to the parallelization of the multifrontal method for the $QR$ factorization of sparse matrices specifically designed for multicore based systems. High efficiency is achieved through a fine-grained partitioning of data and a dynamic scheduling of computational tasks relying on a dataflow parallel programming model. Experimental results show that an implementation of the proposed approach achieves higher performance and better scalability than existing equivalent software.

Randomized Sparse Direct Solvers

We propose randomized direct solvers for large sparse linear systems, which integrate randomization into rank structured multifrontal methods. The use of randomization highly simplifies various essential steps in structured solutions, where fast operations on skinny matrix-vector products replace traditional complex ones on dense or structured matrices. The new methods thus significantly enhance the flexibility and efficiency of using structured methods in sparse solutions. We also consider a variety of techniques, such as some graph methods, the inclusion of additional structures, the concept of reduced matrices, information reuse, and adaptive schemes. The methods are applicable to various sparse matrices with certain rank structures. Particularly, for discretized matrices whose factorizations yield dense fill-in with some off-diagonal rank patterns, the factorizations cost about $O(n)$ flops in two dimensions (2D), and about $O(n)$ to $O(n^{4/3})$ flops in three dimensions (3D). The solution costs and memory sizes are nearly $O(n)$ in both 2D and 3D. These counts are obtained based on two optimization strategies and a sparse rank relaxation idea. The methods are especially useful for approximate solutions and preconditioning. Numerical tests on both discretized PDEs and more general problems are used to demonstrate the efficiency and accuracy. The ideas here also have the potential to be generalized to matrix-free sparse direct solvers based on matrix-vector multiplications in future developments.

IMF: An Incomplete Multifrontal $LU$-Factorization for Element-Structured Sparse Linear Systems

We propose an incomplete multifrontal $LU$-factorization (IMF) preconditioner that extends supernodal multifrontal methods to incomplete factorizations. It can be used as a preconditioner in a Krylov-subspace method to solve large-scale sparse linear systems with an element structure, e.g., those arising from a finite element discretization of a partial differential equation. The fact that the element matrices are dense is exploited to increase the computational performance and the robustness of the factorization through efficient partial pivoting. IMF is compared with the multilevel ARMS2, the level of fill-in ILU, and the threshold-based ILUTP preconditioners. Our experiments indicate that IMF is competitive with ARMS2 on saddle-point problems arising in the solution of the steady-state Navier--Stokes equation. Experiments with element-structured matrices arising from structural engineering applications, found in the University of Florida Sparse Matrix Collection, illustrate the robustness of IMF. Finally, the computational performance of IMF clearly surpasses that of the related ARMS2 preconditioner.

A bubble-inspired algorithm for finite element mesh partitioning

SUMMARY This paper presents a bubble-inspired algorithm for partitioning finite element mesh into subdomains. Differing from previous diffusion BUBBLE and Center-oriented Bubble methods, the newly proposed algorithm employs the physics of real bubbles, including nucleation, spherical growth, bubble–bubble collision, reaching critical state, and the final competing growth. The realization of foaming process of real bubbles in the algorithm enables us to create partitions with good shape without having to specify large number of artificial controls. The minimum edge cut is simply achieved by increasing the volume of each bubble in the most energy efficient way. Moreover, the order, in which an element is gathered into a bubble, delivers the minimum number of surface cells at every gathering step; thus, the optimal numbering of elements in each subdomain has naturally achieved. Because finite element solvers, such as multifrontal method, must loop over all elements in the local subdomain condensation phase and the global interface solution phase, these two features have a huge payback in terms of solver efficiency. Experiments have been conducted on various structured and unstructured meshes. The obtained results are consistently better than the classical kMetis library in terms of the edge cut, partition shape, and partition connectivity. Copyright © 2012 John Wiley & Sons, Ltd.