Sparse Matrix-vector Product Research Articles

Parallel implementation of an efficient preconditioned linear solver for grid-based applications in chemical physics. III: Improved parallel scalability for sparse matrix–vector products

The linear solve problems arising in chemical physics and many other fields involve large sparse matrices with a certain block structure, for which special block Jacobi preconditioners are found to be very efficient. In two previous papers [W. Chen, B. Poirier, Parallel implementation of efficient preconditioned linear solver for grid-based applications in chemical physics. I. Block Jacobi diagonalization, J. Comput. Phys. 219 (1) (2006) 185–197; W. Chen, B. Poirier, Parallel implementation of efficient preconditioned linear solver for grid-based applications in chemical physics. II. QMR linear solver, J. Comput. Phys. 219 (1) (2006) 198–209], a parallel implementation was presented. Excellent parallel scalability was observed for preconditioner construction, but not for the matrix–vector product itself. In this paper, we introduce a new algorithm with (1) greatly improved parallel scalability and (2) generalization for arbitrary number of nodes and data sizes.

Analyzing the execution of sparse matrix-vector product on the Finisterrae SMP-NUMA system

In this paper, the sparse matrix-vector product (SpMV) is evaluated on the FinisTerrae SMP-NUMA supercomputer. Its architecture particularities make the tuning of SpMV especially relevant due to the significant impact on the performance. First, we have estimated the influence of data and thread allocation. Moreover, because of the indirect and irregular memory access patterns of SpMV, we have also studied the influence of the memory hierarchy in the performance. According to the behavior observed in the study, a set of optimizations specially tuned for FinisTerrae were successfully applied to SpMV. Noticeable improvements are obtained in comparison with the SpMV naive implementation.

Parallel symmetric sparse matrix–vector product on scalar multi-core CPUs

We present a massively parallel implementation of symmetric sparse matrix–vector product for modern clusters with scalar multi-core CPUs. Matrices with highly variable structure and density arising from unstructured three-dimensional FEM discretizations of mechanical and diffusion problems are studied. A metric of the effective memory bandwidth is introduced to analyze the impact on performance of a set of simple, well-known optimizations: matrix reordering, manual prefetching, and blocking. A modification to the CRS storage improving the performance on multi-core Opterons is shown. The performance of an entire SMP blade rather than the per-core performance is optimized. Even for the simplest 4 node mechanical element our code utilizes close to 100% of the per-blade available memory bandwidth. We show that reducing the storage requirements for symmetric matrices results in roughly two times speedup. Blocking brings further storage savings and a proportional performance increase. Our results are compared to existing state-of-the-art implementations of SpMV, and to the dense BLAS2 performance. Parallel efficiency on 5400 Opteron cores of the Cray XT4 cluster is around 80–90% for problems with approximately 25 3 mesh nodes per core. For a problem with 820 million degrees of freedom the code runs with a sustained performance of 5.2 TeraFLOPs, over 20% of the theoretical peak.

A matrix approach to tomographic reconstruction and its implementation on GPUs

Electron tomography allows elucidation of the molecular architecture of complex biological specimens. Weighted backprojection (WBP) is the standard reconstruction method in the field. In this work, three-dimensional reconstruction with WBP is addressed from a matrix perspective by formulating the problem as a set of sparse matrix-vector products, with the matrix being constant and shared by all the products. This matrix approach allows efficient implementations of reconstruction algorithms. Although WBP is computationally simple, the resolution requirements may turn the tomographic reconstruction into a computationally intensive problem. Parallel systems have traditionally been used to cope with such demands. Recently, graphics processor units (GPUs) have emerged as powerful platforms for scientific computing and they are getting increasing interest. In combination with GPU computing, the matrix approach for WBP exhibits a significant acceleration factor compared to the standard implementation.

GPU ACCELERATED UNCONDITIONALLY STABLE CRANK-NICOLSON FDTD METHOD FOR THE ANALYSIS OF THREE-DIMENSIONAL MICROWAVE CIRCUITS

The programmable graphics processing unit (GPU) is employed to accelerate the unconditionally stable Crank-Nicolson flnite-difierence time-domain (CN-FDTD) method for the analysis of microwave circuits. In order to e-ciently solve the linear system from the CN-FDTD method at each time step, both the sparse matrix vector product (SMVP) and the arithmetic operations on vectors in the bi-conjugate gradient stabilized (Bi-CGSTAB) algorithm are performed with multiple processors of the GPU. Therefore, the GPU based BI-CGSTAB algorithm can signiflcantly speed up the CN-FDTD simulation due to parallel computing capability of modern GPUs. Numerical results demonstrate that this method is very efiective and a speedup factor of 10 can be achieved.

Increasing the Locality of Iterative Methods and Its Application to the Simulation of Semiconductor Devices

Irregular codes are present in many scientific applications, such as finite element simulations. In these simulations the solution of large sparse linear equation systems is required, which are often solved using iterative methods. The main kernel of the iterative methods is the sparse matrix—vector multiplication which frequently demands irregular data accesses. Therefore, techniques that increase the performance of this operation will have a great impact on the global performance of the iterative method and, as a consequence, on the simulations. In this paper a technique for improving the locality of sparse matrix codes is presented. The technique consists of reorganizing the data guided by a locality model instead of restructuring the code or changing the sparse matrix storage format. We have applied our proposal to different iterative methods provided by two standard numerical libraries. Results show an impact on the overall performance of the considered iterative method due to the increase in the locality of the sparse matrix—vector product. Noticeable reductions in the execution time have been achieved both in sequential and in parallel executions. This positive behavior allows the reordering technique to be successfully applied to real problems. We have focused on the simulation of semiconductor devices and in particular on the BIPS3D simulator. The technique was integrated into the simulator. Both sequential and parallel executions have been analyzed extensively in this paper. Noticeable reductions in the execution time required by the simulations are observed when using our reordered matrices in comparison with the original simulator.

A massively parallel exponential integrator for advection-diffusion models

This work considers the Real Leja Points Method (ReLPM), [M. Caliari, M. Vianello, L. Bergamaschi, Interpolating discrete advection-diffusion propagators at spectral Leja sequences, J. Comput. Appl. Math. 172 (2004) 79–99], for the exponential integration of large-scale sparse systems of ODEs, generated by Finite Element or Finite Difference discretizations of 3-D advection-diffusion models. We present an efficient parallel implementation of ReLPM for polynomial interpolation of the matrix exponential propagators exp ( Δ t A ) v and φ ( Δ t A ) v , φ ( z ) = ( exp ( z ) − 1 ) / z . A scalability analysis of the most important computational kernel inside the code, the parallel sparse matrix–vector product, has been performed, as well as an experimental study of the communication overhead. As a result of this study an optimized parallel sparse matrix–vector product routine has been implemented. The resulting code shows good scaling behavior even when using more than one thousand processors. The numerical results presented on a number of very large test cases gives experimental evidence that ReLPM is a reliable and efficient tool for the simulation of complex hydrodynamic processes on parallel architectures.

A 3-D Spectral-Element Time-Domain Method for Electromagnetic Simulation

A spectral-element time-domain (SETD) method is proposed to solve 3-D transient electromagnetic problems based on Gauss-Lobatto-Legendre polynomials. It has the advantages of spectral accuracy and block-diagonal mass matrix. With the inexpensive inversion of the block-diagonal mass matrix, the proposed method requires only a trivial sparse matrix-vector product at each time step, thus significantly reducing CPU time and memory requirement. Galerkin's method is used for spatial discretization, and a fourth-order Runge-Kutta scheme is employed for the time integration. The perfectly matched layer (PML) is employed to truncate the boundary in unbounded problems. The pseudospectral time-domain method is used to simplify the treatment of the PML inside the proposed SETD method. Numerical examples are shown to verify the efficiency and the spectral accuracy with the order of basis functions

Parallel implementation of efficient preconditioned linear solver for grid-based applications in chemical physics. II: QMR linear solver

Linear systems in chemical physics often involve matrices with a certain sparse block structure. These can often be solved very effectively using iterative methods (sequence of matrix–vector products) in conjunction with a block Jacobi preconditioner [B. Poirier, Numer. Linear Algebra Appl. 7 (2000) 715]. In a two-part series, we present an efficient parallel implementation, incorporating several additional refinements. The present study (paper II) indicates that the basic parallel sparse matrix–vector product operation itself is the overall scalability bottleneck, faring much more poorly than the specialized, block Jacobi routines considered in a companion paper (paper I). However, a simple dimensional combination scheme is found to alleviate this difficulty.

Optimizing Sparse Matrix–Vector Product Computations Using Unroll and Jam

Large-scale scientific applications frequently compute sparse matrix–vector products in their computational core. For this reason, techniques for computing sparse matrix– vector products efficiently on modern architectures are important. In this paper we describe a strategy for improving the performance of sparse matrix–vector product computations using a loop transformation known as unrollandjam. We describe a novel sparse matrix representation that enables us to apply this transformation. Our approach is best suited for sparse matrices that have rows with a small number of predictable lengths. This work was motivated by sparse matrices that arise in SAGE, an application from Los Alamos National Laboratory. We evaluate the performance benefits of our approach using sparse matrices produced by SAGE for a pair of sample inputs. We show that our strategy is effective for improving sparse matrix–vector product performance using these matrices on MIPS R12000, Alpha Ev67, IBM Power 3, and Itanium 2 processors. Our measurements show that for this class of sparse matrices, our strategy improves sparse matrix–vector product performance from a low of 41% on MIPS to well over a factor of 2 on Itanium.

Implicit SUPG solution of Euler equations using edge-based data structures

In this work we present an implicit, edge-based implementation of the semi-discrete SUPG formulation with shock-capturing for the Euler equations in conservative variables. By disassembling the resulting finite element matrices into their edge contributions, sparse matrix coefficients, residuals and matrix-vector products needed in Krylov-update techniques are computed based on edge data structures. The resulting solution method requires less memory and CPU time than element-based implementations.

Modeling and improving locality for the sparse-matrix–vector product on cache memories

A model for representing and improving the locality exhibited by the execution of sparse irregular problems is developed in this work. We focus on the product of a sparse matrix by a dense vector ( SpM× V). We consider the cache memory as a representative level of the memory hierarchy. Locality is evaluated through four functions based on two parameters called entry matches and line matches. In order to increase the locality, two algorithms are applied: one based on the construction of minimum spanning trees and the other on the nearest-neighbor heuristic. These techniques were tested and compared with some standard ordering algorithms.

Towards a fast parallel sparse symmetric matrix–vector multiplication

The sparse matrix–vector product is an important computational kernel that runs ineffectively on many computers with super-scalar RISC processors. In this paper we analyse the performance of the sparse matrix–vector product with symmetric matrices originating from the FEM and describe techniques that lead to a fast implementation. It is shown how these optimisations can be incorporated into an efficient parallel implementation using message-passing. We conduct numerical experiments on many different machines and show that our optimisations speed up the sparse matrix–vector multiplication substantially.

Modeling data locality for the sparse matrix–vector product using distance measures

In this work, we model the data locality in the execution of codes with irregular accesses. We focus on the product of a sparse matrix by a dense vector ( SpM× V). In the model, locality is established taking into account pairs of rows or columns of sparse matrices. In order to evaluate this locality three functions are introduced based on two parameters: number of entry matches and number of block matches. The model is generalized considering windows of locality (groups of consecutive rows/columns of the matrix). We show results for a broad set of matrices measuring the goodness of our predictions of locality.

Edge‐based finite element techniques for non‐linear solid mechanics problems

AbstractEdge‐based data structures are used to improve computational efficiency of inexact Newton methods for solving finite element non‐linear solid mechanics problems on unstructured meshes. Edge‐based data structures are employed to store the stiffness matrix coefficients and to compute sparse matrix–vector products needed in the inner iterative driver of the inexact Newton method. Numerical experiments on three‐dimensional plasticity problems have shown that memory and computer time are reduced, respectively, by factors of 4 and 6, compared with solutions using element‐by‐element storage and matrix–vector products. Copyright © 2001 John Wiley & Sons, Ltd.

Modeling set associative caches behavior for irregular computations

While much work has been devoted to the study of cache behavior during the execution of codes with regular access patterns, little attention has been paid to irregular codes. An important portion of these codes are scientific applications that handle compressed sparse matrices. In this work a probabilistic model for the prediction of the number of misses on a K -way associative cache memory considering sparse matrices with a uniform or banded distribution is presented. Two different irregular kernels are considered: the sparse matrix-vector product and the transposition of a sparse matrix. The model was validated with simulations on synthetic uniform matrices and banded matrices from the Harwell-Boeing collection.

Implementation of a Portable Nested Data-Parallel Language

This paper gives an overview of the implementation of NESL, a portable nested data-parallel language. This language and its implementation are the first to fully support nested data structures as well as nested data-parallel function calls. These features allow the concise description of parallel algorithms on irregular data, such as sparse matrices and graphs. In addition, they maintain the advantages of data-parallel languages: a simple programming model and portability. The current NESL implementation is based on an intermediate language called VCODE and a library of vector routines called CVL. It runs on the Connection Machines CM-2 and CM-5, the Cray Y-MP C90, and serial workstations. We compare initial benchmark results of NESL, with those of machine-specific code on these machines for three algorithms: least-squares line-fitting, median finding, and a sparse-matrix vector product. These results show that NESL′s performance is competitive with that of machine-specific codes for regular dense data, and is often superior for irregular data.

Implementation of a portable nested data-parallel language

This paper gives an overview of the implementation of NESL, a portable nested data-parallel language. This language and its implementation are the first to fully support nested data structures as well as nested data-parallel function calls. These features allow the concise description of parallel algorithms on irregular data, such as sparse matrices and graphs. In addition, they maintain the advantages of data-parallel languages: a simple programming model and portability. The current NESL implementation is based on an intermediate language called VCODE and a library of vector routines called CVL. It runs on the Connection Machine CM-2, the Cray Y-MP C90, and serial machines. We compare initial benchmark results of NESL with those of machine-specific code on these machines for three algorithms: least-squares line-fitting, median finding, and a sparse-matrix vector product. These results show that NESL's performance is competitive with that of machine-specific codes for regular dense data, and is often superior for irregular data.

EFFICIENT ICCG ON A SHARED MEMORY MULTIPROCESSOR

Different approaches are discussed for exploiting parallelism in the ICCG (Incomplete Cholesky Conjugate Gradient) method for solving large sparse symmetric positive definite systems of equations on a shared memory parallel computer. Techniques for efficiently solving triangular systems and computing sparse matrix-vector products are explored. Three methods for scheduling the tasks in solving triangular systems are implemented on the Sequent Balance 21000. Sample problems that are representative of a large class of problems solved using iterative methods are used. We show that a static analysis to determine data dependences in the triangular solve can greatly improve its parallel efficiency. We also show that ignoring symmetry and storing the whole matrix can reduce solution time substantially.

A new storage scheme for an efficient implementation of the sparse matrix-vector product

The most critical point in the implementation on a vector computer of iterative methods to solve sparse linear algebraic systems is the sparse matrix-vector product. Therefor we have compared its performance under three different sparse matrix storage schemes found in literature. These schemes have been tested on a FPS M64/330 using matrices similar to those arising in finite element or finite difference discretizations. The results are reported and discussed. On the basis of them a new storage scheme is proposed and tested in similar situations. The obtained performance is never worse and often much better than those of the other schemes we have analyzed.