Large-scale Sparse Linear Systems Research Articles

Compressed basis GMRES on high-performance graphics processing units

Krylov methods provide a fast and highly parallel numerical tool for the iterative solution of many large-scale sparse linear systems. To a large extent, the performance of practical realizations of these methods is constrained by the communication bandwidth in current computer architectures, motivating the investigation of sophisticated techniques to avoid, reduce, and/or hide the message-passing costs (in distributed platforms) and the memory accesses (in all architectures). This article leverages Ginkgo’s memory accessor in order to integrate a communication-reduction strategy into the (Krylov) GMRES solver that decouples the storage format (i.e., the data representation in memory) of the orthogonal basis from the arithmetic precision that is employed during the operations with that basis. Given that the execution time of the GMRES solver is largely determined by the memory accesses, the cost of the datatype transforms can be mostly hidden, resulting in the acceleration of the iterative step via a decrease in the volume of bits being retrieved from memory. Together with the special properties of the orthonormal basis (whose elements are all bounded by 1), this paves the road toward the aggressive customization of the storage format, which includes some floating-point as well as fixed-point formats with mild impact on the convergence of the iterative process. We develop a high-performance implementation of the “compressed basis GMRES” solver in the Ginkgo sparse linear algebra library using a large set of test problems from the SuiteSparse Matrix Collection. We demonstrate robustness and performance advantages on a modern NVIDIA V100 graphics processing unit (GPU) of up to 50% over the standard GMRES solver that stores all data in IEEE double-precision.

A GPU implementation of the PCG method for large-scale image-based finite element analysis in heterogeneous periodic media

Image-based physics simulations in heterogeneous media at a microscopic scale are a growing trend in various fields within and around Scientific Computing. We consider the scope of Numerical Homogenization, where micro-computed tomography is used to obtain digital models of physical samples, naturally leading to pixel and voxel-based solutions. In this context, the Finite Element Method (FEM) is commonly employed to solve the governing differential equations, via a system of algebraic equations. As image dimensions increase, the memory allocation due to the matrix associated with the FEM quickly becomes unfeasible, even in sparse format. Assembly-free strategies are adopted to reduce memory usage, with the caveat of increased computational cost. The Preconditioned Conjugate Gradient (PCG) method is widely employed to solve this sort of large-scale sparse linear systems , and is fitting to be adapted for assembly-free implementations. This work focuses on a massively parallel PCG solver applied to finite element analyses of heat conduction and linear elasticity on image-based models. Memory-efficiency is one of our main concerns, in an attempt to make feasible the employment of personal-use GPUs for large-scale simulations. The resulting solver is validated with an analytical benchmark, and by comparing the obtained results for a microtomographic model of a cast iron sample against experimental values found in the literature. Time and memory metrics are presented and discussed. It is shown that the developed program allows for homogenization studies of nearly 500 million degrees-of-freedom to be conducted in personal computers equipped with CUDA-enabled devices of 8 GB RAM, taking seconds or a few minutes per system solution with the PCG method. Up to 400× speed-up was observed in comparison to an analogous solver running in a 16-thread CPU. Our GPU implementation makes it possible to conduct, in a matter of minutes, homogenization studies that would take hours, or even days, in personal CPUs.

Computational techniques for [formula omitted] optimal frequency-limited model order reduction of large-scale sparse linear systems

We consider the problem of frequency limited H2 optimal model order reduction for large-scale sparse linear systems. A set of first-order H2 optimality conditions are derived for the frequency limited model order reduction problem. These conditions involve the solution of two frequency limited Sylvester equations that are known to be computationally complex. We discuss a framework for solving these matrix equations efficiently. The idea is also extended to the frequency limited H2 optimal model order reduction of index-1 descriptor systems. Numerical experiments are carried out using Python programming language and the results are presented to demonstrate the approximation accuracy and computational efficiency of the proposed technique.

Algebraic Multigrid Block Preconditioning for Multi-Group Radiation Diffusion Equations

The paper focuses on developing and studying efficient block preconditioners based on classical algebraic multigrid for the large-scale sparse linear systems arising from the fully coupled and implicitly cell-centered finite volume discretization of multi-group radiation diffusion equations, whose coefficient matrices can be rearranged into the $(G+2)\times(G+2)$ block form, where $G$ is the number of energy groups. The preconditioning techniques are based on the monolithic classical algebraic multigrid method, physical-variable based coarsening two-level algorithm and two types of block Schur complement preconditioners. The classical algebraic multigrid is applied to solve the subsystems that arise in the last three block preconditioners. The coupling strength and diagonal dominance are further explored to improve performance. We use representative one-group and twenty-group linear systems from capsule implosion simulations to test the robustness, efficiency, strong and weak parallel scaling properties of the proposed methods. Numerical results demonstrate that block preconditioners lead to mesh- and problem-independent convergence, and scale well both algorithmically and in parallel.

Two-Level Nyström--Schur Preconditioner for Sparse Symmetric Positive Definite Matrices

Randomized methods are becoming increasingly popular in numerical linear algebra. However, few attempts have been made to use them in developing preconditioners. Our interest lies in solving large-scale sparse symmetric positive definite linear systems of equations where the system matrix is preordered to doubly bordered block diagonal form (for example, using a nested dissection ordering). We investigate the use of randomized methods to construct high quality preconditioners. In particular, we propose a new and efficient approach that employs Nystr\"om's method for computing low rank approximations to develop robust algebraic two-level preconditioners. Construction of the new preconditioners involves iteratively solving a smaller but denser symmetric positive definite Schur complement system with multiple right-hand sides. Numerical experiments on problems coming from a range of application areas demonstrate that this inner system can be solved cheaply using block conjugate gradients and that using a large convergence tolerance to limit the cost does not adversely affect the quality of the resulting Nystr\"om--Schur two-level preconditioner.

Comment on “Efficient Secure Outsourcing of Large-Scale Sparse Linear Systems of Equations”

The scheme [IEEE Trans. Big Data, 4 (1), 2018, 26-39] is flawed because the sparsity of coefficient matrix is neglected. In the discussed scenario, we find it is unnecessary to outsource such a problem because the client can solve it locally. Even if the matrix is not sparse, the proposed paradigm which requires a great number of interactions, is rarely adopted because of the possibly incurred massive communication cost.

αSetup-AMG: an adaptive-setup-based parallel AMG solver for sequence of sparse linear systems

The algebraic multigrain (AMG) is one of the most frequently used algorithms for the solution of large-scale sparse linear systems in many realistic simulations of science and engineering applications. However, as the concurrency of supercomputers increasing, the AMG solver increasingly leads to poor parallel scalability due to its coarse-level construction in the setup phase. In this paper, to improve the parallel scalability of the traditional AMG to solve the sequence of sparse linear systems arising from PDE-based simulations, we propose a new AMG procedure called αSetup-AMG based on an adaptive setup strategy. The main idea behind αSetup-AMG is the introduction of a setup condition in the coarsening process so that the setup is constructed as it needed instead of constructing in advance via an independent phase in the traditional procedure. As a result, αSetup-AMG requires fewer setup cost and level numbers for the sequence of linear systems. The numerical results on thousands of cores for a radiation hydrodynamics simulation in the inertial confinement fusion (ICF) application show the significant improvement in the efficiency of the αSetup-AMG solver.

Efficient Secure Outsourcing of Large-Scale Sparse Linear Systems of Equations

Solving large-scale sparse linear systems of equations (SLSEs) is one of the most common and fundamental problems in big data, but it is very challenging for resource-limited users. Cloud computing has been proposed as a timely, efficient, and cost-effective way of solving such expensive computing tasks. Nevertheless, one critical concern in cloud computing is data privacy. Specifically, clients’ SLSEs usually contain private information that should remain hidden from the cloud for ethical, legal, or security reasons. Many previous works on secure outsourcing of linear systems of equations (LSEs) have high computational complexity, and do not exploit the sparsity in the LSEs. More importantly, they share a common serious problem, i.e., a huge number of memory I/O operations. This problem has been largely neglected in the past, but in fact is of particular importance and may eventually render those outsourcing schemes impractical. In this paper, we develop an efficient and practical secure outsourcing algorithm for solving large-scale SLSEs, which has low computational and memory I/O complexities and can protect clients’ privacy well. We implement our algorithm on Amazon Elastic Compute Cloud, and find that the proposed algorithm offers significant time savings for the client (up to 74 percent) compared to previous algorithms.

Insertion of PETSc in the OpenFOAM Framework

OpenFOAM is a widely used open source framework for simulation in several areas of computational fluid dynamics and engineering. As a partial differential equation (PDE)-based framework, OpenFOAM suffers from a performance bottleneck in solving large-scale sparse linear systems of equations. To address the problem, this article proposes a novel OpenFOAM-PETSc framework by inserting PETSc, a dedicated numerical solving package, into the OpenFOAM to speed up the process of solving linear equation systems. The design of the OpenFOAM-PETSc framework is described, and the implementation of an efficient matrix conversion algorithm is given as a case study. Validation tests on a high-performance computing cluster show that OpenFOAM-PETSc reduces the time of solving PDEs by about 27% in the lid-driven cavity flow case and by more than 50% in flow around the cylinder case in comparison with OpenFOAM, without compromising the scalability. In addition, this article also gives a preliminary performance analysis of different numerical solution methods, which may provide guidelines for further optimizations.

An adaptive AMG preconditioning strategy for solving large-scale sparse linear systems

A series of sparse linear systems must be solved in applications that are based on the implicit solution of time-dependent partial differential equations (PDEs). Preconditioned iterative methods are usually employed to solve such sparse linear systems. AMG is one of the most popular preconditioners in real applications. However, it results in poor parallel scalability, owing to its setup phase. In this paper, by utilizing the differences and similarities in property among the systems in series, an adaptive AMG preconditioning strategy is presented to improve the parallel scalability. The results obtained for a radiation hydrodynamics computation within an ICF simulation demonstrate the efficiency and improvement of the adaptive strategy. For a typical model, the new strategy improves the parallel efficiency from 47\% to 61\%, and reduces the CPU time from 19.7 h to 14.5 h.

Performance modeling of hyper-scale custom machine for the principal steps in block Wiedemann algorithm

Solving large-scale sparse linear systems over GF(2) plays a key role in fluid mechanics, simulation and design of materials, petroleum seismic data processing, numerical weather prediction, computational electromagnetics, and numerical simulation of unclear explosions. Therefore, developing algorithms for this issue is a significant research topic. In this paper, we proposed a hyper-scale custom supercomputer architecture that matches specific data features to process the key procedure of block Wiedemann algorithm and its parallel algorithm on the custom machine. To increase the computation, communication, and storage performance, four optimization strategies are proposed. This paper builds a performance model to evaluate the execution performance and power consumption for our custom machine. The model shows that the optimization strategies result in a considerable speedup, even three times faster than the fastest supercomputer, TH2, while consuming less power.

Scalable implicit incompressible resistive MHD with stabilized FE and fully-coupled Newton–Krylov-AMG

The computational solution of the governing balance equations for mass, momentum, heat transfer and magnetic induction for resistive magnetohydrodynamics (MHD) systems can be extremely challenging. These difficulties arise from both the strong nonlinear, nonsymmetric coupling of fluid and electromagnetic phenomena, as well as the significant range of time- and length-scales that the interactions of these physical mechanisms produce. This paper explores the development of a scalable, fully-implicit stabilized unstructured finite element (FE) capability for 3D incompressible resistive MHD. The discussion considers the development of a stabilized FE formulation in context of the variational multiscale (VMS) method, and describes the scalable implicit time integration and direct-to-steady-state solution capability. The nonlinear solver strategy employs Newton–Krylov methods, which are preconditioned using fully-coupled algebraic multilevel preconditioners. These preconditioners are shown to enable a robust, scalable and efficient solution approach for the large-scale sparse linear systems generated by the Newton linearization. Verification results demonstrate the expected order-of-accuracy for the stabilized FE discretization. The approach is tested on a variety of prototype problems, that include MHD duct flows, an unstable hydromagnetic Kelvin–Helmholtz shear layer, and a 3D island coalescence problem used to model magnetic reconnection. Initial results that explore the scaling of the solution methods are also presented on up to 128K processors for problems with up to 1.8B unknowns on a CrayXK7.

The simpler GMRES method combined with finite volume method for simulating viscoelastic flows on triangular grid

An efficient solver integrating the restarted simpler generalized minimal residual method (SGMRES(m)) with finite volume method (FVM) on triangular grid is developed to simulate the viscoelastic fluid flows. In particular, the SGMRES(m) solver is used to solve the large-scale sparse linear systems, which arise from the course of FVM on triangular grid for modeling the Newtonian and the viscoelastic fluid flows. To examine the performance of the solver for the nonlinear flow equations of viscoelastic fluids, we consider two types of numerical tests: the Newtonian flow past a circular cylinder, and the Oldroyd-B fluid flow in a planar channel and past a circular cylinder. It is shown that the numerical results obtained by the SGMRES(m) are consistent with the analytical solutions or empirical values. By comparing CPU time of different solvers, we find our solver is a highly efficient one for solving the flow equations of viscoelastic fluids.

Analysis of Communication-to-Computation Based on Asymptotic Size for Iterative Methods

Iterative method are one of the most efficient numerical algorithms to solve large-scale sparse linear systems arising from scientific computing.The parallel scalability of iterative method can be measured by the communication-to-computation(CtC) during the iterative process.The CtC is high for many iterative methods,as a result,coarse-grain parallelism is needed in order to obtain expected scalability.However,fine-grain parallelism is required for the complicacy of architecture increasing with multi/many-cores.In this paper,we introduce the concept of asymptotic size,which defined as the low bound of the problem size such that satisfying the speedup condition that parallel speedup is more than 1.We hope the asymptotic size can be used to describe the CtC and the ability of fine-grain parallelism for iterative method.Moreover,the theoretical prediction formula of asymptotic size is obtained based on the following parameters: the sparse and communication pattern of matrix,communication parameters of machine,and combination of basic operations of iterative methods.Using asymptotic size,the CtC is analyzed for three popular iterative methods,including Jacobi,CG,BiCGSTAB,on a MPP machine with 128 Double Quad-core nodes.The performance results are given for both MPI-Only and MPI/OpenMP programming model,which show the usefulness of the asymptotic size for describing the CtC of iterative methods.For MPI-Only case,we also give the comparison of the prediction results and experiment results,which show the validation of the formula of asymptotic size.Finally,future research topics for improving the scalability of the iterative methods on more power supercomputers also discussed based on the analysis conclusions.

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 CPU–GPU hybrid approach for the unsymmetric multifrontal method

Multifrontal is an efficient direct method for solving large-scale sparse and unsymmetric linear systems. The method transforms a large sparse matrix factorization process into a sequence of factorizations involving smaller dense frontal matrices. Some of these dense operations can be accelerated by using a graphic processing unit (GPU). We analyze the unsymmetric multifrontal method from both an algorithmic and implementational perspective to see how a GPU, in particular the NVIDIA Tesla C2070, can be used to accelerate the computations. Our main accelerating strategies include (i) performing BLAS on both CPU and GPU, (ii) improving the communication efficiency between the CPU and GPU by using page-locked memory, zero-copy memory, and asynchronous memory copy, and (iii) a modified algorithm that reuses the memory between different GPU tasks and sets thresholds to determine whether certain tasks be performed on the GPU. The proposed acceleration strategies are implemented by modifying UMFPACK, which is an unsymmetric multifrontal linear system solver. Numerical results show that the CPU–GPU hybrid approach can accelerate the unsymmetric multifrontal solver, especially for computationally expensive problems.

Exploiting thread-level parallelism in the iterative solution of sparse linear systems

We investigate the efficient iterative solution of large-scale sparse linear systems on shared-memory multiprocessors. Our parallel approach is based on a multilevel ILU preconditioner which preserves the mathematical semantics of the sequential method in ILUPACK. We exploit the parallelism exposed by the task tree corresponding to the nested dissection hierarchy (task parallelism), employ dynamic scheduling of tasks to processors to improve load balance, and formulate all stages of the parallel PCG method conformal with the computation of the preconditioner to increase data reuse. Results on a CC-NUMA platform with 16 processors reveal the parallel efficiency of this solution.

Towards a scalable fully-implicit fully-coupled resistive MHD formulation with stabilized FE methods

This paper explores the development of a scalable, nonlinear, fully-implicit stabilized unstructured finite element (FE) capability for 2D incompressible (reduced) resistive MHD. The discussion considers the implementation of a stabilized FE formulation in context of a fully-implicit time integration and direct-to-steady-state solution capability. The nonlinear solver strategy employs Newton–Krylov methods, which are preconditioned using fully-coupled algebraic multilevel preconditioners. These preconditioners are shown to enable a robust, scalable and efficient solution approach for the large-scale sparse linear systems generated by the Newton linearization. Verification results demonstrate the expected order-of-accuracy for the stabilized FE discretization. The approach is tested on a variety of prototype problems, including both low-Lundquist number (e.g., an MHD Faraday conduction pump and a hydromagnetic Rayleigh–Bernard linear stability calculation) and moderately-high Lundquist number (magnetic island coalescence problem) examples. Initial results that explore the scaling of the solution methods are presented on up to 4096 processors for problems with up to 64M unknowns on a CrayXT3/4. Additionally, a large-scale proof-of-capability calculation for 1 billion unknowns for the MHD Faraday pump problem on 24,000 cores is presented.

Symmetric Indefinite Preconditioners for Saddle Point Problems with Applications to PDE-Constrained Optimization Problems

We consider large scale sparse linear systems in saddle point form. A natural property of such indefinite 2-by-2 block systems is the positivity of the (1,1) block on the kernel of the (2,1) block. Many solution methods, however, require that the positivity of the (1,1) block is satisfied everywhere. To enforce the positivity everywhere, an augmented Lagrangian approach is usually chosen. However, the adjustment of the involved parameters is a critical issue. We will present a different approach that is not based on such an explicit augmentation technique. For the considered class of symmetric and indefinite preconditioners, assumptions are presented that lead to symmetric and positive definite problems with respect to a particular scalar product. Therefore, conjugate gradient acceleration can be used. An important class of applications are optimal control problems. It is typical for such problems that the cost functional contains an extra regularization parameter. For control problems with elliptic state equations and distributed control, a special indefinite preconditioner for the discretized problem is constructed, which leads to convergence rates of the preconditioned conjugate gradient method that are not only independent of the mesh size but also independent of the regularization parameter. Numerical experiments are presented for illustrating the theoretical results.

Using the parallel algebraic recursive multilevel solver in modern physical applications

This paper discusses the application of a few parallel preconditioning techniques, which are collected in a recently developed suite of codes Parallel Algebraic Recursive Multilevel Solver (pARMS), to tackling large-scale sparse linear systems arising from real-life applications. In particular, we study the effect of different algorithmic variations and parameter choices on the overall performance of the distributed preconditioners in pARMS by means of numerical experiments related to a few realistic applications. These applications include magnetohydrodynamics, nonlinear acoustic field simulation, and tire design.