Sparse Matrix-vector Product Research Articles

Improving performance of iterative solvers with the AXC format using the Intel Xeon Phi

This work is focused on the application of the new AXC format in iterative algorithms on the Intel Xeon Phi coprocessor to solve linear systems by accelerating the sparse matrix-vector (SpMV) product. These algorithms are the Conjugate Gradient (CG) and the Biconjugate Gradient Stabilised (BiCGS) methods used to solve symmetric and non-symmetric linear systems. Two highly efficient formats were selected to compare the AXC format: the Compressed Sparse Row (CSR) format and the Sliced ELLPACK-C- $$\alpha $$ (SELL-C- $$\alpha $$ ) format. The evaluation process consists of two phases. The first phase is focused on the performance comparison of the SpMV kernels alone. The second phase compares the performance of the solvers with the different kernels integrated within them. As this work is oriented to explore the full capabilities of the Intel Xeon Phi architecture, the Intel’s intrinsic instructions set were used to code the kernels for the AXC and SELL-C- $$\alpha $$ formats and compared with the Intel’s MKL optimal implementation of the CSR format. Numerical results demonstrate that the AXC format achieves the higher average performance 15.9 GFLOPS against the 5.5 GFLOPS for the SELL-C- $$\alpha $$ format and 8.4 GFLOPs for the MKL CSR format, and it also have the lowest variability when performing the SpMV product alone. Results also show that the AXC format achieves up to 6.5 $$\times $$ , 1.5 $$\times $$ and 1.9 $$\times $$ greater performance over the MKL CSR format, its closest competitor, for the SpMV product, the CG and the BiCGS algorithms, respectively.

Extension of the parallel Sparse Matrix Vector Product (SpMV) for the implicit coupling of PDEs on non-matching meshes

The Sparse Matrix Vector Product (SpMV) is one of the main operations of iterative solvers, and, in a parallel context, it is also the siege of point-to-point communications between the neighboring MPI processes. The parallel SpMV is built in such a way that it gives, up to round off errors, the same result as its sequential counterpart. In this regards, nodes on the interfaces (or halo nodes if halo are considered) are duplicated nodes of the same original mesh. It is therefore limited to matching meshes. In this work, we generalize the parallel SpMV to glue the solution of non-matching (non-conforming) meshes through the introduction of transmission matrices. This extension of the SpMV thus enables the implicit and parallel solution of partial differential equations on non-matching meshes, as well as the implicit coupling of multiphysics problems, such as fluid-structure interactions. The proposed method is developed similarly to classical parallelization techniques and can therefore be implemented by modifying few subroutines of an already MPI-based code. According to the proposed framework, the classical parallelization technique appears as a particular case of this general setup.

HPC2—A fully-portable, algebra-based framework for heterogeneous computing. Application to CFD

The variety of computing architectures competing in the exascale race makes the portability of codes of major importance. In this work, the HPC2 code is presented as a fully-portable, algebra-based framework suitable for heterogeneous computing. In its application to CFD, the algorithm of the time-integration phase relies on a reduced set of only three algebraic operations: the sparse matrix-vector product, the linear combination of vectors and the dot product. This algebraic approach combined with a multilevel MPI+OpenMP+OpenCL parallelization naturally provides portability. The performance has been studied on different architectures including multicore CPUs, Intel Xeon Phi accelerators and GPUs of AMD and NVIDIA. The multi-GPU scalability is demonstrated up to 256 devices. The heterogeneous execution is tested on a CPU+GPU hybrid cluster. Finally, results of the direct numerical simulation of a turbulent flow in a 3D air-filled differentially heated cavity are presented to show the capabilities of the HPC2 dealing with large-scale CFD simulations.

A Performance Comparison of Continuous and Discontinuous Galerkin Methods with Fast Multigrid Solvers

This study presents a fair performance comparison of the continuous finite element method, the symmetric interior penalty discontinuous Galerkin method, and the hybridized discontinuous Galerkin method. Modern implementations of high-order methods with state-of-the-art multigrid solvers for the Poisson equation are considered, including fast matrix-free implementations with sum factorization on quadrilateral and hexahedral elements. For the hybridized discontinuous Galerkin method, a multigrid approach that combines a grid transfer from the trace space to the space of linear finite elements with algebraic multigrid on further levels is developed. Despite similar solver complexity of the matrix-based HDG solver and matrix-free geometric multigrid schemes with continuous and discontinuous Galerkin finite elements, the latter offer up to order of magnitude faster time to solution, even after including the superconvergence effects. This difference is because of vastly better performance of matrix-free operator evaluation as compared to sparse matrix-vector products. A roofline performance model confirms the advantage of the matrix-free implementation.

Portable implementation model for CFD simulations. Application to hybrid CPU/GPU supercomputers

ABSTRACTNowadays, high performance computing (HPC) systems experience a disruptive moment with a variety of novel architectures and frameworks, without any clarity of which one is going to prevail. In this context, the portability of codes across different architectures is of major importance. This paper presents a portable implementation model based on an algebraic operational approach for direct numerical simulation (DNS) and large eddy simulation (LES) of incompressible turbulent flows using unstructured hybrid meshes. The strategy proposed consists in representing the whole time-integration algorithm using only three basic algebraic operations: sparse matrix–vector product, a linear combination of vectors and dot product. The main idea is based on decomposing the nonlinear operators into a concatenation of two SpMV operations. This provides high modularity and portability. An exhaustive analysis of the proposed implementation for hybrid CPU/GPU supercomputers has been conducted with tests using up to 128 GPUs. The main objective consists in understanding the challenges of implementing CFD codes on new architectures.

Evaluation of the transport matrix method for simulation of ocean biogeochemical tracers

Abstract. Conventional integration of Earth system and ocean models can accrue considerable computational expenses, particularly for marine biogeochemical applications. Offline numerical schemes in which only the biogeochemical tracers are time stepped and transported using a pre-computed circulation field can substantially reduce the burden and are thus an attractive alternative. One such scheme is the transport matrix method (TMM), which represents tracer transport as a sequence of sparse matrix–vector products that can be performed efficiently on distributed-memory computers. While the TMM has been used for a variety of geochemical and biogeochemical studies, to date the resulting solutions have not been comprehensively assessed against their online counterparts. Here, we present a detailed comparison of the two. It is based on simulations of the state-of-the-art biogeochemical sub-model embedded within the widely used coarse-resolution University of Victoria Earth System Climate Model (UVic ESCM). The default, non-linear advection scheme was first replaced with a linear, third-order upwind-biased advection scheme to satisfy the linearity requirement of the TMM. Transport matrices were extracted from an equilibrium run of the physical model and subsequently used to integrate the biogeochemical model offline to equilibrium. The identical biogeochemical model was also run online. Our simulations show that offline integration introduces some bias to biogeochemical quantities through the omission of the polar filtering used in UVic ESCM and in the offline application of time-dependent forcing fields, with high latitudes showing the largest differences with respect to the online model. Differences in other regions and in the seasonality of nutrients and phytoplankton distributions are found to be relatively minor, giving confidence that the TMM is a reliable tool for offline integration of complex biogeochemical models. Moreover, while UVic ESCM is a serial code, the TMM can be run on a parallel machine with no change to the underlying biogeochemical code, thus providing orders of magnitude speed-up over the online model.

Fast Finite Element Analysis Method Using Multiple GPUs for Crustal Deformation and its Application to Stochastic Inversion Analysis with Geometry Uncertainty

Crustal deformation computation using 3-D high-fidelity models has been in heavy demand due to accumulation of observational data. This approach is computationally expensive and more than 105 repetitive computations are required for various application including Monte Carlo simulation, stochastic inverse analysis, and optimization. To handle the massive computation cost, we develop a fast Finite Element (FE) analysis method using multi-GPUs for crustal deformation. We use algorithms appropriate for GPUs and accelerate calculations such as sparse matrix-vector product. By reducing the computation time, we are able to conduct multiple crustal deformation computations in a feasible timeframe. As an application example, we conduct stochastic inverse analysis considering uncertainties in geometry and estimate coseismic slip distribution in the 2011 Tohoku Earthquake, by performing 360,000 crustal deformation computations for different 8 × 107 DOF FE models using the proposed method.

On the performance and energy efficiency of sparse linear algebra on GPUs

In this paper we unveil some performance and energy efficiency frontiers for sparse computations on GPU-based supercomputers. We compare the resource efficiency of different sparse matrix–vector products (SpMV) taken from libraries such as cuSPARSE and MAGMA for GPU and Intel’s MKL for multicore CPUs, and develop a GPU sparse matrix–matrix product (SpMM) implementation that handles the simultaneous multiplication of a sparse matrix with a set of vectors in block-wise fashion. While a typical sparse computation such as the SpMV reaches only a fraction of the peak of current GPUs, we show that the SpMM succeeds in exceeding the memory-bound limitations of the SpMV. We integrate this kernel into a GPU-accelerated Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) eigensolver. LOBPCG is chosen as a benchmark algorithm for this study as it combines an interesting mix of sparse and dense linear algebra operations that is typical for complex simulation applications, and allows for hardware-aware optimizations. In a detailed analysis we compare the performance and energy efficiency against a multi-threaded CPU counterpart. The reported performance and energy efficiency results are indicative of sparse computations on supercomputers.

Merge-based sparse matrix-vector multiplication (SpMV) using the CSR storage format

We present a perfectly balanced, "merge-based" parallel method for computing sparse matrix-vector products (SpMV). Our algorithm operates directly upon the Compressed Sparse Row (CSR) sparse matrix format, a predominant in-memory representation for general-purpose sparse linear algebra computations. Our CsrMV performs an equitable multi-partitioning of the input dataset, ensuring that no single thread can be overwhelmed by assignment to (a) arbitrarily-long rows or (b) an arbitrarily-large number of zero-length rows. This parallel decomposition requires neither offline preprocessing nor specialized/ancillary data formats. We evaluate our method on both CPU and GPU microarchitecture across an enormous corpus of diverse real world matrix datasets. We show that traditional CsrMV methods are inconsistent performers subject to order-of-magnitude slowdowns, whereas the performance response of our method is substantially impervious to row-length heterogeneity.

ViennaCL---Linear Algebra Library for Multi- and Many-Core Architectures

CUDA, OpenCL, and OpenMP are popular programming models for the multicore architectures of CPUs and many-core architectures of GPUs or Xeon Phis. At the same time, computational scientists face the question of which programming model to use to obtain their scientific results. We present the linear algebra library ViennaCL, which is built on top of all three programming models, thus enabling computational scientists to interface to a single library, yet obtain high performance for all three hardware types. Since the respective compute back end can be selected at runtime, one can seamlessly switch between different hardware types without the need for error-prone and time-consuming recompilation steps. We present new benchmark results for sparse linear algebra operations in ViennaCL, complementing results for the dense linear algebra operations in ViennaCL reported in earlier work. Comparisons with vendor libraries show that ViennaCL provides better overall performance for sparse matrix-vector and sparse matrix-matrix products. Additional benchmark results for pipelined iterative solvers with kernel fusion and preconditioners identify the respective sweet spots for CPUs, Xeon Phis, and GPUs.

Time-domain BEM for the wave equation on distributed-heterogeneous architectures: A blocking approach

We parallelize the TD-BEM on clusters of multicore nodes enhanced with multiple GPU.SpMV is replaced by a more efficient operator by reordering the computation.GPU kernels achieve high Flop-rate using blocking schemes.Idle time is drastically reduced by using a greedy balancing algorithm. The problem of time-domain BEM for the wave equation in acoustics and electromagnetism can be expressed as a sparse linear system composed of multiple interaction/convolution matrices. It can be solved by using sparse matrix-vector products which are inefficient to achieve high Flop-rate neither on CPUs nor GPUs. In this paper we extend the approach proposed in a previous work 1 in which we re-order the computation to get a special matrix structure with one dense vector per row. This new structure is called a slice matrix and is computed with a custom matrix/vector product operator. In this study, we present an optimized implementation of this operator on Nvidia GPUs based on two blocking strategies. We explain how we can obtain multiple block-values from a slice and how these can be computed efficiently on GPUs since we target heterogeneous nodes composed of CPUs and GPUs. In order to deal with different efficiencies of the processing units we use a greedy heuristic that dynamically balances work among the workers. We demonstrate the performance of our system by studying the quality of the balancing heuristic and the sequential Flop-rate of the blocked implementations. Finally, we validate our implementation with an industrial test case on 8 heterogeneous nodes, each composed of 12 CPUs and 3 GPUs.

Acceleration of GPU-based Krylov solvers via data transfer reduction

Krylov subspace iterative solvers are often the method of choice when solving large sparse linear systems. At the same time, hardware accelerators such as graphics processing units continue to offer significant floating point performance gains for matrix and vector computations through easy-to-use libraries of computational kernels. However, as these libraries are usually composed of a well optimized but limited set of linear algebra operations, applications that use them often fail to reduce certain data communications, and hence fail to leverage the full potential of the accelerator. In this paper, we target the acceleration of Krylov subspace iterative methods for graphics processing units, and in particular the Biconjugate Gradient Stabilized solver that significant improvement can be achieved by reformulating the method to reduce data-communications through application-specific kernels instead of using the generic BLAS kernels, e.g. as provided by NVIDIA’s cuBLAS library, and by designing a graphics processing unit specific sparse matrix-vector product kernel that is able to more efficiently use the graphics processing unit’s computing power. Furthermore, we derive a model estimating the performance improvement, and use experimental data to validate the expected runtime savings. Considering that the derived implementation achieves significantly higher performance, we assert that similar optimizations addressing algorithm structure, as well as sparse matrix-vector, are crucial for the subsequent development of high-performance graphics processing units accelerated Krylov subspace iterative methods.

Optimized Schwarz method without overlap for the gravitational potential equation on cluster of graphics processing unit

Many engineering and scientific problems need to solve boundary value problems for partial differential equations or systems of them. For most cases, to obtain the solution with desired precision and in acceptable time, the only practical way is to harness the power of parallel processing. In this paper, we present some effective applications of parallel processing based on hybrid CPU/GPU domain decomposition method. Within the family of domain decomposition methods, the so-called optimized Schwarz methods have proven to have good convergence behaviour compared to classical Schwarz methods. The price for this feature is the need to transfer more physical information between subdomain interfaces. For solving large systems of linear algebraic equations resulting from the finite element discretization of the subproblem for each subdomain, Krylov method is often a good choice. Since the overall efficiency of such methods depends on effective calculation of sparse matrix–vector product, approaches that use graphics processing unit (GPU) instead of central processing unit (CPU) for such task look very promising. In this paper, we discuss effective implementation of algebraic operations for iterative Krylov methods on GPU. In order to ensure good performance for the non-overlapping Schwarz method, we propose to use optimized conditions obtained by a stochastic technique based on the covariance matrix adaptation evolution strategy. The performance, robustness, and accuracy of the proposed approach are demonstrated for the solution of the gravitational potential equation for the data acquired from the geological survey of Chicxulub crater.

A Task-Scheduling Approach for Efficient Sparse Symmetric Matrix-Vector Multiplication on a GPU

In this paper, a task-scheduling approach to efficiently calculating sparse symmetric matrix-vector products and designed to run on graphics processing units (GPUs) is presented. The main premise is that, for many sparse symmetric matrices occurring in common applications, it is possible to obtain significant reductions in memory usage and improvements in performance when the matrix is prepared in certain ways prior to computation. The preprocessing proposed in this paper employs task scheduling to overcome the difficulties that have suppressed the development of methods taking advantage of the symmetry of sparse matrices. The performance of the proposed task-scheduling method is verified using a Kepler (Tesla K40c) graphics accelerator, and is compared to the performance of cuSPARSE library functions on a GPU and to functions from the Intel MKL on central processing units (CPUs) executed in the parallel mode. The obtained results indicate that the proposed approach for sparse symmetric matrix-vector products results in up to a 40% reduction in memory usage, as compared to nonsymmetric matrix storage formats, while retaining good throughput. Compared to cuSPARSE and Intel MKL functions for sparse symmetric matrices, the proposed TSMV approach allowed us to achieve a significant speedup (of over one order of magnitude).

GPU accelerated sparse matrix‐vector multiplication and sparse matrix‐transpose vector multiplication

SummaryMany high performance computing applications require computing both sparse matrix‐vector product (SMVP) and sparse matrix‐transpose vector product (SMTVP) for better overall performance. Under such a circumstance, it is critical to maintain a similarly high throughput for these two computing patterns with the underlying sparse matrix encoded in a single storage format. The compressed sparse block (CSB) format proposed by Buluç et al. allows computing both problems on multi‐core CPUs with nearly identical throughputs. On the other hand, a direct porting of CSB to graphics processing units (GPUs), which have been recently recognized as a powerful general purpose computing platform, turns out to be inefficient. In this work, we propose a new data structure, designated as expanded CSB (eCSB), to minimize the throughput gap between SMVP and SMTVP computations on GPUs, while at the same time enable a high computing throughput. We also use a hybrid storage format to store elements in each block, which can be selected dynamically at runtime. Experimental results show that the proposed techniques implemented on a Kepler GPU delivers similar throughput on both SMVP and SMTVP and the throughput is up to 13 times faster than that of the CPU‐based CSB implementation. In addition, our eCSB procedure outperforms the previous GPU results by up to 188% and 914% in computing SMVP and SMTVP, and we validate the effectiveness of eCSB by means of wall‐clock time of bi‐conjugate gradient algorithm; our eCSB is 25% faster than Compressed Sparse Rows (CSR) and 6% faster than HYB, respectively. Copyright © 2014 John Wiley & Sons, Ltd.

Parallel resolution of the 3D Helmholtz equation based on multi‐graphics processing unit clusters

SummaryThe resolution of the 3D Helmholtz equation is required in the development of models related to a wide range of scientific and technological applications. For solving this equation in complex arithmetic, the biconjugate gradient (BCG) method is one of the most relevant solvers. However, this iterative method has a high computational cost because of the large sparse matrix and the vector operations involved. In this paper, a specific BCG method, adapted for the regularities of the Helmholtz equation is presented. This BCG is based on the implementation of a novel format (named ‘Regular Format’) that allows the storage of the large sparse matrix involved in the sparse matrix vector product in a compact form. The contribution of this work is twofold: (1) decreasing the memory requirements of the 3D Helmholtz equation using the ‘Regular Format’ and (2) speeding up the resolution of the equation using high performance computing resources. A hybrid Message Passing Interface (MPI)‐graphics processing unit CUDA GPU parallelization that is capable of solving complex problems in short time has carried out (Fast‐Helmholtz). Fast‐Helmholtz combines optimizations at Message Passing Interface and GPU levels to reduce communications costs and to improve the exploitation of GPU architecture. This strategy makes it possible to extend the dimension of the Helmholtz problem to be solved, thanks to the relevant reduction of memory requirements and runtime. Copyright © 2014 John Wiley & Sons, Ltd.

Compressed Multirow Storage Format for Sparse Matrices on Graphics Processing Units

A new format for storing sparse matrices is proposed for efficient sparse matrix-vector (SpMV) product calculation on modern graphics processing units (GPUs). This format extends the standard compressed row storage (CRS) format and can be quickly converted to and from it. Computational performance of two SpMV kernels for the new format is determined for over 130 sparse matrices on Fermi-class and Kepler-class GPUs and compared with that of five existing generic algorithms and industrial implementations, including Nvidia cuSparse CSR and HYB kernels. We found the speedup of up to $\approx 60%$ over the best of the five alternative kernels.

An efficient GPU acceleration format for finite element analysis

This paper proposes a new Graphics Processing Unit (GPU)-accelerated storage format to speed up Sparse Matrix Vector Products (SMVPs) for Finite Element Method (FEM) analysis of electromagnetic problems. A new format called Modified Compile Time Optimization (MCTO) format is used to reduce much execution time and design for hastening the iterative solution of FEM equations especially when rows have uneven lengths. The MCTO-applied FEM is about 10 times faster than conventional FEM on a CPU, and faster than other row-major ordering formats on a GPU. Numerical results show that the proposed GPU-accelerated storage format turns out to be an excellent accelerator.

A hardware counter‐based toolkit for the analysis of memory accesses in SMPs

SUMMARYIn this paper, a set of three hardware counter (HC)‐based tools to characterise memory access of parallel codes in Symmetric Multiprocessors (SMPs) is presented. This toolkit simplifies accessing and programming HCs, which are included in modern microprocessors. Hardware counters are used to obtain information about memory accesses in a parallel code at very low cost. This information is presented to the user in a friendly way. The first tool can be used to automatically monitor the memory accesses of a system and to analyse a code even if the source is not available. The second tool allows the user to insert in a source code, in a simple and transparent way, the instructions needed to monitor and manage HCs. This way, specific parts of the code can be analysed. The user can either add appropriate directives to a C code or use a graphical interface to select those parts of the code to be analysed. The tool takes this source file and automatically adds the monitoring code. The third tool takes the information gathered by the aforementioned tools, processes it and displays it graphically. This tool shows the information in a comprehensive and simple way, allowing the user to adjust the level of detail. The aim of these tools was to characterise the memory accesses of parallel codes in multicore systems, in which the cache hierarchy can greatly influence the performance. For illustrative purposes, these tools were used to carry out two case studies, a sparse matrix vector product and a dot product. These studies have been made in two different environments. Anyway, they can be used in almost any system as long as the necessary HCs are available.Copyright © 2013 John Wiley & Sons, Ltd.

FastSpMM: An Efficient Library for Sparse Matrix Matrix Product on GPUs

Sparse matrix matrix (SpMM) multiplication is involved in a wide range of scientific and technical applications. The computational requirements for this kind of operation are enormous, especially for large matrices. This paper analyzes and evaluates a method to efficiently compute the SpMM product in a computing environment that includes graphics processing units (GPUs). Some libraries to compute this matricial operation can be found in the literature. However, our strategy (FastSpMM) outperforms the existing approaches because it combines the use of the ELLPACK-R storage format with the exploitation of the high ratio computation/memory access of the SpMM operation and the overlapping of CPU–GPU communications/computations by Compute Unified Device Architecture streaming computation. In this work, FastSpMM is described and its performance evaluated with regard to the CUSPARSE library (supplied by NVIDIA), which also includes routines to compute SpMM on GPUs. Experimental evaluations based on a representative set of test matrices show that, in terms of performance, FastSpMM outperforms the CUSPARSE routine as well as the implementation of the SpMM as a set of sparse matrix vector products.