Solution Of Sparse Linear Systems Research Articles

Alya toward exascale: algorithmic scalability using PSCToolkit

In this paper, we describe an upgrade of the Alya code with up-to-date parallel linear solvers capable of achieving reliability, efficiency and scalability in the computation of the pressure field at each time step of the numerical procedure for solving a Large Eddy Simulation formulation of the incompressible Navier–Stokes equations. We developed a software module in the Alya’s kernel to interface the libraries included in the current version of PSCToolkit, a framework for the iterative solution of sparse linear systems, on parallel distributed-memory computers, by Krylov methods coupled to Algebraic MultiGrid preconditioners. The Toolkit has undergone various extensions within the EoCoE-II project with the primary goal of facing the exascale challenge. Results on a realistic benchmark for airflow simulations in wind farm applications show that the PSCToolkit solvers significantly outperform the original versions of the Conjugate Gradient method available in the Alya’s kernel in terms of scalability and parallel efficiency and represent a very promising software layer to move the Alya code toward exascale.

On the Numerical Solution of Sparse Linear Systems Emerging in Finite Volume Discretizations of 2D Boussinesq-Type Models on Unstructured Grids

This work aims to supplement the realization and validation of a higher-order well-balanced unstructured finite volume (FV) scheme, that has been relatively recently presented, for numerically simulating weakly non-linear weakly dispersive water waves over varying bathymetries. We investigate and develop solution strategies for the sparse linear system that appears during this FV discretisation of a set of extended Boussinesq-type equations on unstructured meshes. The resultant linear system of equations must be solved at each discrete time step as to recover the actual velocity field of the flow and advance in time. The system’s coefficient matrix is sparse, un-symmetric and often ill-conditioned. Its characteristics are affected by physical quantities of the problem to be solved, such as the undisturbed water depth and the mesh topology. To this end, we investigate the application of different well-known iterative techniques, with and without the usage of preconditioners and reordering, for the solution of this sparse linear system. The iiterative methods considered are the GMRES and the BiCGSTAB, three preconditioning techniques, including different ILU factorizations and two different reordering techniques are implemented and discussed. An optimal strategy, in terms of computational efficiency and robustness, is finally proposed which combines the use of the BiCGSTAB method with the ILUT preconditioner and the Reverse Cuthill–McKee reordering.

GPU-accelerated N-1 static security analysis based on fine-grained parallelism HELM

With the increment of scale of power system, the accuracy and efficiency requirements of N-1 static security analysis (SSA) keep increasing. In order to satisfy these requirements, this paper proposed a GPU-accelerated N-1 SSA method to reduce computation burden and ensure the accuracy of alternating current power flow (ACPF) and connectivity test simultaneously. First, in the algorithm, ACPF is performed by holomorphic embedding load flow method (HELM). The construction and solution of sparse linear system both adopts fine-grained parallelism. According to the mathematical lemma of combining linear equations, this paper performs fine-grained parallelism Padé approximation of HELM. Second, considering the high time complexity of the original search-based algorithm and adjacency matrix method which is not suitable for sparse graph, this paper proposes a GPU-accelerated Union Set algorithm. Its main advantage is lower time complexity compared to other graph theory methods. And most importantly, it is very suitable for parallelization. Shared memory can be used in this algorithm well. Therefore, it is not complicated to perform a block-level parallelism in GPU-accelerated Union Set. Finally, the overall framework of GPU-accelerated N-1 SSA method is shown in this paper. It has been tested on large-scale power systems of up to 2869 buses. The computation results are compared with MATPOWER and the execution times are compared with other mainstream method to show the improvement in performance.

Localised and shape-aware functions for spectral geometry processing and shape analysis: A survey & perspectives

Basis functions provide a simple and effective way to interpolate signals on a surface with an arbitrary dimension, topology, and discretisation. However, the definition of local and shape-aware basis functions is still an open problem, which has been addressed through optimisation methods that are time-consuming, over-constrained, or admit more than one solution. In this context, we review the definition, properties, computation, and applications of the class of Laplacian spectral basis functions, which are defined as solutions to the harmonic equation, the diffusion equation, and PDEs involving the Laplace–Beltrami operator. The resulting functions are efficiently computed through the solution of sparse linear systems and satisfy different properties, such as smoothness, locality, and multi-scale shape encoding. Finally, the discussion of the properties of the Laplacian spectral basis functions is integrated with an analysis of their behaviour with respect to different measures (i.e., conformal and area-preserving metrics, transformation distance) and of their applications to spectral geometry processing and shape analysis.

Accelerating advanced preconditioning methods on hybrid architectures

Many problems, in diverse areas of science and engineering, involve the solution of largescale sparse systems of linear equations. In most of these scenarios, they are also a computational bottleneck, and therefore their efficient solution on parallel architectureshas motivated a tremendous volume of research.This dissertation targets the use of GPUs to enhance the performance of the solution of sparse linear systems using iterative methods complemented with state-of-the-art preconditioned techniques. In particular, we study ILUPACK, a package for the solution of sparse linear systems via Krylov subspace methods that relies on a modern inverse-based multilevel ILU (incomplete LU) preconditioning technique.We present new data-parallel versions of the preconditioner and the most important solvers contained in the package that significantly improve its performance without affecting its accuracy. Additionally we enhance existing task-parallel versions of ILUPACK for shared- and distributed-memory systems with the inclusion of GPU acceleration. The results obtained show a sensible reduction in the runtime of the methods, as well as the possibility of addressing large-scale problems efficiently.

An a posteriori-based adaptive preconditioner for controlling a local algebraic error norm

This paper introduces an adaptive preconditioner for iterative solution of sparse linear systems arising from partial differential equations with self-adjoint operators. This preconditioner allows to control the growth rate of a dominant part of the algebraic error within a fixed point iteration scheme. Several numerical results that illustrate the efficiency of this adaptive preconditioner with a PCG solver are presented and the preconditioner is also compared with a previous variant in the literature.

Iteration-fusing conjugate gradient for sparse linear systems with MPI + OmpSs

In this paper, we target the parallel solution of sparse linear systems via iterative Krylov subspace-based method enhanced with a block-Jacobi preconditioner on a cluster of multicore processors. In order to tackle large-scale problems, we develop task-parallel implementations of the preconditioned conjugate gradient method that improve the interoperability between the message-passing interface and OmpSs programming models. Specifically, we progressively integrate several communication-reduction and iteration-fusing strategies into the initial code, obtaining more efficient versions of the method. For all these implementations, we analyze the communication patterns and perform a comparative analysis of their performance and scalability on a cluster consisting of 32 nodes with 24 cores each. The experimental analysis shows that the techniques described in the paper outperform the classical method by a margin that varies between 6 and 48%, depending on the evaluation.

Adaptive solution of linear systems of equations based on a posteriori error estimators

In this paper, we discuss a new adaptive approach for iterative solution of sparse linear systems arising from partial differential equations (PDEs) with self-adjoint operators. The idea is to use the a posteriori estimated local distribution of the algebraic error in order to steer and guide the solve process in such way that the algebraic error is reduced more efficiently in the consecutive iterations. We first explain the motivation behind the proposed procedure and show that it can be equivalently formulated as constructing a special combination of preconditioner and initial guess for the original system. We present several numerical experiments in order to identify when the adaptive procedure can be of practical use.

Accelerating the task/data-parallel version of ILUPACK’s BiCG in multi-CPU/GPU configurations

ILUPACK is a valuable tool for the solution of sparse linear systems via iterative Krylov subspace-based methods. Its relevance for the solution of real problems has motivated several efforts to enhance its performance on parallel machines. In this work we focus on exploiting the task-level parallelism derived from the structure of the BiCG method, in addition to the data-level parallelism of the internal matrix computations, with the goal of boosting the performance of a GPU (graphics processing unit) implementation of this solver. First, we revisit the use of dual-GPU systems to execute independent stages of the BiCG concurrently on both accelerators, while leveraging the extra memory space to improve the data access patterns. In addition, we extend our ideas to compute the BiCG method efficiently in multicore platforms with a single GPU. In this line, we study the possibilities offered by hybrid CPU-GPU computations, as well as a novel synchronization-free sparse triangular linear solver. The experimental results with the new solvers show important acceleration factors with respect to the previous data-parallel CPU and GPU versions.

Exact Solution of Sparse Linear Systems via Left-Looking Roundoff-Error-Free LU Factorization in Time Proportional to Arithmetic Work

The roundoff-error-free (REF) LU factorization, along with the REF forward and backward substitution algorithms, allows a rational system of linear equations to be solved exactly and efficiently. T...

An efficient GPU version of the preconditioned GMRES method

In a large number of scientific applications, the solution of sparse linear systems is the stage that concentrates most of the computational effort. This situation has motivated the study and development of several iterative solvers, among which preconditioned Krylov subspace methods occupy a place of privilege. In a previous effort, we developed a GPU-aware version of the GMRES method included in ILUPACK, a package of solvers distinguished by its inverse-based multilevel ILU preconditioner. In this work, we study the performance of our previous proposal and integrate several enhancements in order to mitigate its principal bottlenecks. The numerical evaluation shows that our novel proposal can reach important run-time reductions.

Adaptive precision in block‐Jacobi preconditioning for iterative sparse linear system solvers

SummaryWe propose an adaptive scheme to reduce communication overhead caused by data movement by selectively storing the diagonal blocks of a block‐Jacobi preconditioner in different precision formats (half, single, or double). This specialized preconditioner can then be combined with any Krylov subspace method for the solution of sparse linear systems to perform all arithmetic in double precision. We assess the effects of the adaptive precision preconditioner on the iteration count and data transfer cost of a preconditioned conjugate gradient solver. A preconditioned conjugate gradient method is, in general, a memory bandwidth‐bound algorithm, and therefore its execution time and energy consumption are largely dominated by the costs of accessing the problem's data in memory. Given this observation, we propose a model that quantifies the time and energy savings of our approach based on the assumption that these two costs depend linearly on the bit length of a floating point number. Furthermore, we use a number of test problems from the SuiteSparse matrix collection to estimate the potential benefits of the adaptive block‐Jacobi preconditioning scheme.

Parallel Implementation of the Coordinates-partitioning Based Aggregation-type Algebraic Multigrid Preconditioners

The coordinates-partitioning based aggregation-type algebraic multigrid preconditioners have been proven to be very efficient in the solution of sparse linear systems with conjugate gradient iterations. In this paper, a parallel algorithm for the setup is provided and the parallelization of the preconditioning process is also considered. The parallel algorithm is based on a good property of the coordinates-partitioning based aggregation, that is, the aggregation process is performed in a fashion from the coarsest level to the finest step by step. Thus, the original adjacent graph is partitioned into a number of sub-graphs, where each sub-graph is related to a node in the coarsest level and is assigned to a processor. The aggregation process can then proceed forward on each processor independently from the assigned sub-graph. When this kind of multigrid preconditioner is applied in Krylov subspace iterations, only the computation on the coarsest level and the matrix-vector multiplications related to the smoothing on each level require communication for V- and W-cycle versions, and require only some extra communication related to dot products for the K-cycle version. The size of the derived linear system on the coarsest level can be controlled by the adjustable arguments and this system can be solved again in parallel with some preconditioned Krylov subspace iterations. The structure information related to the matrix-vector multiplications is invariable to the iterations and can be derived in the setup and be stored, and is used in the latter iterations. Finally, the parallel algorithm is validated for the Gauss-Seidel smoother and some kinds of popular multigrid cycles in solving sparse linear systems from some two-dimensional model partial equations with preconditioned conjugate gradient iterations. The results show that the parallel efficiencies of both the setup and the iteration processes are satisfied.

SparseX

The Sparse Matrix-Vector Multiplication (SpMV) kernel ranks among the most important and thoroughly studied linear algebra operations, as it lies at the heart of many iterative methods for the solution of sparse linear systems, and often constitutes a severe performance bottleneck. Its optimization, which is intimately associated with the data structures used to store the sparse matrix, has always been of particular interest to the applied mathematics and computer science communities and has attracted further attention since the advent of multicore architectures. In this article, we present SparseX, an open source software package for SpMV targeting multicore platforms, that employs the state-of-the-art Compressed Sparse eXtended (CSX) sparse matrix storage format to deliver high efficiency through a highly usable “BLAS-like” interface that requires limited or no tuning. Performance results indicate that our library achieves superior performance over competitive libraries on large-scale problems.

Localized solutions of sparse linear systems for geometry processing

Computing solutions to linear systems is a fundamental building block of many geometry processing algorithms. In many cases the Cholesky factorization of the system matrix is computed to subsequently solve the system, possibly for many right-hand sides, using forward and back substitution. We demonstrate how to exploit sparsity in both the right-hand side and the set of desired solution values to obtain significant speedups. The method is easy to implement and potentially useful in any scenarios where linear problems have to be solved locally. We show that this technique is useful for geometry processing operations, in particular we consider the solution of diffusion problems. All problems profit significantly from sparse computations in terms of runtime, which we demonstrate by providing timings for a set of numerical experiments.

Communication in task‐parallel ILU‐preconditioned CG solvers using MPI + OmpSs

SummaryWe target the parallel solution of sparse linear systems via iterative Krylov subspace–based methods enhanced with incomplete LU (ILU)‐type preconditioners on clusters of multicore processors. In order to tackle large‐scale problems, we develop task‐parallel implementations of the classical iteration for the CG method, accelerated via ILUPACK and ILU(0) preconditioners, using MPI + OmpSs. In addition, we integrate several communication‐avoiding strategies into the codes, including the butterfly communication scheme and Eijkhout's formulation of the CG method. For all these implementations, we analyze the communication patterns and perform a comparative analysis of their performance and scalability on a cluster consisting of 16 nodes, with 16 cores each.

Sparse Matrix-Vector Multiplication on GPGPUs

The multiplication of a sparse matrix by a dense vector (SpMV) is a centerpiece of scientific computing applications: it is the essential kernel for the solution of sparse linear systems and sparse eigenvalue problems by iterative methods. The efficient implementation of the sparse matrix-vector multiplication is therefore crucial and has been the subject of an immense amount of research, with interest renewed with every major new trend in high-performance computing architectures. The introduction of General-Purpose Graphics Processing Units (GPGPUs) is no exception, and many articles have been devoted to this problem.With this article, we provide a review of the techniques for implementing the SpMV kernel on GPGPUs that have appeared in the literature of the last few years. We discuss the issues and tradeoffs that have been encountered by the various researchers, and a list of solutions, organized in categories according to common features. We also provide a performance comparison across different GPGPU models and on a set of test matrices coming from various application domains.

Variable-Size Batched Gauss-Huard for Block-Jacobi Preconditioning

In this work we present new kernels for the generation and application of block-Jacobi precon-ditioners that accelerate the iterative solution of sparse linear systems on graphics processing units (GPUs). Our approach departs from the conventional LU factorization and decomposes the diagonal blocks of the matrix using the Gauss-Huard method. When enhanced with column pivoting, this method is as stable as LU with partial/row pivoting. Due to extensive use of GPU registers and integration of implicit pivoting, our variable size batched Gauss-Huard implementation outperforms the batched version of LU factorization. In addition, the application kernel combines the conventional two-stage triangular solve procedure, consisting of a backward solve followed by a forward solve, into a single stage that performs both operations simultaneously.

Numerical recovery strategies for parallel resilient Krylov linear solvers

SummaryAs the computational power of high‐performance computing systems continues to increase by using a huge number of cores or specialized processing units, high‐performance computing applications are increasingly prone to faults. In this paper, we present a new class of numerical fault tolerance algorithms to cope with node crashes in parallel distributed environments. This new resilient scheme is designed at application level and does not require extra resources, that is, computational unit or computing time, when no fault occurs. In the framework of iterative methods for the solution of sparse linear systems, we present numerical algorithms to extract relevant information from available data after a fault, assuming a separate mechanism ensures the fault detection. After data extraction, a well‐chosen part of missing data is regenerated through interpolation strategies to constitute meaningful inputs to restart the iterative scheme. We have developed these methods, referred to as interpolation–restart techniques, for Krylov subspace linear solvers. After a fault, lost entries of the current iterate computed by the solver are interpolated to define a new initial guess to restart the Krylov method. A well‐suited initial guess is computed by using the entries of the faulty iterate available on surviving nodes. We present two interpolation policies that preserve key numerical properties of well‐known linear solvers, namely, the monotonic decrease of the A‐norm of the error of the conjugate gradient or the residual norm decrease of generalized minimal residual algorithm for solving. The qualitative numerical behavior of the resulting scheme has been validated with sequential simulations, when the number of faults and the amount of data losses are varied. Finally, the computational costs associated with the recovery mechanism have been evaluated through parallel experiments. Copyright © 2016 John Wiley & Sons, Ltd.

Efficient implementation of Jacobi iterative method for large sparse linear systems on graphic processing units

In this paper, an original Jacobi implementation is considered for the solution of sparse linear systems of equations. The proposed algorithm helps to optimize the parallel implementation on GPU. The performance analysis of GPU-based (using CUDA) algorithm of the implementation of this algorithm is compared to the corresponding serial CPU-based algorithm. Numerical experiments performed on a set of matrices arising from the finite element discretization of various equations (3D Laplace equation, 3D gravitational potential equation, 3D Heat equation) with different meshes, illustrate the performance, robustness and efficiency of our algorithm, with a speed up to 23$$\times $$× in double-precision arithmetics.