PACK SELL: A S PARSE M ATRIX F ORMAT FOR P RECISION -AGNOSTIC H IGH -P ERFORMANCE S P MV
arXiv:2604.13433v1 [cs.DC] 15 Apr 2026
A P REPRINT Kengo Suzuki Academic Center for Computing and Media Studies, Kyoto University, Kyoto, Japan [email protected]
Takeshi Iwashita Academic Center for Computing and Media Studies, Kyoto University, Kyoto, Japan [email protected]
April 16, 2026
A BSTRACT We propose a new sparse matrix format, PackSELL, designed to support diverse data representations and enable efficient sparse matrix–vector multiplication (SpMV) on GPUs. Building on sliced ELLPACK (SELL), PackSELL incorporates delta encoding of column indices and a novel packing scheme that stores each index-delta–value pair in a single word, thereby reducing memory footprint and data movement. This design further enables fine-grained control over the bit allocation between deltas and values, allowing flexible data representations, including non-IEEE formats. Experimental results show that, when configured for half precision (FP16), the PackSELL-based SpMV kernel outperforms the cuSPARSE SELL-based kernel by up to 1.63×. Moreover, with configurations using customized formats, PackSELL achieves FP32-level accuracy while exceeding the performance of FP16 cuSPARSE. These benefits extend to sparse linear solvers; for example, a mixed-precision preconditioned conjugate gradient (PCG) solver using PackSELL achieves up to a 2.09× speedup over the standard full-precision PCG. Keywords SpMV · mixed precision · FP16 · non-IEEE formats · Krylov subspace methods · GPUs
1
Introduction
Sparse matrix-vector multiplication (SpMV) is a fundamental kernel in a wide range of scientific applications and has been extensively studied (Gao et al., 2024). In this study, we consider the efficient execution of SpMV on GPUs for real matrices, y = Ax, A = [ai,j ] ∈ Rn×m , x ∈ Rm , y ∈ Rn , (1) for different low-precision data representations. In particular, we develop a new sparse matrix format that efficiently and uniformly supports various data representations, including integer and floating-point formats not defined in IEEE 754. Recent hardware trends have widened the gap between memory bandwidth and computing performance, especially with advancements in low-precision arithmetic units (Lindquist, 2023; Abdelfattah et al., 2021). As a result, reducing data movement has become critical for memory-bound kernels such as SpMV. Consequently, mixed-precision techniques, particularly those using low precision, have attracted attention (Abdelfattah et al., 2021). A prominent example is iterative solvers for sparse linear systems, where the IEEE 754 single-precision (FP32) format has become widely used in addition to the double-precision (FP64) format (Buttari et al., 2008; Ikuno et al., 2012; Nakajima et al., 2021; Lindquist et al., 2022; Zhao et al., 2022, 2023; Amestoy et al., 2024; Spyropoulos and Antonopoulos, 2025). More recently, state-of-the-art methods exploit the half-precision floating-point (FP16) format (Suzuki and Iwashita, 2025). Furthermore, non-IEEE floating-point (Ichimura et al., 2018; Hunhold and Quinlan, 2025) and few-bit integer (or fixed-point) (Jerez et al., 2015; Iwashita et al., 2020; Suzuki et al., 2025) representations have also been investigated. These trends suggest the growing significance of supporting various low-bit formats in memory-bound kernels, including SpMV. However, existing high-performance SpMV kernels are still largely designed for classical higher-precision
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
IEEE formats such as FP64 and FP32, and the investigation of other formats remains limited. This is mainly due to the tight coupling between data representation and the efficiency of memory access: conventional sparse matrix formats implicitly assume data types aligned to hardware-friendly boundaries, typically multiples of bytes. Although several custom formats are designed using memory accessors to better utilize memory bandwidth (Anderson and Gregg, 2016; Mukunoki et al., 2023; Graillat et al., 2024a), they are still constrained by such alignment requirements, and the variability of available data representations is limited. Consequently, many potentially useful non-byte-aligned formats remain difficult to use efficiently. We argue that this limitation stems from treating data representation as a fixed attribute when designing sparse matrix formats. That is, the limitation can be overcome by jointly designing data representation together with other attributes of sparse matrix formats, such as index encoding and memory layout. To this end, we propose a new sparse matrix format based on the Sliced ELLPACK (SELL) format (Monakov et al., 2010), called delta–value packing-based SELL (PackSELL). The proposed format packs each nonzero element together with its column index encoded as the difference (delta) from the preceding element into a single word; these packed data are reconstructed on the fly during operations such as SpMV. By flexibly adjusting the bit allocation between deltas and nonzero values, PackSELL supports a wide range of data representations, including non-byte-aligned formats, while reducing memory footprint through delta encoding and packing. Therefore, PackSELL is expected to enable both flexible data representations and high-performance SpMV. We evaluate a PackSELL-based GPU SpMV kernel across a range of bit allocations between values and deltas. In an FP16-equivalent setting, our method outperforms conventional SpMV kernels, including NVIDIA cuSPARSE and the state-of-the-art DASP (Lu and Liu, 2023). We further examine the trade-off between performance and accuracy by varying the bit allocation, showing that FP32-level accuracy can be achieved with FP16-level performance. Finally, we integrate our kernel into mixed-precision Krylov subspace methods and demonstrate overall speedups and a practical scenario where non-IEEE formats realized by PackSELL provide performance benefits. The contributions of this study are summarized as follows: • A novel sparse matrix format, PackSELL, enabling high-performance SpMV on GPUs across various data formats. • Performance analysis of SpMV kernels and Krylov subspace methods using PackSELL. • Identification of a scenario where non-IEEE representations benefit mixed-precision Krylov subspace methods.
2
Related Work
Since it is practical, SpMV has been researched on various architectures, including GPUs (Filippone et al., 2017; Gao et al., 2024). Previous studies can be broadly classified into two approaches: improving execution efficiency (e.g., load balancing and memory access optimization) and reducing memory footprint via data compression, although some studies address both simultaneously (Aliaga et al., 2022). Both approaches are effective, but they generally target different types of matrices. The former focuses on matrices with irregular sparsity patterns, where workload imbalance is a major bottleneck, especially on GPU SIMT architecture. To address this, formats such as compressed sparse row (CSR) and coordinate (COO), as well as their variants, are widely adopted, and numerous optimization techniques have been proposed (Kreutzer et al., 2014; Anzt et al., 2014; Bell and Garland, 2009; Ashari et al., 2014; Liu and Vinter, 2015; Steinberger et al., 2017; Anzt et al., 2020; Niu et al., 2021). Adaptive techniques that hybridize multiple sparse matrix formats or select multiple sub-kernels based on matrix characteristics have further improved performance (Lu and Liu, 2023). However, these strategies often provide limited benefits for matrices with relatively regular sparsity patterns, as suggested in (Zhang et al., 2025). For such matrices, straightforward implementations using hardware-friendly formats such as SELL can already achieve near-peak performance, as bounded by memory bandwidth. In this case, reducing memory footprint is critical to boost performance. Moreover, since such matrices typically arise in applications such as sparse linear solvers, this direction is also important in practice. One strategy in this direction is to compress indices. Traditional techniques include blocking methods, such as block sparse row (BSR) (Barrett et al., 1994) and block-based SELL variants, which reduce index size by grouping nonzero elements into blocks. While effective for structured matrices, these methods can introduce excessive padding and degrade performance for other matrices. Column-wise blocking has also been proposed for a similar purpose (Nagasaka et al., 2016), and recent work continues to refine such strategies (Cong et al., 2025). Other techniques compress indices, for example, using dictionary compression (Murakami et al., 2026) or encoding differences between neighboring 2
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
Slice Offsets (offset)
0 1 2 3 4 5 6 0 a b 1 c
d
2
e
3
f
4
h
g
Sample Matrix
A P REPRINT
a b
0 1
0 4 8 10
c d
0 4
Values (val)
e -
5 -1
f g
1 3
h
1
a c b d e f - g h Column Indices (col) 0 0 1 4 5 1 -1 3 1 -1
Compressed Representation
SELL Format
Figure 1: Example of the SELL format with slice size C = 2. elements (Maggioni and Berger-Wolf, 2014). Although the latter is closely related to our delta encoding, these methods treat index and value representations separately. In addition to index compression, several studies have explored compressing matrix values. Assuming many elements have similar values, lossless techniques such as delta encoding and run-length encoding have been applied to nonzero elements (Wolfgang et al., 2024; Galanopoulos et al., 2025). However, these approaches primarily target CPUs, and efficient GPU implementations remain challenging. With the growing adoption of mixed-precision algorithms, lossy approaches based on reduced precision have also gained attention. Simplified strategies use lower-precision IEEE 754 formats (FP32 and FP16), while advanced ones employ non-IEEE representations with custom memory accessors (Kawai and Nakajima, 2022; Mukunoki et al., 2023; Graillat et al., 2024b). Some approaches further combine multiple data types within a matrix based on error estimation (Graillat et al., 2024a). While these techniques can effectively reduce memory footprint when accuracy requirements permit, they often require byte-aligned formats, introduce padding to ensure alignment, or impose a constraint on the number of elements processed per memory access. These constraints limit the variety of data representations and GPU efficiency. Overall, previous methods improve SpMV performance for matrices with regular sparsity patterns primarily through data reduction. However, these methods typically treat values and indices separately, imposing alignment constraints on their individual storage and limiting the flexibility in data representation. By contrast, this study proposes a sparse matrix format that jointly compresses matrix values and indices for such matrices. By integrating the design of data representation, index encoding, and memory layout, the proposed format enables efficient GPU execution and flexible data representations, including non-byte-aligned formats.
3
SELL format
This section reviews the SELL format (Monakov et al., 2010), the basis of the proposed PackSELL format. SELL is a state-of-the-art sparse matrix format that groups consecutive rows into slices and aligns the nonzero elements within each slice. This alignment enables a uniform workload and efficient memory access, making SELL well suited for matrices with relatively regular sparsity patterns on SIMD and SIMT architectures. Specifically, as shown in Figure 1, SELL represents a sparse matrix using three arrays, denoted here as val, col, and offset. The arrays val and col store the nonzero values and their corresponding column indices, respectively, in column-major order. However, rows are partitioned into slices of consecutive C rows, and padding values may also be stored so that all rows within a slice have the same number of elements. This padding ensures that all rows in a slice have the same computational workload, which facilitates SIMD and SIMT execution. The starting positions of the slices in val and col are managed in the offset array. When nearby rows have significantly different numbers of nonzero elements, SELL may introduce substantial padding, particularly for large values of C. To mitigate this issue, SELL is often combined with row permutation, where rows are sorted in descending order by the number of nonzero elements per row within blocks of σ rows. This strategy increases the likelihood that rows within each slice have similar numbers of nonzero elements, reducing padding. This technique is known as SELL-C-σ (Kreutzer et al., 2014). Depending on the application, the row permutation can be applied either explicitly or implicitly. If the correctness of the algorithm is unaffected by the row permutation, explicitly reordering the matrix is typically preferable, as it simplifies the implementation and subsequent computations. In this case, the reordered matrix is handled directly 3
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
using the standard SELL format. Conversely, if the original ordering must be preserved, reordering should be applied implicitly during storage and reverted during computation. A straightforward way to implement implicit reordering is to introduce an additional unsigned integer array, perm, that stores permutation information. While several designs are possible, in this study, we store permutation data for every σ rows to save memory usage; specifically, the array stores values in the range 0 to σ − 1 corresponding to the permutation within each block of σ rows. The data type is chosen to be the minimum required to represent σ − 1; for example, an 8-bit unsigned integer suffices when σ ≤ 256. In this implicit case, a simple SpMV algorithm using SELL-C-σ is given below; it employs 0-based indexing: 1: for i = 0, . . . , n − 1 do 2: i′ := ⌊i/σ⌋ · σ + permi 3: k := ⌊i/C⌋ 4: l := i mod C 5: s := offsetk 6: w := (offsetk+1 − s)/C 7: t := 0 8: for j = 0, . . . , w − 1 do 9: p := s + j · C + l 10: t := t + valp · xcolp 11: end for 12: yi′ := t 13: end for Here, a subscripted array name, such as valp , denotes the element at the corresponding index.
4
PackSELL format
SELL-based SpMV performs particularly well for large-scale matrices with relatively regular sparsity patterns, often achieving performance close to the memory bandwidth limit. In such cases, reducing data movement becomes a key direction for further improving performance. In addition, such matrices frequently arise in scientific simulations, where mixed-precision (low-precision) algorithms are actively studied. Motivated by these observations, we propose a new sparse matrix format, PackSELL, built upon the SELL format. Like SELL, PackSELL considers slices of size C and stores nonzero elements with padding so that all rows within each slice have the same number of elements, and in turn, the same workload. To further reduce memory footprint and increase the flexibility of data representation, PackSELL introduces two key techniques: delta encoding for column indices and packing each pair of an index delta and a nonzero value into a single word. Both contribute to reduced memory footprint, and the packing technique provides flexible data representations, covering even non-IEEE 754 formats. 4.1
Delta Encoding for Column Indices
The first key feature is delta encoding of column indices, in which the column index of each nonzero element is represented as the difference from the preceding nonzero element in the same row. Formally, for the i-th row, the column index of a nonzero element ai,j is encoded as a non-negative delta di,j : j − max Si,j if Si,j := {k < j | ai,k ̸= 0} ̸= ∅, di,j = (2) j − di otherwise, where di is an offset used for the leftmost nonzero element, which has no preceding element. While the choice of di should generally ensure that j − di remains small for efficient encoding of the leftmost element, a desirable choice may vary depending on the application and sparsity pattern. When assuming that the matrix has a banded structure, commonly obtained via (reverse) Cuthill–McKee ordering, a natural choice would be i − kleft if kleft < i, di = (3) 0 otherwise, where kleft denotes the lower bandwidth of the matrix. With this definition, di can be computed directly from the row index by storing only kleft . In this work, we employ a slight variation of this definition (see Section 4.3 for details). 4
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
(a) Flag = 0
A P REPRINT
(b) Flag =1 W - 1 bits
0
Delta
V bits
D bits
Value
Delta
Flag
1 Flag
V-bit FP
V-bit integer
… Any V-bit formats
Figure 2: Structure of a word in the PackSELL format. 4.2
Delta–Value Packing
While delta encoding can reduce index size, PackSELL further decreases the memory footprint by packing each delta together with its corresponding matrix value into a single W -bit word. Each word consists of a D-bit unsigned integer for the delta and V bits for the value. By adjusting the allocation of D and V , PackSELL supports a wide range of data formats for matrix values. When required (e.g., in SpMV), the packed delta and value are efficiently unpacked on the fly. 4.2.1
Word Structure
A possible problem of this packing scheme is that some deltas may overflow D bits, especially when D is small. Depending on the sparsity pattern and matrix size, large deltas may arise and force an increase in D to prevent the overflow. This, in turn, reduces V and limits the variety of possible data representations. To address this, PackSELL introduces a flag bit at the least significant bit (LSB) and employs two encoding types depending on the flag value, as illustrated in Figure 2; the bit allocation satisfies W = V + D + 1. When a delta is smaller than 2D , it is stored together with its corresponding matrix value using D and V bits, respectively. Specifically, the word is encoded with flag = 1 as shown in Figure 3(a), which assumes W = 32. The value is converted to a V -bit representation and stored in the upper V bits, while the delta is shifted left by one bit to accommodate the flag in the LSB. The value can be represented in any V -bit data format, including non-IEEE 754 floating-point or integer formats. By contrast, when a delta exceeds 2D − 1, all W − 1(= V + D) bits are used to store the delta (shifted left by one bit as delta << 1) with flag = 0. As a result, the maximum representable delta becomes 2W −1 − 1, equivalent to that of a W -bit signed integer. In this case, however, the corresponding matrix value cannot be represented. This issue is addressed by introducing a dummy element in the sparse matrix format; see Section 4.2.2. This design enables a branch-free unpacking process, as shown in Figure 3(b), which is especially critical for highperformance SpMV on GPUs. First, the flag is extracted to determine the presence of a matrix value. Second, the shift width to recover the delta is adjusted accordingly: V (= 32 − D − 1) bits when flag = 1, and 0 bits when flag = 0. The delta is then extracted by applying left shifting followed by right shifting, based on this shift width. Similarly, the upper V bits corresponding to the value are extracted and multiplied by the flag to drop to zero when the value is absent. Finally, the extracted bits are reinterpreted as the target working datatype. 4.2.2
Examples of Packing and Unpacking
Although arbitrary V -bit data representations can be used in the packing scheme, efficient computation generally requires compatibility with hardware-supported formats. That is, stored values should be easily convertible to such formats. In practice, formats compatible with IEEE 754 floating-point formats are particularly desirable, as they can be efficiently processed by modern arithmetic units. We illustrate two detailed examples of packing and unpacking for W = 32, which are used in Section 5. First, we explain the case where V = 16 and FP16 is used directly. Here, the upper 16 bits of the 32-bit word store an FP16 value. On lines 3–5 in Figure 3(a), after format conversion to FP16, the value is directly placed in the upper bits as uint32_t _value = static_cast < uint32_t >( std :: bit_cast < uint16_t >( static_cast < _half >( value )) ) << 16;
Unpacking is also straightforward. On line 8 in Figure 3(b), the extracted upper 16 bits are reinterpreted as an FP16 value using the CUDA FP16 datatype __half: 5
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
template < typename T > uint32_t packing ( T value , uint32_t delta , uint8_t D ) { uint32_t _value = /* Convert value to a V - bit representation and store it in the upper V bits of uint32_t . */ ; return _value | delta << 1 | 0 x00000001 ; } (a) Packing template < typename T > std :: pair <T , uint32_t > unpacking ( uint32_t pack , uint8_t D ) { uint32_t flag = pack & 0 x00000001 ; uint32_t shift = (31 - D ) * flag ; uint32_t delta = ( pack << shift ) >> ( shift + 1); uint32_t _value = pack & ~((1 u << ( D + 1)) - 1) * flag ; T value = /* Reinterpret _value as type T . */ return { value , delta }; } (b) Unpacking
Figure 3: Pseudocode for the packing and unpacking processes in CUDA/C++. __half value = std :: bit_cast < __half >( static_cast < uint16_t >( pack >> 16));
Next, we consider using an E8MY floating-point format for V bits, which consists of one sign bit, 8 exponent bits, and Y (= 22 − D) mantissa bits. E8MY is compatible with FP32 and preserves key properties such as the exponent bias and implicit leading 1. In this case, to convert FP32 values into E8MY , we use a combination of rounding to the nearest integer and bitwise truncation on lines 3–5 in Figure 3(a): int exp ; frexpf ( value , & exp ); // assuming a float value float scale = ldexpf (1.0 , exp - 24 + D + 1); float __value = std :: round ( value / scale ) * scale ; uint32_t _value = ( std :: bit_cast < uint32_t >( __value ) >> ( D + 1)) << ( D + 1);
Unpacking is simpler than the FP16 case, as the extracted V bits can be directly interpreted as an FP32 value. Line 8 in Figure 3(b) corresponds to float value = std :: bit_cast < float >( _value );
4.3
SELL-style Alignment
Our delta–value packing scheme allows PackSELL to handle large deltas of up to 2W −1 − 1 with flag = 0. However, in this case, the corresponding matrix value cannot be stored together with the delta. To address this, when the delta between two consecutive nonzero elements exceeds the range of the D-bit unsigned integer, PackSELL inserts a dummy element at the same column index as the target element. Since the dummy element has no effective value, it can be encoded with flag = 0. The actual nonzero element then has a delta of 0 relative to the dummy element and can be encoded with flag = 1. This design allows D to be reduced, potentially to as small as 1 bit, at the cost of additional dummy elements. That is, it increases the choice of V for representing matrix values, thereby improving the variation of data representations. After inserting dummy elements for all large deltas, PackSELL applies SELL-style alignment with slice size C to both the original and dummy elements without distinction. Because deltas and values are packed into W -bit words, the val and col arrays in SELL are reduced to a single array, pack, in PackSELL, along with offset for slice boundaries. For this reason, the word size W should be selected considering the memory access efficiency for pack (e.g., 64, 32). 6
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Slice Offsets (offset) a 0 1 b 1 1 c 0 1
0
4 0 d 0 1
5 0 e 0 1
0 6 10 12 Packs of Deltas and Values (pack) a 0 1 c 0 1 b 1 1
4 0
0 d 0 1 …
f 1 1 g 2 1 h 1 1
Flags
…
5 0 f 1 1 e 0 1 g 2 1 h 1 1
Compressed Representation
0
PackSELL Format
Figure 4: Example of the PackSELL format with slice size C = 2 and bit-width for deltas D = 2. For each row, di = 0 for simplicity.
Figure 4 illustrates an example of PackSELL with D = 2 for the same sample matrix as shown in Figure 1, where di is assumed to be 0 for simplicity. In this example, only deltas in the range of 0 to 3 are directly representable. For example, since the delta between entries c and d is 4, a dummy element is inserted immediately before d. In addition, padding elements introduced by SELL-style alignment are treated as zero-valued entries and encoded with flag = 0 and delta = 0, which remains consistent with the delta–value packing scheme. Similar to SELL-C-σ, PackSELL supports row permutation within blocks of σ rows to reduce padding by balancing the number of stored elements across rows within each slice. Most aspects of the row permutation follow those of SELL-C-σ, including the memory layout and the procedure for restoring the original ordering. However, since PackSELL introduces dummy elements for large deltas and applies SELL-style padding after their insertion, the row permutation also accounts for these inserted dummy elements. When using implicit permutation, the offset for the leftmost nonzero elements, di , should be identical within each block of σ rows. Otherwise, to obtain di after the permutation, all di values or the inverse permutation must be stored even when using the definition (3), leading to additional memory accesses. To avoid this, we define di uniformly within each block as: ⌊i/σ⌋ · σ − kleft if kleft < ⌊i/σ⌋ · σ, di = (4) 0 otherwise. 4.4
PackSELL-based SpMV
SpMV in PackSELL can be implemented similarly to the SELL-based implementation in Section 3, with the addition of on-the-fly unpacking. While further hardware-aware optimizations may be possible, a straightforward implementation is shown below and is used in the numerical evaluations: 1: for i = 0, . . . , n − 1 do 2: c := di ▷ Use (4) 3: i′ := ⌊i/σ⌋ · σ + permi 4: k := ⌊i/C⌋ 5: l := i mod C 6: s := offsetk 7: w := (offsetk+1 − s)/C 8: t := 0 9: for j = 0, . . . , w − 1 do 10: p := s + j · C + l 11: Unpack packp to obtain (v, d) 12: c := c + d 13: t := t + v · xc 14: end for 15: yi ′ = t 16: end for Efficient unpacking is critical for high performance. In this study, we implement it as explained in Section 4.2. 7
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
5
A P REPRINT
Numerical Experiments
We evaluate PackSELL from two perspectives: the performance of standalone SpMV kernels (Section 5.1) and their performance when integrated into Krylov subspace methods for solving sparse linear systems (Section 5.2). Although PackSELL supports multiple word sizes W (e.g., 16, 32, and 64), we focus on W = 32. This setting is comparable to common low-precision SpMV implementations that use FP32 or FP16 values with 32-bit (signed) integer indices. All experiments were conducted on a GPU node of the Gardenia supercomputer at Kyoto University. Although each node has four NVIDIA A100 (80GB SXM) GPUs, we used a single GPU. The A100 GPU provides a memory bandwidth of 2,039 GB/s, and the installed driver version was 580.105.08. The evaluated kernels and solvers were developed in CUDA and C++ using the NVIDIA CUDA Compiler (version 13.0.2) and the GNU C++ Compiler (version 13.3.0), with the options -O3, -sm_80, -std=c++20, and -Xcompiler="-O3 -fopenmp -std=c++20". The C++20 standard was required to support advanced bit manipulation. We also used the cuSPARSE library from the same CUDA Toolkit. 5.1
Evaluation via Standalone SpMV Kernels
This subsection evaluates standalone SpMV kernels. To reflect typical use cases of GPUs for low-precision SpMV, we used relatively large matrices from the SuiteSparse Matrix Collection (Davis and Hu, 2011). Specifically, we selected most large-scale real (and integer) matrices with at least 65,536 (= 216 ) columns. Since precision is a key focus of this study, we excluded binary matrices having only 0 and 1 values. Very small matrices were also excluded, as they are inefficient on modern GPUs and their column indices can be represented with 16-bit or 8-bit integers. Among the 438 matrices satisfying these conditions, three matrices (MOLIERE_2016, GAP-urand, and GAP-kron) were excluded due to memory limitations. Consequently, 435 matrices were used in the evaluation. All performance results in this subsection are averages of 10,000 runs after 100 warm-up runs. In addition, FLOPS is measured assuming two floating-point operations per nonzero element, excluding padding and dummy elements. 5.1.1
FP16 SpMV
First, we examine PackSELL for FP16 SpMV, where the input and output vectors are also stored in FP16. Based on the algorithm in Section 4.4, we developed an FP16 SpMV kernel in PackSELL with W = 32 and D = 15; FP16 values were directly embedded into the remaining V (= 16) bits, as explained in Section 4.2.2. The kernel was parallelized such that each thread processes one row, with 256 CUDA threads per block. This kernel was compared with five FP16 SpMV kernels: DASP (Lu and Liu, 2023) and four NVIDIA cuSPARSE kernels1 for COO, CSR, SELL, and BSR formats, hereafter denoted as cuCOO, cuCSR, cuSELL, and cuBSR, respectively. All cuSPARSE kernels were used with cusparseSpMV_preprocess() for potential optimization. BSR used either 2×2 or 4×4 blocks, whichever performed better. For PackSELL and cuSELL, we fixed C = 32 (warp size) and σ = 256, and explicitly reordered matrix rows, since cuSELL, the closest counterpart to PackSELL, does not support implicit permutation. Reordering the output vector to restore the original ordering would incur additional memory accesses and penalize cuSELL. Figure 5 summarizes the achieved FLOPS of each kernel. The horizontal axis shows the relative standard deviation of the number of nonzero elements per row (RSD), where smaller values indicate a more uniform distribution of nonzero elements across all rows. The upper plot shows results for individual matrices, and the lower plot reports the 75th percentile for seven bins (RSD = 0 or in [10e , 10e+1 ) for e ∈ [−2, 3]) to highlight representative performance. To further illustrate the relationship between matrix characteristics and performance, Figure 6 presents detailed results for 11 matrices listed in Table 1, selected to represent diverse sparsity structures. These results first reconfirm the strengths and limitations of SELL-based formats. In general, SELL is suitable for matrices with regular sparsity patterns (i.e., small RSD). For such matrices, cuSELL approached the roofline peak more closely than other kernels, such as cuCSR and cuCOO. By contrast, for matrices with high RSD values, such as Ga41As41H72, language, and degme, SELL incurs substantial padding and underperforms CSR and COO, which require no padding. Accordingly, cuSELL and PackSELL underperformed cuCSR and cuCOO in these cases, which indicates a limitation of SELL-based approaches. Achieving higher performance for such matrices remains future work. Nevertheless, among the 435 matrices tested, cuSELL outperformed cuCSR for 247 matrices; this number increases to 292 when including cases where cuSELL achieved at least 90% of the performance of cuCSR. These results indicate that improving upon cuSELL is practically meaningful. We refer to these 292 matrices as SELL-suitable matrices and focus on them in the following. 1
Although cuSPARSE provides an interface for computing y = αAx + βy for scalars α and β, we observed better performance when setting β = 0.
8
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
GFLOPS
800
cuCSR
cuSELL
A P REPRINT
PackSELL
600 400 200 0
GFLOPS
800
cuBSR cuCOO
600
cuCSR cuSELL
DASP PackSELL
400 200 0 0
10−1
10−2
101
100
103
102
104
Relative standard deviation of the number of nonzeros per row (RSD)
Figure 5: Achieved FLOPS for six SpMV kernels. Dotted horizontal lines indicate the upper bound based only on nonzero elements (FP16 values and 32-bit column indices).
cuBSR
cuCOO
cuCSR
cuSELL
DASP
PackSELL
GFLOPS
750 500 250 0 mri2
para... Curl... GL7d... Flan... cont... HV15... Spie... Ga41... lang... degm...
Figure 6: Detailed results for matrices listed in Table 1.
For many SELL-suitable matrices, PackSELL outperformed cuSELL and exceeded the upper bound of cuCSR and cuSELL by reducing the memory footprint. In particular, PackSELL achieved larger gains for matrices with small RSD and many nonzero elements, such as CurlCurl_4, Flan_1565, and HV15R. Figure 7 further explains these results showing the memory usage of PackSELL relative to SELL. In most cases, PackSELL reduced the memory footprint compared to SELL via delta–value packing, with few dummy elements. In particular, when many nonzero elements are densely clustered, as in the three examples above, D-bit (15-bit) unsigned integers can represent most deltas, enabling compression rates close to the lower bound of 0.75 (= 32 bits / 48 bits). By contrast, when nonzero elements were widely scattered, the number of dummy elements increased, and the performance gain decreased. In such cases, except for small matrices like mri2, where D bits still covered most deltas, PackSELL required a similar amount of storage to cuSELL (e.g., cont11_l and GL7d17), and even underperformed cuSELL due to the increased storage (e.g., parabolic_fem). These results suggest that matrix reordering to improve the locality of nonzero elements is promising for further improvements of PackSELL, which we leave for future work. Figure 8 shows the speedups of PackSELL for the SELL-suitable matrices. PackSELL outperformed cuSELL on 201 of 292 matrices. For matrices with many nonzero elements, it achieved consistent speedups of around 1.5×, matching the ideal gain expected from the reduced data size. In some cases, speedups reached up to 1.63×, likely due to improved memory access patterns from delta–value packing, which may enhance cache efficiency by using a single array for values and deltas. On smaller matrices where even cuSELL struggled to fully utilize memory bandwidth, PackSELL underperformed cuSELL; addressing these cases with hardware-oriented optimizations, such as those to effectively hide memory latency, remains future work. Compared with cuCSR and DASP, PackSELL outperformed on 274 and 265 matrices, respectively, with maximum speedups of 3.01× and 5.09×. For reference, the maximum speedups of PackSELL against cuCOO and cuBSR were 2.76× and 7.57×, respectively. 9
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Table 1: Selected Matrices Matrix / Application mri2 / Computer graphics parabolic_fem / CFD
nnz †
RSD
569,160
0
3,674,625
0.0218
n×m
Sparsity
63,240 ×147,456 525,825 ×525,825
CurlCurl_4 / Model reduction
2,380,515 ×2,380,515
26,515,867
0.0768
GL7d17 / Differential matrix
1,548,650 ×955,128
25,978,098
0.1180
Flan_1565 / Structural analysis
1,564,794 114,165,372 ×1,564,794
0.1550
cont11_l / Linear programming
1,468,599 ×1,961,394
5,382,999
0.2571
HV15R / CFD
2,017,169 283,073,458 ×2,017,169
0.3845
Spielman_k600 / Chimera Laplacian
72,180,402 216,902,404 ×72,180,402
0.6655
268,096 ×268,096
Ga41As41H72 / Quantum chemistry
18,488,476
1.5282
Memory footprint ratio
language 399,130 1,216,334 6.7958 / Weighted graph ×399,130 185,501 degme 8,127,528 33.0696 ×659,415 / Linear programming † nnz denotes the number of nonzero elements. 1.2
Baseline Lower bound
1.0
parabolic_fem cont11_l
0.8
language
Spielman_k600
CurlCurl_4
mri2
degme
106
107
0.6
GL7d17
HV15R
Ga41As41H72
105
Flan_1565
108
109
Number of nonzero elements
Figure 7: Memory footprint ratio of PackSELL to SELL, shown as scatter and letter-value plots. Orange markers indicate matrices listed in Table 1.
5.1.2
SpMV Using Non-IEEE Formats
We next assess the advantage of PackSELL in supporting various data formats. We tested PackSELL-based SpMV while varying D from 1 to 12. The implementation in each case was mostly the same as that in Section 5.1.1, expect that matrix values were represented in E8MY with Y = 22 − D. We compared these kernels against FP32 cuSELL; in both cases, the input and output vectors were stored in FP32. In addition to performance, we evaluated accuracy using the backward error, defined as ∥y − Ax∥ , ∥A∥∥x∥
(5)
where the infinity norm was used. Although acceptable values of the backward error depend on the application, comparing backward errors across different precision settings provides a relative measure of practicality. For more practical evaluations using iterative linear solvers, see the next subsection. 10
Speedup over DASP
Speedup over cuCSR
Speedup over cuSELL
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
1.5 1.0 0.5 3.0 2.0 1.0
4.0 2.0 1.0 105
107
106
109
108
Number of nonzero elements
Speedup over FP32 cuSELL
Figure 8: Speedups of PackSELL over cuSELL, cuCSR, and DASP for the SELL-suitable matrices. 2.0 1.5 1.0 0.5
10−4 PackSELL FP32 cuSELL
10−5 10
−6
1
2
3
4
5
6
7
8
9
10
11
FP16 cuSELL BF16 cuSELL 12
13
14
BF16 cuSELL
10−3
FP16 cuSELL
Backward error
0.0
15
Number of bits used for deltas D (E8MY, Y = 22 - D)
Figure 9: Achieved performance and backward error of cuSELL and PackSELL-based SpMV using E8MY . The reported backward error is the average across all matrices, excluding the result of FP16 cuSELL for t3dh_e, which is omitted due to an outlier caused by underflow. For a focused evaluation, we restricted our analysis to the 292 SELL-suitable matrices. To mitigate overflow and underflow, we applied diagonal scaling G−1 A, where G = diag(g1 , . . . , gn ) with gi = Σj |ai,j |. The performance of PackSELL relative to FP32 cuSELL is summarized in Figure 9 using letter-value plots, along with the average backward error over all matrices. For reference, the figure also includes results for cuSELL in FP16 and bfloat16 (BF16) with FP16 or BF16 inputs and FP32 outputs. As expected, increasing D (and thus decreasing Y ) improved performance by reducing the number of dummy elements for large deltas, at the cost of lower precision. Although the backward error gradually increased as the number of mantissa bits Y decreased, the degradation was moderate, especially when only a few bits were truncated; in most cases, the increase remained within one order of magnitude. In contrast to the comparable backward error, PackSELL with small D outperformed cuSELL in many cases, e.g., 129 out of 292 matrices for D = 1 (E8M21) and 145 matrices for D = 2 (E8M20). It also delivered notable speedups, approaching 2.0× in some instances. Furthermore, under these settings, PackSELL even outperformed FP16/BF16 cuSELL, while achieving much lower backward errors than both. As observed in Section 5.1.1, the performance of PackSELL is affected by the locality of nonzero elements. These results suggest that, for matrices suited to PackSELL 11
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Table 2: Test Matrices: 15 Symmetric Positive Definite Matrices (Top) and 15 Nonsymmetric Matrices (Bottom) Matrix n nnz nnz /n κ2 ‡ Bump_2911 2,911,419 127,729,899 43.87 1.5e+06 Emilia_923 923,136 40,373,538 43.74 1.8e+06 G3_circuit 1,585,478 7,660,826 4.83 2.4e+05 Queen_4147 4,147,110 316,548,962 76.33 9.4e+06 Serena 1,391,349 64,131,971 46.09 2.7e+04 apache2 715,176 4,817,870 6.74 1.3e+06 audikw_1 943,695 77,651,847 82.28 1.4e+07 ecology2 999,999 4,995,991 5.00 6.3e+07 ldoor 952,203 42,493,817 44.63 4.1e+06 thermal2 1,228,045 8,580,313 6.99 4.5e+06 tmt_sym 726,713 5,080,961 6.99 2.5e+08 HPCG_7_7_7 2,097,152 55,742,968 26.58 2.3e+03 HPCG_8_7_7 4,194,304 111,777,784 26.65 3.0e+03 HPCG_8_8_7 8,388,608 224,140,792 26.72 4.5e+03 HPCG_8_8_8 16,777,216 449,455,096 26.79 8.9e+03 Freescale1 3,428,755 17,052,626 4.97 3.8e+07 Transport 1,602,111 23,487,281 14.66 4.0e+05 atmosmodd 1,270,432 8,814,880 6.94 5.2e+03 atmosmodl 1,489,752 10,319,760 6.93 1.1e+03 rajat31 4,690,002 20,316,253 4.33 7.9e+03 ss 1,652,680 34,753,577 21.03 5.0e+04 stokes 11,449,533 349,321,980 30.51 5.1e+06 t2em 921,632 4,590,832 4.98 2.3e+05 tmt_unsym 917,825 4,584,801 5.00 1.5e+09 vas_stokes_1M 1,090,664 34,767,207 31.88 4.3e+07 vas_stokes_2M 2,146,677 65,129,037 30.34 6.3e+07 HPGMP_7_7_7 2,097,152 55,742,968 26.58 1.6e+03 HPGMP_8_7_7 4,194,304 111,777,784 26.65 1.9e+03 HPGMP_8_8_7 8,388,608 224,140,792 26.72 2.2e+03 HPGMP_8_8_8 16,777,216 449,455,096 26.79 4.2e+03 ‡ κ2 denotes an estimated 2-norm condition number of Ḡ−1 AḠ−1 .
with high locality of nonzero elements, even a few D bits suffice to represent deltas and yield both higher accuracy and improved performance. Another notable case is D = 12 (E8M10). Despite allocating 8 bits to the exponent (and using an FP32 input vector), PackSELL outperformed FP16 cuSELL for 175 matrices. For matrices with a wide dynamic range, a larger exponent part can be beneficial. Indeed, FP16 cuSELL suffered severe underflow for the t3dh_e matrix, with a backward error of 8.6 × 10−1 , even with diagonal scaling (this case is excluded from Figure 9 as an outlier). In those cases, BF16 might act as an alternative; however, the performance of BF16 cuSELL is similar to that of FP16 cuSELL, and its 7-bit mantissa can be insufficient in certain applications. These results further emphasize the key benefit of PackSELL: the ability to flexibly adjust exponent and mantissa bits according to the characteristics of matrices and applications. 5.2
Evaluation via Krylov Solvers
This subsection evaluates the practical performance of PackSELL-based SpMV kernels in Krylov subspace methods for solving sparse linear systems Ax = b. Specifically, we consider two mixed-precision solvers: the state-of-the-art FP16-enabled solver F3R (Suzuki and Iwashita, 2025) and an inner-outer conjugate gradient (IO-CG) method. The experiments used 30 square matrices listed in Table 2. Of these, 22 belong to the previously described set of 292 SELL-suitable matrices, while the remaining 8 are derived from the HPCG (Dongarra et al., 2016) and HPGMxP (Yamazaki et al., 2022b) benchmarks. These matrices are denoted as HPCG_x_y_z and HPGMP_x_y_z, where the number of rows n is given by 2x+y+z . For HPGMxP, the parameter to control the asymmetry was set to 0.5. While F3R experiments used all 30 matrices, IO-CG experiments used only the 15 symmetric positive definite (SPD) matrices since CG is designed for SPD systems. For all matrices, diagonal scaling of the form Ḡ−1 AḠ−1 was applied, where p Ḡ = diag(ḡ1 , . . . , ḡn ) with ḡi = |ai,i |. All numerical results below are averages over five runs. In each run, the right-hand side vector b was a random vector with elements uniformly distributed in the range [0, 1), and the initial guess was the zero vector. Convergence was 12
4
A P REPRINT
FP16-F3R PackSELL-F3R
3 2 1
1.3
RE S
1.2
G
16
1.1
Speedup over FP16-F3R
M
1.0
Pa -F c 3R -F kSE 3R LL
0
FP
Speedup over FP64-F3R
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
Figure 10: Performance of different F3R implementations. determined using the following criterion, where x̃ denotes an approximate solution: ∥b − Ax̃∥2 /∥b∥2 < 10−9 . 5.2.1
(6)
FP16-Enabled F3R Solver
F3R is a nested Krylov subspace method that hierarchically combines Krylov solvers as preconditioners for other Krylov iterations. It consists of four layers: three flexible GMRES (FGMRES) layers and an innermost preconditioned Richardson layer. Within this hierarchy, the innermost Richardson and its outer FGMRES employ FP16 SpMV. Since FP16 SpMV accounts for over 85% of all SpMV operations under the default parameter settings, implementing these kernels with PackSELL, as described in Section 5.1.1, provides an evaluation of practical performance. We implemented F3R in three ways. The first is an FP64-only version (FP64-F3R), and the second is a mixed-precision variant that uses FP16 SpMV in the two inner layers mentioned above (FP16-F3R). For their SpMV operations, we used our SELL-based kernels following the algorithm in Section 3, with C = 32 and σ = 256. This is because mixed-precision F3R requires combinations of vector and matrix value types unsupported by cuSPARSE. Additionally, we adopted the implicit σ permutation to preserve the original ordering, which is also not available in cuSPARSE. The third implementation has the same structure as FP16-F3R but uses the FP16 SpMV kernel in PackSELL (V = 16, D = 15) described in Section 5.1.1 (PackSELL-F3R). This implementation also employs implicit permutation with σ = 256. In all solvers, all other parameters followed the default settings; for example, the innermost Richardson used the approximate inverse preconditioner SD-AINV (Suzuki et al., 2022); see (Suzuki and Iwashita, 2025; Suzuki, 2025) for details. The performance of these solvers is presented in Figure 10. The left plot shows the speedups of PackSELL-F3R over FP16-F3R and FP64-F3R, with dotted lines connecting results for each problem. The right plot summarizes the speedup over FP64-F3R, with reference results for an FP64 GMRES solver combined with SD-AINV and restarted every 100 iterations. Except for one case (rajat31), F3R outperformed GMRES, and PackSELL further improved performance. Since FP16 values are directly embedded in PackSELL, FP16-F3R and PackSELL-F3R exhibit identical convergence. Thus, the performance benefits in SpMV translate into gains in the overall performance. For all problems except Freescale1, PackSELL-F3R outperformed FP16-F3R by 1.06–1.32×. These improvements resulted in speedups of up to 2.64× over FP64-F3R and up to 28.14× over GMRES, indicating that PackSELL is effective in this practical scenario. 5.2.2
Inner-Outer CG Solver
We demonstrate that SpMV with custom data formats enabled by PackSELL increases the flexibility of the design of mixed-precision solvers. To illustrate this, we consider IO-CG, a variant of the mixed-precision solver presented in (Buttari et al., 2008). In IO-CG, min iterations of preconditioned CG (PCG) serve as a preconditioner for the outer flexible CG (FCG) (Notay, 2000). While this scheme is relatively uncommon, it can serve as a good example of how PackSELL can increase the flexibility of solver design. We developed four IO-CG variants: one FP64-only solver and three mixed-precision solvers. In all mixed-precision variants, the outer FCG uses FP64, while the inner PCG primarily employs FP32. These three variants differ in the SpMV kernel used for the coefficient matrix A: FP32 SELL, FP16 SELL, and PackSELL, in which the values are encoded in E8MY (Y = 22 − D) with D ranging from 1 to 12. We refer to these solvers as FP64-, FP32-, FP16-, 13
min = 50
min = 20
2.0
A P REPRINT
min = 80
1.5 1.0 0.5
G
G
-C IO
M
Y-
G
-C -IO
16
E8
G
-C FP
-C 32 FP
-IO 64
FP
-IO
G
G
-C IO
M
Y-
G
-C -IO
16
E8
G
-C FP
-C
-IO 32
64 FP
FP
-IO
-C
G
G M
Y-
IO
-C
G 16
-IO
E8
-C
-C -IO
32
64
FP
FP
-IO
G
0.0
FP
Speedup over FP64 PCG
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
Figure 11: Performance of four mixed-precision inner-outer CG variants relative to the standard FP64 PCG solver. Table 3: Formats That Achieved the Best Performance Matrix Bump_2911 Emilia_923 G3_circuit Queen_4147 Serena apache2 audikw_1 ecology2 ldoor thermal2 tmt_sym HPCG_7_7_7 HPCG_8_7_7 HPCG_8_8_7 HPCG_8_8_8
min = 20 E8M11 E8M12 E8M14 E8M13 E8M13 E8M20 E8M13 E8M11 E8M13 E8M11 E8M10 E8M10 E8M10 E8M11 E8M10
min = 50 E8M14 E8M15 E8M13 E8M14 E8M12 E8M15 E8M14 E8M12 E8M20 E8M16 E8M10 E8M11 E8M11 E8M13 E8M12
min = 80 E8M15 E8M19 E8M17 E8M14 E8M11 E8M21 E8M16 E8M12 E8M17 E8M16 E8M20 E8M14 E8M14 E8M13 E8M14
and E8MY -IO-CG, respectively. SpMV (and preconditioning) dominate the cost of PCG, and the cost of index access is also non-negligible when storing A explicitly. Thus, based primarily on SpMV performance, the expected ideal speedups over FP64 PCG for this approach and similar previous methods are roughly 1.5× and 2.0× when using FP32 and FP16 for SpMV, respectively. We applied the SD-AINV preconditioner to all solvers, and tested three settings for the number of inner iterations min = 20, 50, and 80. Figure 11 shows the achieved performance of the four IO-CG variants relative to the standard FP64 PCG solver; for E8MY -IO-CG, it reports the best results obtained by selecting the optimal E8MY format listed in Table 3. First, FP64-IO-CG generally underperformed FP64 PCG due to the increased iterations and operations introduced by the inner-outer scheme. However, this overhead became less significant as min increased; for min = 50 and 80, FP64-IOCG exhibited convergence behavior and performance close to PCG. Due to these properties, FP32-IO-CG outperformed PCG in many cases by reducing data movement. Nevertheless, its speedup over PCG remained around 1.3×, which is not particularly remarkable compared to the ideal gain and the speedups achieved in previous studies (Buttari et al., 2008; Yamazaki et al., 2022a; Guo et al., 2025; Chen et al., 2026). One approach to further improve solver performance is to use FP16. However, as shown in Figure 12, larger min also degraded convergence, particularly in FP16-IO-CG, because the precision became insufficient during many inner iterations. Even when FP32-IO-CG behaved similarly to FP64-IO-CG, FP16-IO-CG required more iterations (e.g., for ldoor); this increase outweighed the benefits of reduced data movement, resulting in performance similar to or inferior to FP32-IO-CG. A well-known non-IEEE alternative is BF16 (E8M7). In our experiments, however, E8M10 showed convergence similar to FP16 (Figure 12). This indicates that the number of mantissa bits, rather than exponent bits, is critical for good convergence and that BF16 is also ineffective in this context. By contrast, as discussed in Section 5.1.2, PackSELL can control the mantissa length of E8MY , achieving FP32level accuracy while outperforming both FP32 and FP16 SpMV kernels. Owing to these advantages, E8MY -IO-CG 14
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
Residual norm
PCG FP64-IO-CG 100
FP32-IO-CG FP16-IO-CG
E8M14-IO-CG E8M10-IO-CG
min = 50
min = 20
A P REPRINT
min = 80
10−3 10−6 10−9 0
500
0
Iteration count
500
0
Iteration count
500
Iteration count
(a) HPCG_8_8_8
Residual norm
PCG FP64-IO-CG 100
FP32-IO-CG FP16-IO-CG
E8M14-IO-CG E8M10-IO-CG
min = 50
min = 20
min = 80
10−3 10−6 10−9 0
2500
5000
Iteration count
0
2500
5000
Iteration count
0
2500
5000
Iteration count
(b) ldoor
Figure 12: History of the relative residual norm. For IO-CG, the iteration count denotes the number of inner iterations.
outperformed FP64 PCG on 13 of the 15 problems, even for min = 80. Indeed, Figure 12 shows that its convergence behavior closely matched that of FP32-IO-CG, for instance at Y = 14. Consequently, unlike FP32-IO-CG and FP16-IO-CG, E8MY -IO-CG attained speedups of up to 2.0× relative to FP64 PCG, which is a notable improvement over previous techniques and approaching the ideal gain. Furthermore, Table 3 shows that the best format tended to require more mantissa bits as min increased, in order to maintain sufficient accuracy over many inner iterations. That is, the ability of PackSELL to finely tune the mantissa length benefits the inner-outer scheme in IO-CG. Admittedly, IO-CG requires more memory than PCG, and several practical considerations remain. However, these results demonstrate that PackSELL can realize solvers that outperform PCG with gains close to the ideal, which are difficult to achieve when relying on standard formats such as FP16 and BF16 in this IO-CG scheme. Therefore, PackSELL and its core idea of packing value and index information to enable flexible data representation are expected to be a useful basis for future research on advanced mixed-precision solvers, including more practical variants of IO-CG.
6
Conclusions
We propose a new sparse matrix format, PackSELL, based on SELL, for high-performance SpMV on GPUs across various data representations. It incorporates delta encoding of column indices and a new delta–value packing scheme to reduce memory footprint and data movement. This packing scheme also allows flexible data representations by adjusting the bit allocation between deltas and values. The numerical results on an NVIDIA A100 GPU show that SELL is effective for 292 of the 435 practical matrices tested, indicating that improving SELL-based SpMV is meaningful. Although hardware-oriented optimizations particularly for smaller matrices remain future work, PackSELL performs well across a variety of data representations, including custom non-IEEE formats. For example, on 201 of these 292 matrices, it outperforms the cuSPARSE SELL kernel in FP16 SpMV by up to 1.63×. Furthermore, evaluations with Krylov subspace methods confirm its practical effectiveness. In the FP16-enabled F3R solver, using PackSELL for FP16 SpMV improves the overall performance by up to 1.32×. PackSELL is also effective when FP32 is costly and FP16 lacks sufficient accuracy, due to its ability to support custom data formats. An inner15
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
outer CG method using PackSELL with the E8MY format outperforms FP64 PCG by up to 2.09×. Mixed-precision solvers using non-IEEE formats remain relatively underexplored because of the limited support of these formats in high-performance SpMV kernels. Thus, PackSELL not only improves FP16-enabled solvers, but also provides a practical means to further investigate advanced mixed-precision solvers that utilize custom data representations. Future work includes three directions. First, evaluating PackSELL-based SpMV on other hardware platforms may reveal different performance trends. While a similar behavior is expected in memory-bound scenarios, other hardware features such as cache size may change the importance of hardware-oriented optimizations. Second, extending the delta–value packing scheme to CSR-based formats may address the cases where SELL needs excessive padding, though this may require careful handling of load balancing and inter-thread communication. Third, applying PackSELL to other sparse matrix kernels, such as sparse triangular solves, is promising because some of their implementations are similar to SpMV.
Acknowledgment This work was supported by JSPS KAKENHI Grant Number JP25K24388.
References Ahmad Abdelfattah, Hartwig Anzt, Erik G Boman, Erin Carson, Terry Cojean, Jack Dongarra, Alyson Fox, Mark Gates, Nicholas J Higham, Xiaoye S Li, Jennifer Loe, Piotr Luszczek, Srikara Pranesh, Siva Rajamanickam, Tobias Ribizel, Barry F Smith, Kasia Swirydowicz, Stephen Thomas, Stanimire Tomov, Yaohung M Tsai, and Ulrike Meier Yang. 2021. A Survey of Numerical Linear Algebra Methods Utilizing Mixed-Precision Arithmetic. The International Journal of High Performance Computing Applications 35, 4 (July 2021), 344–369. doi:10.1177/ 10943420211003313 José I. Aliaga, Hartwig Anzt, Thomas Grützmacher, Enrique S. Quintana-Ortí, and Andrés E. Tomás. 2022. Compression and Load Balancing for Efficient Sparse Matrix-vector Product on Multicore Processors and Graphics Processing Units. Concurrency and Computation 34, 14 (June 2022), e6515. doi:10.1002/cpe.6515 Patrick Amestoy, Alfredo Buttari, Nicholas J. Higham, Jean-Yves L’Excellent, Theo Mary, and Bastien Vieublé. 2024. Five-Precision GMRES-Based Iterative Refinement. SIAM J. Matrix Anal. Appl. 45, 1 (March 2024), 529–552. doi:10.1137/23M1549079 Andrew Anderson and David Gregg. 2016. Vectorization of Multibyte Floating Point Data Formats. In Proc. 2016 Int. Conf. Parallel Archit. Compil. (PACT ’16). Association for Computing Machinery, New York, NY, USA, 363–372. doi:10.1145/2967938.2967966 Hartwig Anzt, Terry Cojean, Chen Yen-Chen, Jack Dongarra, Goran Flegar, Pratik Nayak, Stanimire Tomov, Yuhsiang M. Tsai, and Weichung Wang. 2020. Load-Balancing Sparse Matrix Vector Product Kernels on GPUs. ACM Trans. Parallel Comput. 7, 1 (March 2020), 1–26. doi:10.1145/3380930 Hartwig Anzt, Stanimire Tomov, and Jack Dongarra. 2014. Implementing a Sparse Matrix Vector Product for the SELL-C/SELL-C-σ Formats on NVIDIA GPUs. Technical Report. University of Tennessee. Arash Ashari, Naser Sedaghati, John Eisenlohr, Srinivasan Parthasarath, and P. Sadayappan. 2014. Fast Sparse MatrixVector Multiplication on GPUs for Graph Applications. In SC14 Int. Conf. High Perform. Comput. Netw. Storage Anal. IEEE, New Orleans, LA, USA, 781–792. doi:10.1109/SC.2014.69 Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, and Henk van der Vorst. 1994. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, Philadelphia, PA. doi:10.1137/1.9781611971538 Nathan Bell and Michael Garland. 2009. Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors. In Proc. Conf. High Perform. Comput. Netw. Storage Anal. (SC ’09). Association for Computing Machinery, New York, NY, USA, 1–11. doi:10.1145/1654059.1654078 Alfredo Buttari, Jack Dongarra, Jakub Kurzak, Piotr Luszczek, and Stanimir Tomov. 2008. Using Mixed Precision for Sparse Matrix Computations to Enhance the Performance While Achieving 64-Bit Accuracy. ACM Trans. Math. Softw. 34, 4 (July 2008), 1–22. doi:10.1145/1377596.1377597 Yanxiang Chen, Pablo De Oliveira Castro, Paolo Bientinesi, Niclas Jansson, and Roman Iakymchuk. 2026. Enabling Mixed-Precision in Spectral Element Codes. Future Generation Computer Systems 174 (Jan. 2026), 107990. doi:10.1016/j.future.2025.107990 16
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Xing Cong, FuKai Sun, YiFan Chen, Chenhao Xie, Yi Liu, and Depei Qian. 2025. CB-SpMV:A Data Aggregating and Balance Algorithm for for Cache-Friendly Block-Based SpMV on GPUs. In Proc. 39th ACM Int. Conf. Supercomput. ACM, Salt Lake City USA, 149–160. doi:10.1145/3721145.3725746 Timothy A Davis and Yifan Hu. 2011. The University of Florida Sparse Matrix Collection. ACM Trans. Math. Softw. 38, 1 (2011), 1–25. doi:10.1145/2049662.2049663 Jack Dongarra, Michael A Heroux, and Piotr Luszczek. 2016. High-Performance Conjugate-Gradient Benchmark: A New Metric for Ranking High-Performance Computing Systems. The International Journal of High Performance Computing Applications 30, 1 (Feb. 2016), 3–10. doi:10.1177/1094342015593158 Salvatore Filippone, Valeria Cardellini, Davide Barbieri, and Alessandro Fanfarillo. 2017. Sparse Matrix-Vector Multiplication on GPGPUs. ACM Trans. Math. Softw. 43, 4 (Jan. 2017), 30:1–30:49. doi:10.1145/3017994 Dimitrios Galanopoulos, Panagiotis Mpakos, Petros Anastasiadis, Nectarios Koziris, and Georgios Goumas. 2025. DIV: An Index & Value Compression Method for SpMV on Large Matrices. In Proc. 39th ACM Int. Conf. Supercomput. ACM, Salt Lake City USA, 705–717. doi:10.1145/3721145.3725767 Jianhua Gao, Bingjie Liu, Weixing Ji, and Hua Huang. 2024. A Systematic Literature Survey of Sparse Matrix-Vector Multiplication. arXiv:2404.06047 [cs] doi:10.48550/arXiv.2404.06047 Stef Graillat, Fabienne Jézéquel, Théo Mary, and Roméo Molina. 2024a. Adaptive Precision Sparse Matrix–Vector Product and Its Application to Krylov Solvers. SIAM J. Sci. Comput. 46, 1 (2024), C30–C56. doi:10.1137/ 22M1522619 Stef Graillat, Fabienne Jézéquel, Theo Mary, Roméo Molina, and Daichi Mukunoki. 2024b. Reduced-Precision and Reduced-Exponent Formats for Accelerating Adaptive Precision Sparse Matrix–Vector Product. In Euro-Par 2024: Parallel Processing, Jesus Carretero, Sameer Shende, Javier Garcia-Blas, Ivona Brandic, Katzalin Olcoz, and Martin Schreiber (Eds.). Vol. 14803. Springer Nature Switzerland, Cham, 17–30. doi:10.1007/978-3-031-69583-4_2 Yichen Guo, Eric de Sturler, and Tim Warburton. 2025. An Adaptive Mixed Precision and Dynamically Scaled Preconditioned Conjugate Gradient Algorithm. arXiv:2505.04155 [math] doi:10.48550/arXiv.2505.04155 Laslo Hunhold and James Quinlan. 2025. Evaluation of Bfloat16, Posit, and Takum Arithmetics in Sparse Linear Solvers. In 2025 IEEE 32nd Symp. Comput. Arith. ARITH. IEEE, 61–68. doi:10.1109/ARITH64983.2025.00019 Tsuyoshi Ichimura, Kohei Fujita, Takuma Yamaguchi, Akira Naruse, Jack C. Wells, Thomas C. Schulthess, Tjerk P. Straatsma, Christopher J. Zimmer, Maxime Martinasso, Kengo Nakajima, Muneo Hori, and Lalith Maddegedara. 2018. A Fast Scalable Implicit Solver for Nonlinear Time-Evolution Earthquake City Problem on Low-Ordered Unstructured Finite Elements with Artificial Intelligence and Transprecision Computing. In SC18 Int. Conf. High Perform. Comput. Netw. Storage Anal. IEEE, Dallas, TX, USA, 627–637. doi:10.1109/SC.2018.00052 Soichiro Ikuno, Yuki Kawaguchi, Norihisa Fujita, Taku Itoh, Susumu Nakata, and Kota Watanabe. 2012. Iterative Solver for Linear System Obtained by Edge Element: Variable Preconditioned Method With Mixed Precision on GPU. IEEE Trans. Magn. 48, 2 (Feb. 2012), 467–470. doi:10.1109/TMAG.2011.2175375 Takeshi Iwashita, Kengo Suzuki, and Takeshi Fukaya. 2020. An Integer Arithmetic-Based Sparse Linear Solver Using a GMRES Method and Iterative Refinement. In 2020 IEEEACM 11th Workshop Latest Adv. Scalable Algorithms Large-Scale Syst. ScalA. IEEE, 1–8. doi:10.1109/ScalA51936.2020.00006 Juan Luis Jerez, George A. Constantinides, and Eric C. Kerrigan. 2015. A Low Complexity Scaling Method for the Lanczos Kernel in Fixed-Point Arithmetic. IEEE Trans. Comput. 64, 2 (2015), 303–315. doi:10.1109/TC.2013.162 Masatoshi Kawai and Kengo Nakajima. 2022. Low/Adaptive Precision Computation in Preconditioned Iterative Solvers for Ill-Conditioned Problems. In Int. Conf. High Perform. Comput. Asia-Pac. Reg. (HPCAsia ’22). Association for Computing Machinery, New York, NY, USA, 30–40. doi:10.1145/3492805.3492813 Moritz Kreutzer, Georg Hager, Gerhard Wellein, Holger Fehske, and Alan R. Bishop. 2014. A Unified Sparse Matrix Data Format for Efficient General Sparse Matrix-Vector Multiply on Modern Processors with Wide SIMD Units. SIAM J. Sci. Comput. 36, 5 (Jan. 2014), C401–C423. arXiv:1307.6209 [cs] doi:10.1137/130930352 Neil Lindquist. 2023. Reducing Communication in the Solution of Linear Systems. Ph. D. Dissertation. The University of Tennessee, Knoxville, TN, USA. https://trace.tennessee.edu/handle/20.500.14382/29814 Neil Lindquist, Piotr Luszczek, and Jack Dongarra. 2022. Accelerating Restarted GMRES With Mixed Precision Arithmetic. IEEE Trans. Parallel Distrib. Syst. 33, 4 (April 2022), 1027–1037. doi:10.1109/TPDS.2021.3090757 Weifeng Liu and Brian Vinter. 2015. CSR5: An Efficient Storage Format for Cross-Platform Sparse Matrix-Vector Multiplication. In Proc. 29th ACM Int. Conf. Supercomput. (ICS ’15). Association for Computing Machinery, New York, NY, USA, 339–350. doi:10.1145/2751205.2751209 17
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Yuechen Lu and Weifeng Liu. 2023. DASP: Specific Dense Matrix Multiply-Accumulate Units Accelerated General Sparse Matrix-Vector Multiplication. In Proc. Int. Conf. High Perform. Comput. Netw. Storage Anal. ACM, Denver CO USA, 1–14. doi:10.1145/3581784.3607051 Marco Maggioni and Tanya Berger-Wolf. 2014. CoAdELL: Adaptivity and Compression for Improving Sparse MatrixVector Multiplication on GPUs. In 2014 IEEE Int. Parallel Distrib. Process. Symp. Workshop. IEEE, 933–940. doi:10.1109/IPDPSW.2014.106 Alexander Monakov, Anton Lokhmotov, and Arutyun Avetisyan. 2010. Automatically Tuning Sparse Matrix-Vector Multiplication for GPU Architectures. In High Performance Embedded Architectures and Compilers, David Hutchison, Takeo Kanade, Josef Kittler, Jon M. Kleinberg, Friedemann Mattern, John C. Mitchell, Moni Naor, Oscar Nierstrasz, C. Pandu Rangan, Bernhard Steffen, Madhu Sudan, Demetri Terzopoulos, Doug Tygar, Moshe Y. Vardi, Gerhard Weikum, Yale N. Patt, Pierfrancesco Foglia, Evelyn Duesterwald, Paolo Faraboschi, and Xavier Martorell (Eds.). Vol. 5952. Springer Berlin Heidelberg, Berlin, Heidelberg, 111–125. doi:10.1007/978-3-642-11515-8_10 Daichi Mukunoki, Masatoshi Kawai, and Toshiyuki Imamura. 2023. Sparse Matrix-Vector Multiplication with ReducedPrecision Memory Accessor. In 2023 IEEE 16th Int. Symp. Embed. MulticoreMany-Core Syst.–Chip MCSoC. IEEE, 608–615. doi:10.1109/MCSoC60832.2023.00094 Shun Murakami, Kazunori Yoneda, Takashi Iwamura, Masahiro Watanabe, and Yasushi Inoguchi. 2026. CoD-SELL: A Non-Zero Location Dictionary Compression Sparse Matrix Format for SpMV on GPU. IEEE Access 14 (2026), 17058–17068. doi:10.1109/ACCESS.2026.3659140 Yusuke Nagasaka, Akira Nukada, and Satoshi Matsuoka. 2016. Adaptive Multi-level Blocking Optimization for Sparse Matrix Vector Multiplication on GPU. Procedia Computer Science 80 (2016), 131–142. doi:10.1016/j.procs. 2016.05.304 Kengo Nakajima, Takseshi Ogita, and Masatoshi Kawai. 2021. Efficient Parallel Multigrid Methods on Manycore Clusters with Double/Single Precision Computing. In 2021 IEEE Int. Parallel Distrib. Process. Symp. Workshop IPDPSW. IEEE, 760–769. doi:10.1109/IPDPSW52791.2021.00114 Yuyao Niu, Zhengyang Lu, Meichen Dong, Zhou Jin, Weifeng Liu, and Guangming Tan. 2021. TileSpMV: A Tiled Algorithm for Sparse Matrix-Vector Multiplication on GPUs. In 2021 IEEE Int. Parallel Distrib. Process. Symp. IPDPS. IEEE, 68–78. doi:10.1109/IPDPS49936.2021.00016 Yvan Notay. 2000. Flexible Conjugate Gradients. SIAM J. Sci. Comput. 22, 4 (Jan. 2000), 1444–1460. doi:10.1137/ S1064827599362314 Antony Spyropoulos and Christos Antonopoulos. 2025. Numerical Study of Mixed Precision GMRES(m) Preconditioned by Deflation. Numer Algor 101 (May 2025), 2631–2657. doi:10.1007/s11075-025-02102-z Markus Steinberger, Rhaleb Zayer, and Hans-Peter Seidel. 2017. Globally Homogeneous, Locally Adaptive Sparse Matrix-Vector Multiplication on the GPU. In Proc. Int. Conf. Supercomput. (ICS ’17). Association for Computing Machinery, New York, NY, USA, 1–11. doi:10.1145/3079079.3079086 Kengo Suzuki. 2025. suzuki-hpc/F3R: v1.0.2. doi:10.5281/zenodo.16882405 Kengo Suzuki, Takeshi Fukaya, and Takeshi Iwashita. 2022. A New AINV Preconditioner for the CG Method in Hybrid CPU-GPU Computing Environment. Journal of Information Processing 30, 0 (2022), 755–765. doi:10. 2197/ipsjjip.30.755 Kengo Suzuki, Takeshi Fukaya, and Takeshi Iwashita. 2025. An Integer Arithmetic-Based AMG Preconditioned FGMRES Solver. ACM Trans. Math. Softw. 51, 1 (March 2025), 1:1–1:25. doi:10.1145/3704726 Kengo Suzuki and Takeshi Iwashita. 2025. A Nested Krylov Method Using Half-Precision Arithmetic. In Proc. Int. Conf. High Perform. Comput. Netw. Storage Anal. (SC ’25). Association for Computing Machinery, New York, NY, USA, 711–727. doi:10.1145/3712285.3759807 Seth Wolfgang, Skyler Ruiter, Marc Tunnell, Timothy Triche, Erin Carrier, and Zachary DeBruine. 2024. ValueCompressed Sparse Column (VCSC): Sparse Matrix Storage for Single-cell Omics Data. In 2024 IEEE Int. Conf. Big Data BigData. IEEE, Washington, DC, USA, 4952–4958. doi:10.1109/BigData62323.2024.10825091 Ichitaro Yamazaki, Erin Carson, and Brian Kelley. 2022a. Mixed Precision S-Step Conjugate Gradient with Residual Replacement on GPUs. In 2022 IEEE Int. Parallel Distrib. Process. Symp. IPDPS. IEEE, 886–896. doi:10.1109/ IPDPS53621.2022.00091 Ichitaro Yamazaki, Christian Glusa, Jennifer Loe, Piotr Luszczek, Sivasankaran Rajamanickam, and Jack Dongarra. 2022b. High-Performance GMRES Multi-Precision Benchmark: Design, Performance, and Challenges. In 2022 IEEEACM Int. Workshop Perform. Model. Benchmarking Simul. High Perform. Comput. Syst. PMBS. IEEE, 112–122. doi:10.1109/PMBS56514.2022.00015 18
PackSELL: A Sparse Matrix Format for Precision-Agnostic High-Performance SpMV
A P REPRINT
Lingqi Zhang, Jiajun Huang, Sheng Di, Satoshi Matsuoka, and Mohamed Wahib. 2025. Can Tensor Cores Benefit Memory-Bound Kernels? (NO!). In Proc. 17th Workshop Gen. Purp. Process. Using GPU (GPGPU ’25). Association for Computing Machinery, New York, NY, USA, 28–34. doi:10.1145/3725798.3725803 Yingqi Zhao, Takeshi Fukaya, and Takeshi Iwashita. 2023. Numerical Behavior of Mixed Precision Iterative Refinement Using the BiCGSTAB Method. Journal of Information Processing 31 (2023), 860–874. doi:10.2197/ipsjjip.31. 860 Yingqi Zhao, Takeshi Fukaya, Linjie Zhang, and Takeshi Iwashita. 2022. Numerical Investigation into the Mixed Precision GMRES(m) Method Using FP64 and FP32. J. Inf. Process. 30 (2022), 525–537. doi:10.2197/ipsjjip. 30.525
19