Hybrid GPU Matrix Formats

• Hybrid matrix formats for GPUs are advanced storage strategies that adaptively combine dense, sparse, and low-rank representations to overcome the limitations of single-format approaches.
• They leverage structured decomposition and dynamic runtime selection—using methods like ELLPACK+COO, H-matrices, and AR-CSR—to optimize arithmetic intensity and memory layout.
• Empirical results demonstrate that multi-format adaptivity can boost throughput significantly (up to 10× speedup), benefiting scientific and engineering computations on GPUs.

Hybrid matrix formats for GPUs are storage and computation strategies that combine multiple data representations—typically mixing dense, sparse, low-rank, or block formats—to maximize arithmetic intensity, memory coalescing, load balance, and throughput in key linear algebra operations. These approaches are central in high-performance GPU computing for scientific and engineering applications, enabling efficient solution of problems ranging from boundary integral equations to electronic structure simulations and large-scale iterative solvers. Hybrid formats mitigate the limitations of single-format approaches by adaptively selecting or combining formats based on sparsity, structure, or hardware characteristics.

1. Fundamental Principles of Hybrid Matrix Formats

Hybrid matrix formats on GPUs leverage the strengths of multiple storage representations to overcome the inefficiencies that arise when sparse or structured matrices are mapped to GPUs using naive schemes. The essential principle is decomposition of the matrix into segments (by block, row, or cluster) that are each mapped to the representation most conducive to GPU arithmetic or memory hierarchy (Oberhuber et al., 2010, Heller et al., 2012, Harbrecht et al., 2018, Stylianou et al., 2022, Fattebert et al., 24 Jan 2024). For example:

• "Hybrid" ELLPACK+COO stores the regular part of a sparse matrix in a column-major ELLPACK block for optimal coalescing, and irregularly sized row tails in COO to avoid excessive padding.
• $H$-matrices decompose dense operator matrices into a mix of small dense near-field blocks and far-field blocks that are recompressed as low-rank products.
• Adaptive Row-grouped CSR (AR-CSR) partitions rows into groups, minimizing padding and maximizing coalescing by redistributing threads according to row workload.

By exploiting such structural heterogeneity, hybrid formats reduce global memory traffic, minimize padding, and balance the arithmetic workload across GPU thread groups.

2. Core Data Structures and GPU Implementations

Hybrid matrix formats are instantiated via flexible data structures that support efficient memory layout and kernel dispatch. Key design patterns include:

• Interleaved arrays for regular blocks and separate scatter-gather arrays for irregular regions (as in ELLPACK+COO hybrid) (Oberhuber et al., 2010).
• Batched storage of low-rank factors ($U, V$) and dense blocks for $H$-matrices on GPUs, using flat arrays and pointer tables to organize contiguous batched operations (Harbrecht et al., 2018).
• Group info tables and adaptive per-group chunk sizing in AR-CSR, supporting thread-to-row assignments that ensure coalescing and balanced workload among threads (Heller et al., 2012).
• Opaque format-agnostic matrix handles (as in BML) or dynamic containers (as in Morpheus) that trigger fast on-device format conversions and runtime selection (Stylianou et al., 2022, Fattebert et al., 24 Jan 2024).

The implementation relies on fine-tuned CUDA kernels or OpenMP offloading targeting regular regions for vector-friendly BLAS routines and irregular/small regions for custom scatter or atomic-update kernels.

3. Representative Hybrid Matrix Formats

The following table summarizes several prominent hybrid matrix formats used in GPU linear algebra, with their key characteristics and effective regimes:

Format/Library	Composition	Best suited for
ELLPACK + COO (Hybrid)	ELL for first K entries per row, COO for tail	Low-variance row lengths (Oberhuber et al., 2010)
$H$-matrix	Dense near, low-rank far (ACA)	Boundary integral, kernel matrices (Harbrecht et al., 2018)
Adaptive Row-Grouped CSR	CSR partitioned into load-balanced groups	High row-nnz variability (Heller et al., 2012)
BML/PROGRESS	Dense, ELLPACK-R, CSR, ELLBLOCK	DFT, Tight-Binding, O(N) solvers (Fattebert et al., 24 Jan 2024)
Morpheus DynamicMatrix	CSR, COO, DIA, hybrid split	Runtime-adaptive SpMV (Stylianou et al., 2022)

In the ELLPACK+COO format, a cutoff $K_{1}$ is selected to trade off between ELLPACK padding and COO atomic overhead. For $H$-matrices, the block-cluster tree structure encodes both the dense and low-rank segments, where admissibility is determined by geometric considerations and ACA constructs the low-rank decompositions directly on the GPU.

4. Algorithms and Runtime Format Selection

Hybrid formats often make dynamic decisions at setup or runtime regarding which format to use for each matrix segment:

• AR-CSR partitions are formed such that each group has similar total nnz; within a group, thread allocation is adapted to per-row workload, and per-thread chunk sizes are minimized to maximize memory utilization and minimize padding (Heller et al., 2012).
• Morpheus employs a profiling-based auto-tuner: it times small numbers of SpMV trials per candidate format (CSR, COO, DIA), then selects the format with minimal runtime for each segment (Stylianou et al., 2022).
• BML+PROGRESS maintains the original format in memory but converts as needed before dispatching compute-intensive kernels to vendor libraries (e.g., ELLPACK→CSR for sparse-matrix multiply with ROCSPARSE, or to Dense for MAGMA/CUBLAS) (Fattebert et al., 24 Jan 2024).

A plausible implication is that as format conversion overheads diminish (e.g., sub-millisecond on modern GPUs), runtime adaptivity becomes increasingly attractive, especially for iterative solvers or applications with mixed regularity.

5. Performance Characteristics and Scalability

Empirical results demonstrate substantial performance advantages for hybrid formats in GPU SpMV and other linear algebra tasks:

• $H$-matrix GPU solvers achieve strong scaling efficiency of ≈67% when increasing from 128 to 1024 GPUs, solving 1.5M-unknown BEM systems in <6 minutes; arithmetic throughput is dictated by block size and batching efficiency, with bottlenecks arising when blocks are too small for high occupancy (Harbrecht et al., 2018).
• AR-CSR outperforms cuSPARSE-CSR, Hybrid, and Row-grouped CSR on 1168/1600 UF matrices, with up to 10× speedup for matrices with high std/mean row-nnz (>1.5), and low-overhead for regular matrices (≤1.2× nnz storage overhead) (Heller et al., 2012).
• In Morpheus, DIA format yields up to 4.5× speedup over CSR for large 3D-Poisson-like matrices in SpMV, and the hybrid local/remote split maintains optimal SpMV throughput even as per-GPU problem size shrinks with strong scaling (Stylianou et al., 2022).
• BML+PROGRESS achieve 5–10× speedups for dense SP2 mat-muls vs. traditional diagonalization, with sparse-matrix mat-mul switching to ROCSPARSE CSR or hypre as appropriate, depending on sparsity and hardware (Fattebert et al., 24 Jan 2024).

A plausible implication is that multi-format adaptivity, especially when coupled with batching and arithmetic intensity optimization, is the predominant path to practical peak performance on modern GPU clusters.

6. Integration, API Strategies, and Lessons Learned

Modern GPU hybrid matrix libraries emphasize unified abstractions and minimal user intervention:

• OpenMP 4.5+ (as in BML) enables data migration and targeting across CPU, GPUs (Nvidia/AMD/Intel), and supports format-agnostic APIs with transparent invocation of vendor libraries for high-performance regions (Fattebert et al., 24 Jan 2024).
• Morpheus exposes a DynamicMatrix container and automatic or user-activated format switching, hiding all format selection and conversion details behind a small number of API calls; conversions proceed entirely on-device via COO intermediates (Stylianou et al., 2022).
• For distributed-memory environments, $H$-matrix codebases replicate critical data (trees, input vectors) across GPUs and use CUDA-aware MPI reductions for collectives, avoiding complex communication patterns (Harbrecht et al., 2018).

Observed best practices include:

• Selecting block or leaf sizes to maximize arithmetic intensity and occupancy.
• Grouping like-sized subblocks for batched DGEMV/GEMM calls.
• Exploiting space-filling curves or Morton ordering to simplify cluster construction for recursive formats.
• Providing minimal callback-based interfaces for application integration (e.g., for customized quadrature evaluation in BEM solvers).

The ability to seamlessly transition between multiple hybrid storage backends, as exemplified by BML, Morpheus, and recent distributed $H$-matrix solvers, is now fundamental for architecturally portable, high-performance GPU-based linear algebra.

7. Comparative Assessment and Format Selection Guidelines

Correct hybrid format selection is dependent on matrix structure and computational regime. Key selection guidelines, as reported:

• ELLPACK+COO hybrids are favored for matrices with low coefficient of variation (std/mean <0.3) in per-row nnz; choose K₁ to contain ≳90% of the nonzeros in ELLPACK for maximal coalescing (Oberhuber et al., 2010).
• For irregular/heterogeneous matrices (row-nnz Gini >0.3 or std/mean >0.5), AR-CSR or similar group-adaptive formats provide highest throughput by load-balancing and minimizing padding (Heller et al., 2012).
• $H$-matrix and low-rank hybrid structures dominate in boundary integral and kernel-based problems where asymptotically smooth blocks permit rank compression and dense batched operations in the near field (Harbrecht et al., 2018).
• Runtime/adaptive approaches, leveraging on-device profiling or cost models, increasingly supersede static format mapping as conversion overheads decline and solver iterations rise (Stylianou et al., 2022, Fattebert et al., 24 Jan 2024).

The following table presents format adequacy as a function of structural characteristics:

Matrix property	Preferred format(s)
Uniform row-nnz	ELLPACK, Hybrid (ELL+COO)
High row-nnz variance	AR-CSR, segment hybrid
Block structure	Block-ELL, ELLBLOCK
Asymptotically smooth	$H$-matrix (low-rank, dense)
Unknown, mixed, runtime	Dynamic/auto-tuned (Morpheus)

The adaptability and composability of hybrid matrix formats underpin their preeminence in scientific GPU computing, and they offer a path to scalable, performant solutions across a diversity of problem domains and hardware platforms (Harbrecht et al., 2018, Heller et al., 2012, Stylianou et al., 2022, Fattebert et al., 24 Jan 2024, Oberhuber et al., 2010).