PaperlayerView exact revision on arXiv

Source-linked AI summary

A unified sparse matrix data format for efficient general sparse matrix-vector multiply on modern processors with wide SIMD units

Moritz Kreutzer, Georg Hager, Gerhard Wellein, Holger Fehske, Alan R. Bishop

arXiv:1307.6209v2cs.MScs.DC

TL;DR

Sparse matrix-vector multiplication is difficult to optimize across heterogeneous hardware because storage formats and SIMD usage are architecture-sensitive. The paper proposes SELL-C-σ, analyzes its parameters and transfer behavior, and evaluates it across CPUs, Xeon Phi, and GPGPU systems. Across most matrices, it provides best or highly competitive performance with one format, within a single-chip scope and stated modeling assumptions.

  • Problem

    spMVM is often the most time-consuming kernel, yet performance depends on irregular sparsity, indirect vector accesses, and hardware-sensitive storage formats.

  • Method

    SELL-C-σ combines Sliced ELLPACK with chunk-size and row-sorting parameters to support SIMD/SIMT execution across heterogeneous processors.

  • Results

    For most matrices, one SELL-C-σ format incurs no significant performance loss versus hardware-specific formats and significantly outperforms CRS on most Intel Xeon Phi matrices.

  • Takeaways & Limitations

    A unified format simplifies matrix storage and supports straightforward use of hybrid programming models across heterogeneous compute systems.

  • Takeaways & Limitations

    The study is limited to single-chip systems and OpenMP threading, while its performance upper bound assumes memory bandwidth is the only limiting factor and ignores latency effects.

Abstract

from arXiv · show

Sparse matrix-vector multiplication (spMVM) is the most time-consuming kernel in many numerical algorithms and has been studied extensively on all modern processor and accelerator architectures. However, the optimal sparse matrix data storage format is highly hardware-specific, which could become an obstacle when using heterogeneous systems. Also, it is as yet unclear how the wide single instruction multiple data (SIMD) units in current multi- and many-core processors should be used most efficiently if there is no structure in the sparsity pattern of the matrix. We suggest SELL-C-sigma, a variant of Sliced ELLPACK, as a SIMD-friendly data format which combines long-standing ideas from General Purpose Graphics Processing Units (GPGPUs) and vector computer programming. We discuss the advantages of SELL-C-sigma compared to established formats like Compressed Row Storage (CRS) and ELLPACK and show its suitability on a variety of hardware platforms (Intel Sandy Bridge, Intel Xeon Phi and Nvidia Tesla K20) for a wide range of test matrices from different application areas. Using appropriate performance models we develop deep insight into the data transfer properties of the SELL-C-sigma spMVM kernel. SELL-C-sigma comes with two tuning parameters whose performance impact across the range of test matrices is studied and for which reasonable choices are proposed. This leads to a hardware-independent ("catch-all") sparse matrix format, which achieves very high efficiency for all test matrices across all hardware platforms.

1. Introduction and Related Work.

Sparse matrix-vector multiplication is a central, often time-consuming kernel whose performance depends strongly on sparsity patterns, hardware, and storage format. The paper proposes SELL-C-σ as a unified SIMD-oriented format for heterogeneous systems.

  • Motivation: spMVM frequently dominates numerical algorithms built on sparse linear algebra, including eigenvalue and sparse-system solvers.Its efficient implementation is therefore highly important across science and engineering applications.
  • Motivation: Irregular nonzero patterns complicate performance prediction through indirect RHS-vector accesses, while storage formats remain hardware-sensitive.Different formats optimize CPUs, vector computers, and GPGPUs, making heterogeneous execution more difficult.
  • Contribution: SELL-C-σ builds on Sliced ELLPACK to provide one storage format for diverse processor designs, extending its prior use beyond GPGPUs.The format targets SIMD vectorization and heterogeneous compute clusters.
  • Contribution: Choosing chunk size C to match SIMD width and sorting rows within scope σ reduces overhead and improves performance when σ is not too large.On GPGPUs, C corresponds to the number of threads per warp.
  • Contribution: Performance modeling explains how C and σ interact with chunk occupancy and nonzeros per row, while experiments assess portability across Sandy Bridge, Xeon Phi, and Nvidia K20.Fixed C and σ deliver best or highly competitive performance for most tested matrix types.
  • Scope: The study is limited to single-chip systems and OpenMP threading, and assumes general matrices without exploiting special substructures.MPI, hybrid MPI+X parallelization, blocking, and unrolling remain outside the paper’s scope.

2. Hardware and test matrices.

The evaluation spans CPU, many-core, and GPGPU architectures with distinct SIMD or SIMT capabilities and uses matrix collections chosen to expose layout-sensitive behavior. Hardware bandwidth and cache characteristics provide baselines for interpreting spMVM performance.

  • Hardware: The test systems represent cache-based x86 multicore processors, massively parallel Intel Phi accelerators, and Nvidia GPGPU architectures.These platforms cover contrasting core counts, clock speeds, SIMD widths, and threading models.
  • Performance baselines: Bandwidth baselines use array-copy and read-accumulate benchmarks to represent unfavorable and upper-limit memory-transfer cases.These baselines are intended to provide sensible references for different matrix types.
  • Hardware: Intel Sandy Bridge uses 256-bit AVX registers, whereas Intel Xeon Phi uses 512-bit SIMD units and requires vectorization for reasonable performance.The Phi can process eight double-precision or 16 single-precision fused multiply-add operations per instruction.
  • Hardware: Nvidia K20 executes most hardware operations in SIMT groups of 32 threads called warps, while reductions can constrain read-only performance.The study models a warp as a SIMD execution unit with width 2048 bits for double precision.
  • Test matrices: Detailed analysis uses four matrices selected as corner cases, then validates broader applicability with 14 matrices from the Williams group collection.The selected matrices represent characteristics that strongly influence data-layout efficiency.
  • Test matrices: The benchmark descriptions include dimensions, nonzero counts, average nonzeros per row, density, and chunk occupancy with and without sorting.Chunk occupancy is reported for C = 16 and σ values of 1 and 256.

3. Matrix formats and spMVM kernels.

This section compares sparse matrix storage formats and their spMVM kernels, emphasizing SIMD/SIMT suitability, padding, row sorting, and indirect vector access. SELL-C-σ combines chunking and limited row sorting to reduce storage overhead while preserving locality across diverse architectures.

  • 3. Matrix formats and spMVM kernels.: CRS stores nonzeros row by row, but SIMD vectorization incurs horizontal-reduction and remainder-loop overhead that grows with SIMD width.The compiler can vectorize bulk iterations, although alignment peeling and scalar cleanup reduce efficiency, especially for short rows.
  • 3.2. Analysis of the CRS format.: On wide-SIMD processors, average nonzeros per row must substantially exceed SIMD width because masked remainder handling and alignment constraints add overhead.Intel Phi processes 16 single-precision or eight double-precision values per 512-bit instruction.
  • 3.2. Analysis of the CRS format.: CRS is unsuitable for GPGPUs: inner-loop parallelization incurs warp reductions, while outer-loop parallelization destroys coalesced memory access.The passage concludes that CRS is a poor choice on GPGPUs under either execution strategy.
  • 3.3. Sliced ELLPACK and SELL-C-σ.: ELLPACK enables coalesced row-wise thread access through column-major storage and equal row lengths, while Sliced ELLPACK reduces its padding overhead using row chunks.Within each chunk, rows are padded only to the longest row and stored consecutively in column-major order.
  • 3.3. Sliced ELLPACK and SELL-C-σ.: SELL-C-σ aligns chunk size C with SIMD-register or warp width and sorts rows locally within scope σ to increase β without globally disrupting RHS-vector locality.Global sorting can approach β ≈ 1 but may destroy physical-problem locality and increase code balance.
  • 3.4. Analysis of the SELL-C-σ format.: Row reordering generally requires corresponding column-index permutation because iterative solvers may operate in permuted index space and reordering can escalate matrix bandwidth.These changes can affect RHS-vector access patterns and locality.
  • 3.4. Analysis of the SELL-C-σ format.: On GPGPUs, one thread per row can stop at the actual row length, but low-β chunks still occupy warp resources until their longest-running thread finishes.This optimization makes SELL-C-σ equivalent to Sliced ELLR-T with T = 1.
  • 3.5. General performance issues of spMVM.: Irregular RHS-vector access can sharply reduce performance by defeating cache locality or load coalescing, while x86 lacks gather instructions and MIC gather benefits depend on row locality.Bandwidth-reduction transformations such as RCM are outside the work’s scope.

4. Performance models.

The performance models express spMVM cost through memory traffic, accounting for matrix-format padding, RHS locality, and LHS updates. Under stated assumptions, they estimate an upper bound on attainable performance and relate SELL-C-σ to CRS.

  • Code balance: The code balance counts bytes transferred over the memory interface per floating-point operation in spMVM.For square matrices, the model derives this quantity from the kernel's data traffic.
  • Traffic components: Matrix traffic includes values and column indices, while LHS traffic contributes 16 bytes/Nnzr per non-zero row element.The model assumes double-precision matrix and vector data with four-byte integer indices.
  • SELL-C-σ model: SELL-C-σ padding increases average matrix traffic by the reciprocal chunk occupancy β, while padded entries do not add RHS traffic.The corresponding CRS model is recovered by setting β = 1.
  • Performance bound: The roofline model predicts maximum spMVM performance from achievable memory bandwidth and code balance.For sufficiently large Nnzr, the model yields the best attainable performance level considered in this work.
  • Model assumptions: The resulting estimates are upper bounds because they assume bandwidth-only limitation, infinitely fast cache access, optimal replacement, and no latency effects.These assumptions may not hold in reality, although the bound applies when matrix data comes from main memory.

5. Performance results and analysis.

Across Intel Sandy Bridge, Intel Xeon Phi, and Nvidia K20, SELL-C-σ generally delivers the strongest spMVM performance when chunk size and sorting scope are tuned, while irregular matrices expose limits from load imbalance, short rows, and cache behavior.

  • Unified data layout performance: SELL-C-σ with optimal sorting attains best performance on all architectures, with the largest sorting impact for matrices having poor chunk occupancy.RM07R and kkt power are cited as examples.
  • Unified data layout performance: 1.5× to 4×: vectorized SELL-C-σ boosts Intel Phi performance over vectorized CRS for matrices that fit into its LLC.SELL-C-σ also substantially outperforms CRS for cache-resident matrices on Intel platforms.
  • Performance limits: 7.2 GF/s, 27.5 GF/s, and 25.2 GF/s: modeled memory-bound limits for Intel SNB, Intel Phi, and Nvidia K20, respectively, when β = 1.For matrices with Nnzr > 50, Intel Phi and Nvidia K20 reach around 80% of their limits, while Intel SNB exceeds 90%.
  • Performance limits: Low Nnzr reduces performance by roughly twofold or more across architectures, as observed for Hamrle3 and kkt power in memory-bound cases.The performance model attributes this decline to short rows; cache-bound matrices fall outside the model’s assumptions.
  • Pathological and irregular matrices: Accelerator performance declines for rail4284, webbase-1M, dense2, and scircuit because limited parallelism, load imbalance, irregular access, or low chunk occupancy constrain execution.Multiple threads per row can improve parallelism in problematic accelerator cases, while HYB dominates SELL-C-σ for scircuit on Nvidia K20.
  • Overall conclusion: SELL-C-σ achieves high performance across the tested architectures and matrix types, supporting its role as a unified sparse-matrix storage format.The format combines Sliced ELLPACK with SIMD vectorization and is reported as suitable for Intel SNB, Intel Xeon Phi, and Nvidia Kepler.

6. Conclusions and outlook.

SELL-C-σ is presented as a single sparse-matrix format for heterogeneous systems, maintaining near hardware-specific performance across most investigated matrices. The authors identify future work on accelerator support, MPI integration, and automatic tuning.

  • SELL-C-σ enables heterogeneous algorithms to store matrices in a single format without significant performance loss for most investigated matrices.
  • On Intel Xeon Phi, SELL-C-σ significantly outperforms CRS on most matrices.
  • Future work targets small-row-count or high-row-density matrices through multiple threads per row, plus MPI support and automatic tuning.

Appendix A. Description of the corner case benchmark matrices.

The appendix describes a detailed storage-format performance analysis based on four matrices from the University of Florida Sparse Matrix Collection.

  • The detailed performance analysis uses four matrices from the University of Florida Sparse Matrix Collection.
  • The appendix introduces matrix descriptions drawn from the same University of Florida collection.
  • These matrices form the benchmark basis for analyzing various sparse storage formats.

(a) RM07R

RM07R comes from a computational-fluid-dynamics finite-volume discretization of a three-dimensional viscous case with frozen turbulence.

  • RM07R arises from a CFD finite-volume discretization.
  • The matrix represents a three-dimensional viscous case.
  • The modeled turbulence is frozen.

(b) kkt power

The kkt power matrix arises from nonlinear optimization for finding an optimal power flow.

  • kkt power comes from a nonlinear optimization problem.
  • The optimization is associated with Karush-Kuhn-Tucker conditions.
  • The application goal is finding the optimal power flow.

(c) Hamrle3

Hamrle3 is a matrix from a very large electrical network simulation.

  • Hamrle3 originates from a very large electrical network simulation.
  • The matrix represents an electrical network simulation case.
  • Its stated provenance is a large-scale simulation application.

(d) ML Geer

ML Geer is a matrix obtained to determine the deformed configuration of an axial-symmetric porous medium under pore-pressure drawdown.

  • ML Geer was obtained to find the deformed configuration of an axial-symmetric porous medium.
  • The medium is subject to a pore-pressure drawdown.
  • The matrix was generated through a meshless Petrov-Galerkin discretization.

Appendix B. Code listings.

The appendix lists SELL-C-σ kernels implemented with Intel MIC and AVX intrinsics, using chunked loops, SIMD loads, gathered right-hand-side values, accumulation, and stores.

  • Intel MIC kernel: The Intel MIC implementation uses SELL-C-σ with C = 16 for 64-bit values and 32-bit indices.
  • Intel MIC kernel: The MIC kernel loads matrix values and indices, gathers right-hand-side values, then multiplies and accumulates them.
  • Intel MIC kernel: The MIC kernel stores two groups of 8 left-hand-side values after accumulation.
  • Intel AVX kernel: The AVX implementation uses SELL-C-σ with C = 4 for 64-bit values and 32-bit indices.
  • Intel AVX kernel: The AVX kernel loops over chunks, loads 4 left-hand-side values and 4 matrix values, assembles right-hand-side values, and accumulates.
Loading 1307.6209v2…