Multifrontal Algorithm Research Articles

A study of concurrent multi-frontal solvers for modern massively parallel architectures

Leveraging Trace Theory, we investigate the efficient parallelization of direct solvers for large linear equation systems. Our focus lies on a multi-frontal algorithm, and we present a methodology for achieving near-optimal scheduling on modern massively parallel machines. By employing trace theory with Diekert Graphs and Foata Normal Form, we rigorously validate the correctness of our proposed solution. To establish a strong link between the mesh and elimination tree of the multi-frontal solver, we conduct extensive testing on matrices derived from the Finite Element Method (FEM). Furthermore, we assess the performance of computations on both GPU and CPU platforms, employing practical implementation strategies.

VoxImp: Impedance Extraction Simulator for Voxelized Structures

An impedance extractor, called VoxImp, is proposed to compute the impedances of the structures discretized by voxels. VoxImp iteratively solves the volume-surface integral equations discretized by a carefully selected set of basis functions. During iterative solution, matrix-vector multiplications are accelerated by the fast Fourier transform, while the rapid convergence of the iterative solution at resonant frequencies is ensured by a novel sparse preconditioner. The memory requirement of the sparse preconditioner is reduced by a sparse LU decomposition obtained by a multifrontal algorithm with compressed frontal matrices. The overall memory requirement of the VoxImp is further reduced by the Tucker decomposition. VoxImp’s accuracy, efficiency, and capability are demonstrated through impedance extraction of various voxelized structures, including an RF coil array discretized by more than four million voxels and fifteen million panels and analyzed on a commodity desktop computer.

Three-Dimensional Modeling of Marine Controlled Source Electromagnetic Using a New High-Order Finite Element Method With an Absorption Boundary Condition

The marine controlled source electromagnetic method (MCSEM) is widely used in marine oil, gas exploration and deep structures investigation because of its low cost and high efficiency. In this paper, a high-precision forward algorithm for 3D modeling for MCSEM is proposed. Unstructured meshes are used to discretize the model domain, and local mesh refinement techniques are applied, which is conducive to simulating complex terrains and targets. The application of electric dipole discrete technology can load wire excitation sources with arbitrarily complex spatial shapes, which can simulate more realistic electromagnetic field distribution of marine controllable sources electromagnetic in actual exploration. Absorption boundary conditions based on real number and exponential stretching techniques are introduced to improve the numerical solution accuracy. To improve the efficiency of solving, the open-source massively parallel solver (MUMPS) based on the multifrontal algorithm is used to solve the finite element equations. Finally, some typical geoelectric models are designed to verify the correctness and effectiveness of the proposed algorithm. The calculation results show that higher order finite elements bring about higher accuracy. In addition, the loading of absorption boundary conditions is simple and effective than conventional Dirichlet boundary conditions.

High performance sparse multifrontal solvers on modern GPUs

We have ported the numerical factorization and triangular solve phases of the sparse direct solver STRUMPACK to GPU. STRUMPACK implements sparse LU factorization using the multifrontal algorithm, which performs most of its operations in dense linear algebra operations on so-called frontal matrices of various sizes. Our GPU implementation off-loads these dense linear algebra operations, as well as the sparse scatter–gather operations between frontal matrices. For the larger frontal matrices, our GPU implementation relies on vendor libraries such as cuBLAS and cuSOLVER for NVIDIA GPUs and rocBLAS and rocSOLVER for AMD GPUs. For the smaller frontal matrices we developed custom CUDA and HIP kernels to reduce kernel launch overhead. Overall, high performance is achieved by identifying submatrix factorizations corresponding to sub-trees of the multifrontal assembly tree which fit entirely in GPU memory. The multi-GPU setting uses SLATE (Software for Linear Algebra Targeting Exascale) as a modern GPU-aware replacement for ScaLAPACK. On 4 nodes of SUMMIT the code runs ∼10× faster when using all 24 V100 GPUs compared to when it only uses the 168 POWER9 cores. On 8 SUMMIT nodes, using 48 V100 GPUs, the sparse solver reaches over 50TFlop/s. Compared to SuperLU, on a single V100, for a set of 17 matrices our implementation is faster for all but one matrix, and is on average 5× (median 4×) faster

Inexact proximal Newton methods for self-concordant functions

We analyze the proximal Newton method for minimizing a sum of a self-concordant function and a convex function with an inexpensive proximal operator. We present new results on the global and local convergence of the method when inexact search directions are used. The method is illustrated with an application to L1-regularized covariance selection, in which prior constraints on the sparsity pattern of the inverse covariance matrix are imposed. In the numerical experiments the proximal Newton steps are computed by an accelerated proximal gradient method, and multifrontal algorithms for positive definite matrices with chordal sparsity patterns are used to evaluate gradients and matrix-vector products with the Hessian of the smooth component of the objective.

A fast, memory efficient and robust sparse preconditioner based on a multifrontal approach with applications to finite‐element matrices

SummaryIn this article, we introduce a fast, memory efficient and robust sparse preconditioner that is based on a direct factorization scheme for sparse matrices arising from the finite‐element discretization of elliptic partial differential equations. We use a fast (but approximate) multifrontal approach as a preconditioner and use an iterative scheme to achieve a desired accuracy. This approach combines the advantages of direct and iterative schemes to arrive at a fast, robust, and accurate preconditioner. We will show that this approach is faster (∼2×) and more memory efficient (∼2–3×) than a conventional direct multifrontal approach. Furthermore, we will demonstrate that this preconditioner is both faster and more effective than other preconditioners such as the incomplete LU preconditioner. Specific speedups depend on the matrix size and improve as the size of the matrix increases. The preconditioner can be applied to both structured and unstructured meshes in a similar manner. We build on our previous work and utilize the fact that dense frontal and update matrices, in the multifrontal algorithm, can be represented as hierarchically off‐diagonal low‐rank matrices. Using this idea, we replace all large dense matrix operations in the multifrontal elimination process with O(N) hierarchically off‐diagonal low‐rank operations to arrive at a faster and more memory efficient factorization scheme. We then use this direct factorization method at low accuracies as a preconditioner and apply it to various real‐life engineering test cases. Copyright © 2016 John Wiley & Sons, Ltd.

Quasi-Optimal Elimination Trees for 2D Grids with Singularities

We construct quasi-optimal elimination trees for 2D finite element meshes with singularities. These trees minimize the complexity of the solution of the discrete system. The computational cost estimates of the elimination process model the execution of the multifrontal algorithms in serial and in parallel shared-memory executions. Since the meshes considered are a subspace of all possible mesh partitions, we call these minimizers quasi-optimal. We minimize the cost functionals using dynamic programming. Finding these minimizers is more computationally expensive than solving the original algebraic system. Nevertheless, from the insights provided by the analysis of the dynamic programming minima, we propose a heuristic construction of the elimination trees that has costONelog⁡Ne, whereNeis the number of elements in the mesh. We show that this heuristic ordering has similar computational cost to the quasi-optimal elimination trees found with dynamic programming and outperforms state-of-the-art alternatives in our numerical experiments.

Application of supernodal sparse factorization and inversion to the estimation of (co)variance components by residual maximum likelihood.

We demonstrated that supernodal techniques were more efficient than traditional methods for factorization and inversion of a coefficient matrix of mixed model equations (MME), which are often required in residual maximum likelihood (REML). Supernodal left-looking and inverse multifrontal algorithms were employed for sparse factorization and inversion, respectively. The approximate minimum degree or multilevel nested dissection was used for ordering. A new computer package, Yet Another MME Solver (YAMS), was developed and compared with FSPAK with respect to computing time and size of temporary memory for 13 test matrices. The matrices were produced by fitting animal models to dairy data and by using simulations from sire, sire-maternal grand sire, maternal and dominance models for phenotypic data and animal model for genomic data. The order of matrices ranged from 32,840 to 1,048,872. The YAMS software factorized and inverted the matrices up to 13 and 10 times faster than FSPAK, respectively, when an appropriate ordering strategy was applied. The YAMS package required at most 282 MB and 512 MB of temporary memory for factorization and inversion, respectively. Processing time per iteration in average information REML was reduced, using YAMS. The YAMS package is freely available on request by contacting the corresponding author.

The closed-form solution of vector finite element integral equations for three dimensional electromagnetic analysis

In this paper, a set of closed-form formulas for vector Finite Element Method (FEM) to analyze three dimensional electromagnetic problems is presented on the basis of Gaussian quadrature integration scheme. By analyzing the open region problems, the first-order Absorbing Boundary Condition (ABC) is considered as the truncation boundary condition and the equation is carried out in a closed-form. Based on the formulas, the hybrid Expanded Cholesky Method (ECM) and MultiFrontal algorithm (MF) is applied to solve finite element equations. Using the closed-form solution, the electromagnetic field of three dimensional targets can be studied sententiously and accurately. Simulation results show that the presented formulas are successfully and concise, which can be easily used to analyze three dimensional electromagnetic problems.

Numerical methods for on-line power system load flow analysis

Newton-Raphson method is the most widely accepted load flow solution algorithm. However LU factorization remains a computationally challenging task to meet the real-time needs of the power system. This paper proposes the application of very fast multifrontal direct linear solvers for solving the linear system sub-problem of power system real-time load flow analysis by utilizing the state-of-the-art algorithms for ordering and preprocessing. Additionally the unsymmetric multifrontal method for LU factorization and highly optimized Intel® Math Kernel Library BLAS has been used. Two state-of-the-art multifrontal algorithms for unsymmetric matrices namely UMFPACK V5.2.0 and sequential MUMPS 4.8.3 (“Multifrontal Massively Parallel Solver”) are customized for the AC power system Newton-Raphson based load flow analysis. The multifrontal solvers are compared against the state-of-the-art sparse Gaussian Elimination based HSL sparse solver MA48. This study evaluates the performance of above multifrontal solvers in terms of number of factors, computational time, number of floating-point operations and memory, in the context of load flow solution on nine systems including very large real power systems. The results of the performance evaluation are reported. The proposed method achieves significant reduction in computational time.

Frequency-Domain Linearized Euler Model for Turbomachinery Noise Radiation Through Engine Exhaust

A numerical model for the exhaust noise radiation problem is presented. In the model, it is assumed that an incoming wave is propagating through the exhaust nozzle, or the fan duct, and radiating outside. The near-field propagation is based on the solution of the linearized Euler equations in the frequency domain: for each wave number, a linearized Euler problem is solved using a finite element method on unstructured grids for arbitrarily shaped axisymmetric geometries. The frequency-domain approach enables the suppression of the Kelvin-Helmholtz instability waves. Moreover, each single calculation, limited to a single frequency, is well suited to the exhaust noise radiation problem in which the incoming wave can be treated as a superposition of elementary duct modes. To reduce the memory requirements, a continuous Galerkin formulation with linear triangular and quadrangular elements is employed and the global matrix inversion is performed with a direct solver based on a parallel memory distributed multifrontal algorithm for sparse matrices. The acoustic near field is then radiated in the far field using the formulation of Ffowcs Williams and Hawkings. Numerical calculations for a validation test case, the Munt problem, and two turbomachinery configurations are compared with analytical solutions and experimental data.

An out-of-core sparse Cholesky solver

Direct methods for solving large sparse linear systems of equations are popular because of their generality and robustness. Their main weakness is that the memory they require usually increases rapidly with problem size. We discuss the design and development of the first release of a new symmetric direct solver that aims to circumvent this limitation by allowing the system matrix, intermediate data, and the matrix factors to be stored externally. The code, which is written in Fortran and called HSL_MA77, implements a multifrontal algorithm. The first release is for positive-definite systems and performs a Cholesky factorization. Special attention is paid to the use of efficient dense linear algebra kernel codes that handle the full-matrix operations on the frontal matrix and to the input/output operations. The input/output operations are performed using a separate package that provides a virtual-memory system and allows the data to be spread over many files; for very large problems these may be held on more than one device. Numerical results are presented for a collection of 30 large real-world problems, all of which were solved successfully.

AN EFFICIENT APPROACH FOR MULTIFRONTAL ALGORITHM TO SOLVE NON-POSITIVE-DEFINITE FINITE ELEMENT EQUATIONS IN ELECTROMAGNETIC PROBLEMS

A new method called Expanded Cholesky Method (ECM) is proposed in this paper. The method can be used to decompose sparse symmetric non-positive-deflnite flnite element (FEM) matrices. There are some advantages of the ECM, such as low storage, simplicity and easy parallelization. Based on the method, multifrontal (MF) algorithm is applied in non-positive-deflnite FEM computation. Numerical results show that the hybrid ECM/MF algorithm is stable and efiective. In comparison with Generalized Minimal Residual Method (GMRES) in FEM electromagnetic computation, hybrid ECM/MF technology has distinct advantages in precision. The proposed method can be used to calculate a class of non-positive- deflnite electromagnetic problems.

THE COMBINATION OF BCGSTAB WITH MULTIFRONTAL ALGORITHM TO SOLVE FEBI-MLFMA LINEAR SYSTEMS ARISING FROM INHOMOGENEOUS ELECTROMAGNETIC SCATTERING PROBLEMS

The hybrid flnite-element/boundary-integral method (FEBI) combined with the multilevel fast multipole algorithm (MLFMA) has been applied to model the three-dimensional scattering problems of inhomogeneous media. The stabilized Bi-conjugate gradi- ent (BCGATAB) iterative solver based on the inner-looking algorithm is proposed to solve the flnal FEBI linear system, and the multifrontal algorithm combined with the approximate minimal degree permutation (AMD) is used for the LU decomposition of the FEM matrix. The ac- curacy and e-ciency of the combined algorithm has been validated in the flnal of the paper. Numerical results show that the proposed method can greatly improve the e-ciency of FEBI for scattering prob- lems of inhomogeneous media.

An efficient out‐of‐core multifrontal solver for large‐scale unsymmetric element problems

AbstractIn many applications where the efficient solution of large sparse linear systems of equations is required, a direct method is frequently the method of choice. Unfortunately, direct methods have a potentially severe limitation: as the problem size grows, the memory needed generally increases rapidly. However, the in‐core memory requirements can be limited by storing the matrix and its factors externally, allowing the solver to be used for very large problems. We have designed a new out‐of‐core package for the large sparse unsymmetric systems that arise from finite‐element problems. The code, which is called HSL_MA78, implements a multifrontal algorithm and achieves efficiency through the use of specially designed code for handling the input/output operations and efficient dense linear algebra kernels. These kernels, which are available as a separate package called HSL_MA74, use high‐level BLAS to perform the partial factorization of the frontal matrices and offer both threshold partial and rook pivoting. In this paper, we describe the design of HSL_MA78 and explain its user interface and the options it offers. We also describe the algorithms used by HSL_MA74 and illustrate the performance of our new codes using problems from a range of practical applications. Copyright © 2008 John Wiley & Sons, Ltd.

A Shared- and distributed-memory parallel general sparse direct solver

An important recent development in the area of solution of general sparse systems of linear equations has been the introduction of new algorithms that allow complete decoupling of symbolic and numerical phases of sparse Gaussian elimination with partial pivoting. This enables efficient solution of a series of sparse systems with the same nonzero pattern but different coefficient values, which is a fairly common situation in practical applications. This paper reports on a shared- and distributed-memory parallel general sparse solver based on these new symbolic and unsymmetric-pattern multifrontal algorithms.

High-resolution finite element simulation of 4:1 planar contraction flow of viscoelastic fluid

In this work, we present high-resolution solutions for viscoelastic 4:1 planar contraction flow problems using a transient finite element method based on the fractional step method (FSM) and stabilization techniques (DEVSS-G/DG) with linear equal-order interpolation function. The Oldroyd-B model was used as the constitutive equation. A parallel multi-frontal algorithm was implemented to enhance computational speed and all solutions were obtained on a parallel machine. The vortex intensity and the re-attachment length of corner vortex show good mesh-convergent behavior and are compared with previous results from the literature. In particular, the present results are in good agreement with the predictions of the high-resolution finite volume method of Alves et al. [15]. This may be the first case that quantitative agreement is obtained between studies using different numerical methods for the benchmark problem of 4:1 planar contraction flow. As there has been little quantitative agreement in the previous investigations and only few simulation results with highly refined meshes exit, this study may well be regarded as accurate and meaningful in the sense that reasonable convergence is achieved for prediction of 4:1 planar contraction flow using transient finite element methods.

Analysis of waveguide problems using a relaxed iterative domain decomposition method combined with multifrontal algorithm

In this paper, an absorbing Fictitious Boundary Condition (FBC) is presented to generate an iterative Domain Decomposition Method (DDM) for analyzing waveguide problems. The relaxed algorithm is introduced to improve the iterative convergence. And the matrix equations are solved using the multifrontal algorithm. The resulting CPU time is greatly reduced. Finally, a number of numerical examples are given to illustrate its accuracy and efficiency.

Improved Symbolic and Numerical Factorization Algorithms for Unsymmetric Sparse Matrices

We present algorithms for the symbolic and numerical factorization phases in the direct solution of sparse unsymmetric systems of linear equations. We have modified a classical symbolic factorization algorithm for unsymmetric matrices to inexpensively compute minimal elimination structures. We give an efficient algorithm to compute a near-minimal data-dependency graph for unsymmetric multifrontal factorization that is valid irrespective of the amount of dynamic pivoting performed during factorization. Finally, we describe an unsymmetric-pattern multifrontal algorithm for Gaussian elimination with partial pivoting that uses the task- and data-dependency graphs computed during the symbolic phase. These algorithms have been implemented in WSMP---an industrial strength sparse solver package---and have enabled WSMP to significantly outperform other similar solvers. We present experimental results to demonstrate the merits of the new algorithms.

Efficient Methods for Out-of-Core Sparse Cholesky Factorization

We consider the problem of sparse Cholesky factorization with limited main memory. The goal is to efficiently factor matrices whose Cholesky factors essentially fill the available disk storage, using very little memory (as little as 16 Megabytes (MBytes)). This would enable very large industrial problems to be solved with workstations of very modest cost. We consider three candidate algorithms. Each is based on a partitioning of the matrix into panels. The first is a robust, out-of-core multifrontal method that keeps the factor, the stack, and the large frontal matrices on disk. The others are left-looking methods. We find that straightforward implementations of all of them suffer from excessive disk I/O for large problems that arise in interior-point algorithms for linear programming. We introduce several improvements to these simple out-of-core methods and find that a left-looking method that nevertheless uses the multifrontal algorithm for portions of the matrix (subtrees of the supernodal elimination tree whose multifrontal stack fits in memory) is very effective. With 32 Mbytes of main memory, it achieves over 77% of its in-core performance on all but one of our 12 test matrices (67% in that one case), even though the size of the factor is, in all cases, hundreds of millions or even billions of bytes.