Sparse LU Factorization Research Articles

A graphics processing unit accelerated sparse direct solver and preconditioner with block low rank compression

We present the GPU implementation efforts and challenges of the sparse solver package STRUMPACK. The code is made publicly available on github with a permissive BSD license. STRUMPACK implements an approximate multifrontal solver, a sparse LU factorization which makes use of compression methods to accelerate time to solution and reduce memory usage. Multiple compression schemes based on rank-structured and hierarchical matrix approximations are supported, including hierarchically semi-separable, hierarchically off-diagonal butterfly, and block low rank. In this paper, we present the GPU implementation of the block low rank (BLR) compression method within a multifrontal solver. Our GPU implementation relies on highly optimized vendor libraries such as cuBLAS and cuSOLVER for NVIDIA GPUs, rocBLAS and rocSOLVER for AMD GPUs and the Intel oneAPI Math Kernel Library (oneMKL) for Intel GPUs. Additionally, we rely on external open source libraries such as SLATE (Software for Linear Algebra Targeting Exascale), MAGMA (Matrix Algebra on GPU and Multi-core Architectures), and KBLAS (KAUST BLAS). SLATE is used as a GPU-capable ScaLAPACK replacement. From MAGMA we use variable sized batched dense linear algebra operations such as GEMM, TRSM and LU with partial pivoting. KBLAS provides efficient (batched) low rank matrix compression for NVIDIA GPUs using an adaptive randomized sampling scheme. The resulting sparse solver and preconditioner runs on NVIDIA, AMD and Intel GPUs. Interfaces are available from PETSc, Trilinos and MFEM, or the solver can be used directly in user code. We report results for a range of benchmark applications, using the Perlmutter system from NERSC, Frontier from ORNL, and Aurora from ALCF. For a high frequency wave equation on a regular mesh, using 32 Perlmutter compute nodes, the factorization phase of the exact GPU solver is about 6.5× faster compared to the CPU-only solver. The BLR-enabled GPU solver is about 13.8× faster than the CPU exact solver. For a collection of SuiteSparse matrices, the STRUMPACK exact factorization on a single GPU is on average 1.9× faster than NVIDIA’s cuDSS solver.

Exactly Solving Sparse Rational Linear Systems via Roundoff-Error-Free Cholesky Factorizations

Exactly solving sparse symmetric positive definite (SPD) linear systems is a key problem in mathematics, engineering, and computer science. This paper derives two new sparse roundoff-error-free (REF) Cholesky factorization algorithms which exactly solve sparse SPD linear systems $A {x} = {b}$, where $A \in \mathbb{Q}^{n x n}$ and ${x}, {b} \in {Q}^{n x p}$. The key properties of these factorizations are that (1) they exclusively use integer-arithmetic and (2) in the bit-complexity model, they solve the linear system $A {x} = {b}$ in time proportional to the cost of the integer-arithmetic operations. Namely, the overhead related to data structures and ancillary operations (those not strictly required to perform the factorization) is subsumed by the cost of the integer-arithmetic operations that are essential/intrinsic to the factorization. Notably, to date our algorithms are the only exact algorithm for solving SPD linear systems with this asymptotically efficient complexity bound. Computationally, we show that the novel factorizations are faster than both sparse rational-arithmetic LDL and sparse exact LU factorization. Altogether, the derived sparse REF Cholesky factorizations present a framework to solve any rational SPD linear system exactly and efficiently.

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

Performance impact of precision reduction in sparse linear systems solvers.

It is well established that reduced precision arithmetic can be exploited to accelerate the solution of dense linear systems. Typical examples are mixed precision algorithms that reduce the execution time and the energy consumption of parallel solvers for dense linear systems by factorizing a matrix at a precision lower than the working precision. Much less is known about the efficiency of reduced precision in parallel solvers for sparse linear systems, and existing work focuses on single core experiments. We evaluate the benefits of using single precision arithmetic in solving a double precision sparse linear system using multiple cores. We consider both direct methods and iterative methods and we focus on using single precision for the key components of LU factorization and matrix–vector products. Our results show that the anticipated speedup of 2 over a double precision LU factorization is obtained only for the very largest of our test problems. We point out two key factors underlying the poor speedup. First, we find that single precision sparse LU factorization is prone to a severe loss of performance due to the intrusion of subnormal numbers. We identify a mechanism that allows cascading fill-ins to generate subnormal numbers and show that automatically flushing subnormals to zero avoids the performance penalties. The second factor is the lack of parallelism in the analysis and reordering phases of the solvers and the absence of floating-point arithmetic in these phases. For iterative solvers, we find that for the majority of the matrices computing or applying incomplete factorization preconditioners in single precision provides at best modest performance benefits compared with the use of double precision. We also find that using single precision for the matrix–vector product kernels provides an average speedup of 1.5 over double precision kernels. In both cases some form of refinement is needed to raise the single precision results to double precision accuracy, which will reduce performance gains.

Algorithmic Advancements and a Comparative Investigation of Left and Right Looking Sparse LU Factorization on GPU Platform for Circuit Simulation

Sparse LU factorization is a key tool in the solution of large linear set of algebraic equations encompassing a wide range of computing applications. Recent advances in this field exploit the massively parallel architecture of the GPUs via left-looking algorithm (LLA) and right-looking algorithm (RLA). In this paper, <i>adaptive cluster</i> mode is proposed to improve the state-of-the-art in LLA for GPU platforms. The proposed method takes into consideration of varying sparsity at different levels during cluster mode execution, to adaptively configure the GPU block size and the number of parallel columns. The new refinements for LLA are also integrated with the dynamic parallelism that is available in modern GPU architectures. The paper also provides a comprehensive performance comparison of the LLA and hybrid RLA along with state-of-the-art advances on the same GPU platform. The results indicate that, when implemented with similar refinements and on a same platform, LLA provides better performance compared to the hybrid-RLA. The results would be useful to the scientific community while making decision on adopting LLA or RLA algorithms for sparse LU factorization.

GPU-Accelerated Adaptive PCBSO Mode-Based Hybrid RLA for Sparse LU Factorization in Circuit Simulation

LU factorization is extensively used in engineering and scientific computations for solution of large set of linear equations. Particularly, circuit simulators rely heavily on sparse version of LU factorization for solution involving circuit matrices. One of the recent advances in this field is exploiting the emerging computing platform of graphics processing units (GPUs) for parallel and sparse LU factorization. In this article, following contributions are made to advance the state of the art in hybrid right-looking algorithm (RLA): 1) a novel GPU kernel based on parallel column and block size optimization (PCBSO) is developed for adaptively allocating the block size while optimizing the number of columns for parallel execution based on the size of their associated submatrices at every level. The proposed approach helps to minimize the resource contention and to improve the computational performance and 2) an algorithm is developed to enable the execution of the new adaptive mode with dynamic parallelism. Also, a comprehensive performance comparison using a set of benchmark circuit examples is presented. The results indicate that, the proposed advancements can improve the results of state-of-the-art right looking sparse LU factorization in GPU by <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$1.54\times $ </tex-math></inline-formula> (Arithmetic Mean).

An on-node scalable sparse incomplete LU factorization for a many-core iterative solver with Javelin

We present a scalable incomplete LU factorization to be used as a preconditioner for solving sparse linear systems with iterative methods in the package called Javelin. Javelin allows for improved parallel factorization on shared-memory many-core systems by packaging the coefficient matrix into a format that allows for high performance sparse matrix-vector multiplication and sparse triangular solves with minimal overheads. The framework achieves these goals by using a collection of traditional permutations, point-to-point thread synchronizations, tasking, and segmented prefix scans in a conventional compressed sparse row (CSR) format. Moreover, this framework stresses the importance of co-designing dependent tasks, such as sparse factorization and triangular solves, on highly-threaded architectures. We compare our method to the past distributed methods for incomplete factorization (Aztec) and current multithreaded packages (WSMP) in order to demonstrate the importance of having highly threaded factorizations on many-core systems. Using these changes, traditional fill-in and drop tolerance methods can be used, while still being able to have observed speedups of up to ~ 42 × on 68 Intel Knights Landing cores and ~ 12 × on 14 Intel Haswell cores. Moreover, this work provides insight into how the new data-structure impacts iteration counts, and provides insight into future improvements, such as point to GPUs.

A Hierarchical Low-Rank Approximation Based Network Solver for EMT Simulation

In electromagnetic transient (EMT) simulation, 80–97% of the computational time is devoted to solving the network equations. A key observation is that the sub-matrix representing the interaction between two far-away groups of buses is usually sparse and can be approximated by a low-rank matrix. Based on this observation, we propose a novel low-rank approximation method which permits ${O(N \log N)}$ -time matrix-vector multiplication for each network solution time step. Comprehensive numerical studies are conducted on a 39-bus system and a 179-bus system from the literature, and large cases created from the two systems. The results demonstrate that the proposed approach is up to $2.8\times$ faster than the state-of-the-art sparse LU factorization based network solution, without compromising simulation accuracy. Since our low-rank approximation is highly parallelizable, further speedup may be possible.

A novel TDFEM scheme based on high‐performance parallel solver PaStiX

AbstractIn this paper, we focus on the theory and fulfillment of parallel time domain finite element method (TDFEM) with PaStiX (Parallel Sparse matriX package) for solving three‐dimensional electromagnetic (EM) radiation problems. The motivation behind the research is required to maintain numerical stability and achieve high speed for parallel sparse LU factorization on parallel platforms. Firstly, we overview the concept of an unconditional stable TDFEM scheme, and numerical optimization for ill‐conditioned equations obtained from the TDFEM with perfectly matched layer (PML) was presented to alleviate nonconvergent results. Secondly, we develop sparse linear system solvers that use DS factorization to solve the preconditioned system. Numerical results were obtained from the cases of the cavity resonance problem and a cone antenna. A noteworthy observation from the time‐consuming results in different solvers is that PaStiX is usually much faster compared with UMFPACK and SuperLU.

Research on the performance of multifrontal method in power system simulations

Efficiency of power system simulation has close relationship to the sparse LU factorization algorithm. Multifrontal method is a notable LU factorization member. This paper investigates the performance of multifrontal method in power system simulations. Performance of multifrontal method heavily relies on the frontal matrix size. Frontal matrix size can be represented by front size. Front size is subject to the treewidth of power network. A treewidth approximation method (TAM) is presented to approximate the treewidth of power network. Two derivatives of TAM are given to tackle the sparse matrices from power flow calculation and small signal stability eigenvalue analysis respectively. Simulation results verify the accuracy and efficiency of TAM and its derivatives. Further results indicate that most power networks have low treewidth, which is prone to degrading the performance of multifrontal method. Multifrontal method is able to gain its performance in power system simulations when the front size of the involved matrix is large enough.

Revisiting the BE SBS algorithm and applying it to solve torsion problems in composite bars: Robustness and efficiency study

The construction of truly competitive boundary-element (BE) codes, capable of analyzing real-life engineering problems, unavoidably requires the devising of coupling strategies to make the modeling of complex heterogeneous domains user-friendly. Thereby, high-performance algorithms for solving the highly sparse resulting system of equations are essential. Moreover, as discontinuous boundary elements are necessary to alleviate the modeling process of coupled domains, efficient (low-order) quadratures for integrating singular and nearly-singular fundamental kernels over the boundary elements must be implemented. This paper newly discusses the construction of the boundary-element subregion-by-subregion (BE SBS) technique based on a BEM formulation for torsion problems in general composite bars. One also presents details of the formulation of the Krylov solvers BiCG and BiCGSTAB-(l), embedded in the coupling algorithm. In addition, the BE SBS matrix structure itself is used to form an efficient sparse incomplete LU factorization (SILU) preconditioner to accelerate the iterative solution process. Torsion problems in bars with complex composite patterns (e.g. with many different materials) are analyzed to attest the efficiency and robustness of the whole boundary-element technique.

A communication-avoiding 3D algorithm for sparse LU factorization on heterogeneous systems

We propose a new algorithm to improve the strong scalability of right-looking sparse LU factorization on distributed memory systems. Our 3D algorithm for sparse LU uses a three-dimensional MPI process grid, exploits elimination tree parallelism, and trades off increased memory for reduced per-process communication. We also analyze the asymptotic improvements for planar graphs (e.g., those arising from 2D grid or mesh discretizations) and certain non-planar graphs (specifically for 3D grids and meshes). For a planar graph with n vertices, our algorithm reduces communication volume asymptotically in n by a factor of Ologn and latency by a factor of Ologn. For non-planar cases, our algorithm can reduce the per-process communication volume by 3× and latency by On13 times. In all cases, the memory needed to achieve these gains is a constant factor. We implemented our algorithm by extending the 2D data structure used in SuperLU_DIST. Our new 3D code achieves empirical speedups up to 27× for planar graphs and up to 3.3× for non-planar graphs over the baseline 2D SuperLU_DIST when run on 24,000 cores of a Cray XC30. We extend the 3D algorithm for heterogeneous architectures by adding the Highly Asynchronous Lazy Offload (Halo) algorithm for co-processor offload [44]. On 4096 nodes of a Cray XK7 with 32,768 CPU cores and 4096 Nvidia K20x GPUs, the 3D algorithm achieves empirical speedups up to 24× for planar graphs and 3.5× for non-planar graphs over the baseline 2D SuperLU_DIST with co-processor acceleration.

Basker: Parallel sparse LU factorization utilizing hierarchical parallelism and data layouts

Abstract Transient simulation in circuit simulation tools, such as SPICE and Xyce, depend on scalable and robust sparse LU factorizations for efficient numerical simulation of circuits and power grids. As the need for simulations of very large circuits grow, the prevalence of multicore architectures enable us to use shared memory parallel algorithms for such simulations. A parallel factorization is a critical component of such shared memory parallel simulations. We develop a parallel sparse factorization algorithm that can solve problems from circuit simulations efficiently, and map well to architectural features. This new factorization algorithm exposes hierarchical parallelism to accommodate irregular structure that arise in our target problems. It also uses a hierarchical two-dimensional data layout which reduces synchronization costs and maps to memory hierarchy found in multicore processors. We present an OpenMP based implementation of the parallel algorithm in a new multithreaded solver called Basker in the Trilinos framework. We present performance evaluations of Basker on the Intel SandyBridge and Xeon Phi platforms using circuit and power grid matrices taken from the University of Florida sparse matrix collection and from Xyce circuit simulation. Basker achieves a geometric mean speedup of 5.91× on CPU (16 cores) and 7.4× on Xeon Phi (32 cores) relative to state-of-the-art solver KLU. Basker outperforms Intel MKL Pardiso solver (PMKL) by as much as 30× on CPU (16 cores) and 7.5× on Xeon Phi (32 cores) for low fill-in circuit matrices. Furthermore, Basker provides 5.4× speedup on a challenging matrix sequence taken from an actual Xyce simulation.

A sparse stiff chemistry solver based on dynamic adaptive integration for efficient combustion simulations

A sparse stiff chemistry solver based on dynamic adaptive hybrid integration (AHI-S) is developed and demonstrated for efficient combustion simulations. In a previous study, a dynamic adaptive method for hybrid integration (AHI) was developed to speed up the time integration of chemically reacting flows with detailed chemistry. The AHI method solves the fast subcomponent of chemistry implicitly and the slow subcomponent of chemistry and transport explicitly, and it was shown that AHI is more accurate and efficient than the operator-splitting schemes when there are significant radical sources from the transport term. In the present study, the AHI method is first improved to minimize the number of nontrivial entries in the Jacobian. Sparse matrix techniques are further integrated into AHI to achieve high computational efficiency. The performance of the new AHI-S solver is investigated in constant-pressure auto-ignition systems using different mechanisms that consist of 9–2878 species. It is shown that the computational cost of the AHI-S solver is overall linearly proportional to the mechanism size and is comparable to that of evaluating reaction rates using CHEMKIN-II subroutines. The AHI-S solver achieves speed-up factors ranging from approximately 10, for the 9-species hydrogen mechanism, to approximately 3000, for the 2878-species biodiesel mechanism, compared with the fully implicit VODE solver with Jacobian evaluated through numerical perturbations and factorized with dense matrix operations. It is further found that for mechanisms with less than approximately 100 species, the time saving of AHI-S is primarily attributed to the reduced size of the implicit core of the governing equations, while for mechanisms with more than 100 species, the computational cost of VODE is dominated by the dense LU factorization, such that the time saving of AHI-S is mostly attributed to the sparse LU factorization. The AHI-S solver is then applied to unsteady perfectly stirred reactors involving extinction and re-ignition. Speed-up factors from 50 to 30,000 are achieved compared with the Strang splitting scheme with the chemistry substeps implicitly integrated with VODE, while speed-up factors of 10–100 are achieved compared with the Strang splitting scheme implemented with the sparse stiff LSODES solver. In the end, the performance of AHI-S is investigated in one-dimensional (1-D) unsteady freely propagating laminar premixed flames for a methane/air mixture, for which the time step size in AHI-S is limited by the fastest transport process. A speed-up factor of approximately 200 is achieved compared with the Strang splitting scheme for fixed time step sizes between 10−8s and 10−6s.

An algebraic multifrontal preconditioner that exploits the low‐rank property

SummaryWe present an algebraic structured preconditioner for the iterative solution of large sparse linear systems. The preconditioner is based on a multifrontal variant of sparse LU factorization used with nested dissection ordering. Multifrontal factorization amounts to a partial factorization of a sequence of logically dense frontal matrices, and the preconditioner is obtained if structured factorization is used instead. This latter exploits the presence of low numerical rank in some off‐diagonal blocks of the frontal matrices. An algebraic procedure is presented that allows to identify the hierarchy of the off‐diagonal blocks with low numerical rank based on the sparsity of the system matrix. This procedure is motivated by a model problem analysis, yet numerical experiments show that it is successful beyond the model problem scope. Further aspects relevant for the algebraic structured preconditioner are discussed and illustrated with numerical experiments. The preconditioner is also compared with other solvers, including the corresponding direct solver. Copyright © 2015 John Wiley &amp; Sons, Ltd.

GPU-Accelerated Sparse LU Factorization for Circuit Simulation with Performance Modeling

The sparse matrix solver by LU factorization is a serious bottleneck in Simulation Program with Integrated Circuit Emphasis (SPICE)-based circuit simulators. The state-of-the-art Graphics Processing Units (GPU) have numerous cores sharing the same memory, provide attractive memory bandwidth and compute capability, and support massive thread-level parallelism, so GPUs can potentially accelerate the sparse solver in circuit simulators. In this paper, an efficient GPU-based sparse solver for circuit problems is proposed. We develop a hybrid parallel LU factorization approach combining task-level and data-level parallelism on GPUs. Work partitioning, number of active thread groups, and memory access patterns are optimized based on the GPU architecture. Experiments show that the proposed LU factorization approach on NVIDIA GTX580 attains an average speedup of 7.02 $\times$ (geometric mean) compared with sequential PARDISO, and 1.55 $\times$ compared with 16-threaded PARDISO. We also investigate bottlenecks of the proposed approach by a parametric performance model. The performance of the sparse LU factorization on GPUs is constrained by the global memory bandwidth, so the performance can be further improved by future GPUs with larger memory bandwidth.

A Novel Factorized Sparse Approximate Inverse Preconditioner with Supernodes

Abstract Krylov methods preconditioned by Factorized Sparse Approximate Inverses (FSAI) are an efficient approach for the solution of symmetric positive definite linear systems on massively parallel computers. However, FSAI often suffers from a high set-up cost, especially in ill-conditioned problems. In this communication we propose a novel algorithm for the FSAI computation that makes use of the concept of supernode borrowed from sparse LU factorizations and direct methods.

The Use of Supernodes in Factored Sparse Approximate Inverse Preconditioning

In recent years the growing popularity of supercomputers has fostered the development of algorithms able to take advantage of the massive parallelism offered by multiple processors. Direct methods, though robust and computationally efficient, hardly exploit high degrees of parallelism. By contrast, Krylov methods preconditioned by Factored Sparse Approximate Inverses (FSAI) provide, at least in principle, a perfectly parallel approach but are often thwarted by an excessive set-up cost. In this paper we extend the concept of supernode from sparse LU factorizations to approximate inverses, and use it to accelerate the computation of an FSAI-type preconditioner. The numerical experiments on real-world problems show that the overall FSAI efficiency can be significantly increased while preserving its intrinsic parallelism.

Dynamic Multiple Work Stealing Strategy for Flexible Load Balancing

SUMMARY Lazy-task creation is an efficient method of overcoming the overhead of the grain-size problem in parallel computing. Work stealing is an effective load balancing strategy for parallel computing. In this paper, we present dynamic work stealing strategies in a lazy-task creation technique for efficient fine-grain task scheduling. The basic idea is to control load balancing granularity depending on the number of task parents in a stack. The dynamic-length strategy of work stealing uses run-time information, which is information on the load of the victim, to determine the number of tasks that a thief is allowed to steal. We compare it with the bottommost first work stealing strategy used in StackThread/MP, and the fixed-length strategy of work stealing, where a thief requests to steal a fixed number of tasks, as well as other multithreaded frameworks such as Cilk and OpenMP task implementations. The experiments show that the dynamic-length strategy of work stealing performs well in irregular workloads such as in UTS benchmarks, as well as in regular workloads such as Fibonacci, Strassen’s matrix multiplication, FFT, and Sparse-LU factorization. The dynamic-length strategy works better than the fixed-length strategy because it is more flexible than the latter; this strategy can avoid

Algorithm 907

KLU is a software package for solving sparse unsymmetric linear systems of equations that arise in circuit simulation applications. It relies on a permutation to Block Triangular Form (BTF), several methods for finding a fill-reducing ordering (variants of approximate minimum degree and nested dissection), and Gilbert/Peierls’ sparse left-looking LU factorization algorithm to factorize each block. The package is written in C and includes a MATLAB interface. Performance results comparing KLU with SuperLU, Sparse 1.3, and UMFPACK on circuit simulation matrices are presented. KLU is the default sparse direct solver in the Xyce TM circuit simulation package developed by Sandia National Laboratories.