Sparse Matrix Storage Schemes Research Articles

A new AXT format for an efficient SpMV product using AVX-512 instructions and CUDA

The Sparse Matrix-Vector (SpMV) product is a key operation used in many scientific applications. This work proposes a new sparse matrix storage scheme, the AXT format, that improves the SpMV performance on vector capability platforms. AXT can be adapted to different platforms, improving the storage efficiency for matrices with different sparsity patterns. Intel AVX-512 instructions and CUDA are used to optimise the performances of the four different AXT subvariants. Performance comparisons are made with the Compressed Sparse Row (CSR) and AXC formats on an Intel Xeon Gold 6148 processor and an NVIDIA Tesla V100 Graphics Processing Units using 26 matrices. On the Intel platform the overall AXT performance is 18% and 44.3% higher than the AXC and CSR respectively, reaching speed-up factors of up to x7.33. On the NVIDIA platform the AXT performance is 44% and 8% higher than the AXC and CSR performances respectively, reaching speed-up factors of up to x378.5.

Performance Analysis of Sparse Matrix-Vector Multiplication (SpMV) on Graphics Processing Units (GPUs)

Graphics processing units (GPUs) have delivered a remarkable performance for a variety of high performance computing (HPC) applications through massive parallelism. One such application is sparse matrix-vector (SpMV) computations, which is central to many scientific, engineering, and other applications including machine learning. No single SpMV storage or computation scheme provides consistent and sufficiently high performance for all matrices due to their varying sparsity patterns. An extensive literature review reveals that the performance of SpMV techniques on GPUs has not been studied in sufficient detail. In this paper, we provide a detailed performance analysis of SpMV performance on GPUs using four notable sparse matrix storage schemes (compressed sparse row (CSR), ELLAPCK (ELL), hybrid ELL/COO (HYB), and compressed sparse row 5 (CSR5)), five performance metrics (execution time, giga floating point operations per second (GFLOPS), achieved occupancy, instructions per warp, and warp execution efficiency), five matrix sparsity features (nnz, anpr, nprvariance, maxnpr, and distavg), and 17 sparse matrices from 10 application domains (chemical simulations, computational fluid dynamics (CFD), electromagnetics, linear programming, economics, etc.). Subsequently, based on the deeper insights gained through the detailed performance analysis, we propose a technique called the heterogeneous CPU–GPU Hybrid (HCGHYB) scheme. It utilizes both the CPU and GPU in parallel and provides better performance over the HYB format by an average speedup of 1.7x. Heterogeneous computing is an important direction for SpMV and other application areas. Moreover, to the best of our knowledge, this is the first work where the SpMV performance on GPUs has been discussed in such depth. We believe that this work on SpMV performance analysis and the heterogeneous scheme will open up many new directions and improvements for the SpMV computing field in the future.

Efficiency analysis of numerical integrations for finite element substructure in real-time hybrid simulation

Finite element (FE) is a powerful tool and has been applied by investigators to real-time hybrid simulations (RTHSs). This study focuses on the computational efficiency, including the computational time and accuracy, of numerical integrations in solving FE numerical substructure in RTHSs. First, sparse matrix storage schemes are adopted to decrease the computational time of FE numerical substructure. In this way, the task execution time (TET) decreases such that the scale of the numerical substructure model increases. Subsequently, several commonly used explicit numerical integration algorithms, including the central difference method (CDM), the Newmark explicit method, the Chang method and the Gui-λ method, are comprehensively compared to evaluate their computational time in solving FE numerical substructure. CDM is better than the other explicit integration algorithms when the damping matrix is diagonal, while the Gui-λ (λ = 4) method is advantageous when the damping matrix is non-diagonal. Finally, the effect of time delay on the computational accuracy of RTHSs is investigated by simulating structure-foundation systems. Simulation results show that the influences of time delay on the displacement response become obvious with the mass ratio increasing, and delay compensation methods may reduce the relative error of the displacement peak value to less than 5% even under the large time-step and large time delay.

Numerical modeling of borehole-surface electromagnetic responses with 3-D finite difference method and comparison with physical simulations

Borehole-surface electromagnetic method (BSEM) is a high-accuracy electromagnetic prospecting method that uses AC-powered vertical finite line source as the excitation source. The observed surface electromagnetic field is inverted for the subsurface conductivity structure. In this paper, the total field is separated into a primary part, due to a horizontally layered host medium, and a secondary part due to 3-D heterogeneities. After solving for the primary field with an analytical method, the frequency-domain second-order differential equation for the secondary field is discretized by finite differences. A sparse matrix storage scheme is employed and a BiCGSSTAB(m) method with a diagonal matrix preconditioner is used to obtain the secondary field as well as the 3-D BSEM response of the model. We compare the result of anomalous responses of three-layered medium derived by 3-D forward modeling with the result of semi-analytical solution. We also perform physical simulation and 3-D numerical forward modeling based on similarity criterion. As a result, the shape of both anomalous response curves are the same, which validates the 3-D numerical simulation method. The anomalous fields of 3-D numerical forward and physical simulation share similar anomalous feature of symmetric bimodal structure that is consistent with its harmonic response curve. It demonstrates that borehole-surface electromagnetic method can be used not only for prospection by employing a multi-frequency response, but also can provide multi-angle information about subsurface anomaly by varying relative depths of vertical finite line source in the borehole.

Efficient multithreaded untransposed, transposed or symmetric sparse matrix–vector multiplication with the Recursive Sparse Blocks format

In earlier work we have introduced the ''Recursive Sparse Blocks'' (RSB) sparse matrix storage scheme oriented towards cache efficient matrix-vector multiplication (SpMV) and triangular solution (SpSV) on cache based shared memory parallel computers. Both the transposed (SpMV_T) and symmetric (SymSpMV) matrix-vector multiply variants are supported. RSB stands for a meta-format: it recursively partitions a rectangular sparse matrix in quadrants; leaf submatrices are stored in an appropriate traditional format - either Compressed Sparse Rows (CSR) or Coordinate (COO). In this work, we compare the performance of our RSB implementation of SpMV, SpMV_T, SymSpMV to that of the state-of-the-art Intel Math Kernel Library (MKL) CSR implementation on the recent Intel's Sandy Bridge processor. Our results with a few dozens of real world large matrices suggest the efficiency of the approach: in all of the cases, RSB's SymSpMV (and in most cases, SpMV_T as well) took less than half of MKL CSR's time; SpMV's advantage was smaller. Furthermore, RSB's SpMV_T is more scalable than MKL's CSR, in that it performs almost as well as SpMV. Additionally, we include comparisons to the state-of-the art format Compressed Sparse Blocks (CSB) implementation. We observed RSB to be slightly superior to CSB in SpMV_T, slightly inferior in SpMV, and better (in most cases by a factor of two or more) in SymSpMV. Although RSB is a non-traditional storage format and thus needs a special constructor, it can be assembled from CSR or any other similar row-ordered representation arrays in the time of a few dozens of matrix-vector multiply executions. Thanks to its significant advantage over MKL's CSR routines for symmetric or transposed matrix-vector multiplication, in most of the observed cases the assembly cost has been observed to amortize with fewer than fifty iterations.

The Study of Thermal-solid Coupling Problems Using Open Source CAE Software

We have been developing the ADVENTURE system, which is an open source CAE software based on parallel finite element method. Two of it‘s main modules are ADVENTURE_Thermal and ADVENTURE_Solid. A thermal-solid coupling system has been developed using these two modules. Three dimensional heat transfer and thermal stresses of a cooling stave in a blast furnace is modeled and analyzed using the developed system. The steady state heat conduction through the model is analysed using the ADVENTURE_Thermal module. Then the temperature distribution is considered as thermal load in the structural analysis. The ADVENTURE_Solid is used for this purpose. Both modules are featured with the parallel finite element method to solve large scale problems. Quite commonly, the matrix that participates in the finite element computation is sparse. Sparse matrix by definition, are populated primarily with many zeroes and thus special storage schemes need to be used to enable efficient storage and computational operations. Skyline storage scheme was used in the modules of old version. In this paper, several sparse matrix storage schemes that show the better performance compare to the Skyline storage are implemented in both modules. Impacts of sparse matrix storage schemes on the total execution time and required memory areevaluated in both thermal and solid analysis of a large scale model.

Advantages of using adaptive remeshing and parallel processing for modelling biodegradation in groundwater

Biodegrading plumes in groundwater are often typified by relatively reactive zones around the fringes and less reactive zones in the core. A high degree of refinement is required at the fringes if a model is to be of use in improving the conceptual understanding of plumes. Two strategies for dealing with the potentially high computational demands are (i) parallel processing, where the workload is shared between multiple processors, and (ii) locally adaptive remeshing, where a refined area of the grid tracks the moving plume fringes through the domain. The partial differential equation toolbox, UG (Unstructured Grids) offers advanced numerical tools including adaptive remeshing, sparse matrix storage schemes, and multigrid solvers. It embraces many of these features within a parallel processing environment. This paper reports on a recent development of UG to simulate field scale reactive biogeochemistry including Monod kinetics, NAPL dissolution, mineral precipitation/dissolution and ion exchange. The non-linear multicomponent reactive transport system is solved with the fully coupled method. Test cases have been used for verification of the new capability. The paper illustrates an application to a 3D field site. It is demonstrated that both adaptive remeshing and parallel processing can improve efficiency and in turn facilitate the incorporation of a more complex set of species and reactions such that understanding of plume processes is enhanced.

Determination of Temperatures and Heat Fluxes on Surfaces and Interfaces of Multidomain Three-Dimensional Electronic Components

A finite element method (FEM) formulation for the prediction of unknown steady boundary conditions in heat conduction for multidomain three-dimensional (3D) solid objects is presented. The FEM formulation is capable of determining temperatures and heat fluxes on the boundaries where such quantities are unknown, provided such quantities are sufficiently overspecified on other boundaries. An inverse finite element program has been previously developed and successfully tested on 3D simple geometries. The finite element code uses an efficient sparse matrix storage scheme that allows treatment of realistic 3D problems on personal computer. The finite element formulation also allows for very straightforward treatment of geometries composed of many different materials. The inverse FEM formulation was applied to the prediction of die-junction temperature distribution in a simple ball grid array electronic package. Examples are presented with simulated measurements, which include random measurement errors. Regularization was applied to control numerical error when large measurement errors were added to the overspecified boundary conditions.

The inverse iteration with the WZ factorization used to the Markovian models

In the article a new approach for solving Markovian linear systems is presented. It is the inverse iteration with the use of the WZ factorization. The paper contains the theoretical base of the method and numerical experiment results where the accuracy and the performance time were investigated. A modified (with the use of the SBLAS library and the Harwell-Boeing sparse matrix storage scheme) method was also presented. The experiments were carried out with the matrices generated for a Markovian model of a leaky bucket with tokens.

An Efficient Finite Element-Based Air Bearing Simulator for Pivoted Slider Bearings using Bi-Conjugate Gradient Algorithms

A computationally efficient air bearing simulator-for pivoted slider bearings has been developed that is based on finite element theory and uses bi-conjugate gradient algorithms in conjunction with a sparse matrix storage scheme. The air bearing simulator involves the simultaneous solution of the Reynolds equation and the slider equilibrium equations. The highly nonlinear problem requires the repeated solution of a nonsymmetric system of equations, typically with a large number of unknowns depending on the complexity of the slider geometry. Iterative solvers, such, as the bi-conjugate gradient algorithms used for this study, require significantly less core memory as compared to direct solvers and reduce the solution time if combined with a suitable preconditions. Of the bi-conjugate gradient algorithm/preconditioner combinations implemented, the Bi-CGSTAB algorithm combined with an ILU preconditioner provided the best performance in terms of smooth convergence and computational efficiency. Presented as a ...

Parallel computation of 3‐D electromagnetic scattering using finite elements

AbstractThe finite element method (FEM) with local absorbing boundary conditions has been recently applied to compute electromagnetic scattering from large 3‐D geometries. In this paper, we present details pertaining to code implementation and optimization. Various types of sparse matrix storage schemes are discussed and their performance is examined in terms of vectorization and net storage requirements. The system of linear equations is solved using a preconditioned biconjugate gradient (BCG) algorithm and a fairly detailed study of existing point and block preconditioners (diagonal and incomplete LU) is carried out. A modified ILU preconditioning scheme is also introducted which works better than the traditional version for our matrix systems. The parallelization of the iterative sparse solver and the matrix generation/assembly as implemented on the KSR1 multiprocessor is described and the interprocessor communication patterns are analysed in detail. Near‐linear speed‐up is obtained for both the iterative solver and the matrix generation/assembly phases. Results are presented for a problem having 224,476 unknowns and validated by comparison with measured data.

Power system modelling and sparse matrix operations using object-oriented programming

This paper reports on power system modelling and sparse matrix operations using object-oriented programming (OOP). It has been claimed that OOP leads to more flexible, modular and reusable code, and that programs can be written more generally. The main focus of the paper is OOP design principles and practical implementations for power systems. Specific examples included demonstration of a power system model design, particularly focusing on OOP mechanisms and object-oriented style of programming. Also demonstrated, using OOP, are tailor-made sparse matrix storage schemes and operations. Numerical tests indicate that the proposed design is efficient compared to standard numerical library routines, and that the particular OOP features that are used to obtain flexibility etc. do not significantly increase computation time. >

Amélioration des performances d'une implantation parallélovectorielle du gradient conjugué par extension du schéma de stockage matriciel

Electromagnetic field computation with the Finite Element (FE) method implies solving of large linear systems of equations. Performances and memory capacities of today computers allow to achieve three-dimensional FE discretizations of electromagnetic problems, but the number of unknowns grows high. So, to improve time to the numerical solution of the linear system(s) thus arising, the use of parallel and/or vector computers has to be envisaged. In this paper, the main constitutive steps of the Pre-conditioned Conjugate Gradient algorithm (PCG) are analysed. After a short recall of our previous work concerning their improvement by use of vector and parallel computations, we show some speedup limitations due to the sparse row-wise matrix storage scheme employed. Then, an extension of this matrix representation is proposed, leading to introduce redundant storage of non-zero coefficients. In spite of the «memory waste» thus implied, it is shown how this extension can be successfully employed to icnrease the speedup due to parallelism and vectorization on the whole algorithm, and in particular to derive a parallel preconditioner

Decomposition of sparse matrices on parallel computers

Sparse matrices are an integral part of finite element analyses and an efficient parallel matrix decomposition strategy is essential to high performance finite element systems. To support this need a study is carried out on the performance of parallel algorithms for decomposing sparse matrices for substructure and finite element problems. Results show that the sparse matrix storage schemes studied have a high degree of parallelism in the process of the sparse matrix decomposition but require careful attention to the management of cache memories for a massively parallel computer in which cache processing yields performance gain.

An evaluation of speedup in conjugate gradient routines with a mathematical vector library

The experience in vectorizing finite element and finite difference electromagnetic analysis codes is reported. Techniques to efficiently use vector machines in problems involving sparse, unstructured matrices are discussed. In this context a sparse matrix storage scheme particularly suited to vector architectures is introduced. The results of experience in the use of a library of conjugate-gradient-type solvers suited to vector architectures implemented according to this matrix structure are reported and discussed. >

Comparing direct and iterative equation solvers in a large structural analysis software system

This paper describes two vectorized implementations of preconditioned conjugate gradient (PCG) solvers. Sparse and diagonal matrix storage schemes are described and compared. A vectorized incomplete Choleski preconditioning is described and compared with Jacobi preconditioning. A modification to the basic no-fill incomplete Choleski method to improve performance and robustness is described. The two PCG solvers are compared with direct Choleski methods using a sparse Choleski solver from SPARSPAK and a vectorized variable-band Choleski solver developed at NASA Langley Research Center. All of the linear equation solvers are implemented in a large structural analysis finite element software system called the Computational Structural Mechanics (CSM) Testbed. The CSM Testbed is used to provide a common software system in which new methods are developed and tested. Several representative two- and three-dimensional structural analysis problems are solved using the various equation solvers. Results are given from runs made on the CONVEX C220 and CRAY 2 computer systems. Comparisons of the convergence rates for the iterative solvers as well as the computation rates, number of operations, and overall CPU time required by all of the equation solvers are given.

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.

Approximation of the Consecutive Ones Matrix Augmentation Problem

In this publication we will prove a number of negative results concerning the approximation of the NP-complete CONSECUTIVE ONES MATRIX AUGMENTATION problem. We will characterize a large class of simple algorithms that do not find a near optimum augmentation of their input matrices. We will show that there are matrices for which these algorithms find augmentations that are even far from optimal. These results are important for the analysis of a sparse matrix storage scheme.

Application of sparse matrix techniques to inter-regional input-output analysis

Combining the concepts described in this paper makes it possible to solve inter-regional input-output systems (and other types of large, sparse, linear systems) with considerable efficiency in storage and computation. The exact number of operations and corresponding savings in computational time and storage depend on the particular zero-non zero structure of each matrix in the system, but in any case the savings can be enormous. A recommended procedure is summarized below. 1. Take the entire inter-regional IO system which includes the representation of each of the regions and of the links among them and express it in the formMx=Bz or, more simply,Mx=b (since bothB andz are known). It may be helpful for this purpose to draw a picture of the large matrixM. 2. PartitionM into blocks, exploiting the structure of the particular system. For example, at the very least, the matrix or matrices corresponding to each region will probably be separate blocks. The analysis required for this step may lead to reformulating the matrix representation of the given economic system by, for example, replacing a set of equations with linear combinations of these same equations, particularly for the equations representing the links among regions. 3. Identify an appropriate partition of the blocks ofM, and a corresponding partition of the vectorsx andb, for performing a block factorization. This solution algorithm, which solves forx inMx=b (whereM now represents the entire system) is not available as a packaged program and so for the foreseeable future must be written in ‘home-made’ computer code that assumes a sparse matrix storage scheme compatible with the package to be used (as in 4 below). The algorithm for block factorization in the 2×2 case relevant to our particular inter-regional IO system is given in Section 8 of this paper. 4. Some steps of the algorithm require solving smaller systems of the formHu=v foru, given some matrixH and some vectorv. Efficient packaged subroutines can be obtained at low cost in order to solve these systems by finding the block triangular form corresponding to a particularH, then performing theLU factorization of the diagonal blocks only. The home-made code will interact with these subroutines.