$g$MAGNUS: Fast SpGEMM on GPUs for Irregular Matrices via Hierarchical Multisplit

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

$g$MAGNUS: Fast SpGEMM on GPUs for Irregular Matrices via Hierarchical Multisplit

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributedcomputing parallelcomputing

distributed computing, parallel computing, cloud

arXiv:2607.22866v1 [cs.DC] 24 Jul 2026

g MAGNUS: Fast SpGEMM on GPUs for Irregular Matrices via Hierarchical Multisplit Jordi Wolfson-Pou

Ahmed Helal

Fabrizio Petrini

Intel Santa Clara, USA [email protected]

challenging: input row sizes vary widely, and some rows produce large intermediate products distributed across the column space. Scaling up these highly irregular matrices amplifies this issue, creating many heavy rows where intermediate products exceed the local-memory hash map capacity, forcing global memory fallbacks. One approach uses global spillover queues when local memory capacity is reached, requiring an unbounded number of passes of the queue, creating load imbalance when some work groups perform many more passes than others [18]. More recent algorithms use backup global memory accumulators accessed only when local memory overflows, creating uncoalesced accesses to global memory [19], [20]. Another approach is to sort and merge the intermediate products, which can be far too expensive for large numbers of intermediate elements. We present gMAGNUS (Matrix Algebra for Gigantic NUmerical Systems on GPUs), inspired by prior CPU work [21], I. I NTRODUCTION designed for massive irregular matrices with heavy rows. In the The sparse general matrix-matrix multiplication (SpGEMM), CPU approach, the goal was to maximize L2 cache reuse during which multiplies two sparse matrices A and B to produce accumulation. Each thread processed a distinct set of heavy a sparse matrix C, is critical to the performance of a wide rows, performing an outer product and an intra-row reordering of range of applications [1], including genome assembly [2]–[4], their intermediate products independently of other threads. This machine learning [5]–[9], algebraic multigrid [10], [11], and reordering introduced sufficient locality for the accumulator to fit graph analytics [12]–[17]. Although GPUs are widely used to in the L2 cache. The goal of gMAGNUS is to generate enough accelerate SpGEMM, their performance is often limited by the locality to perform accumulation entirely in local memory. Achieving this on GPUs requires a markedly different design, difficulty of mapping sparse structures onto the massive compute as work groups and threads must coordinate closely to maintain resources of GPUs. The challenge becomes more pronounced as matrices grow larger (both in the dimensions and nonzeros load balance and memory coalescing. To this end, gMAGNUS per row) and more structurally irregular, introducing high levels first uses an outer product to construct the global intermediate matrix Ĉ, then applies a novel hierarchical multisplit [22] to of imbalance and uncoalesced accesses to global memory. Gustavson’s method is the basis of most SpGEMM algorithms, reorder each row of Ĉ into chunks that form the rows of Ĉreord . processing matrices A and B in row-major order. Nonzero This hierarchical design supports an unbounded number of elements in each row of A are multiplied with corresponding chunks, so it scales to matrices of any size by further subdividing rows of B to form intermediate products, which are then summed rows through a small number of passes over Ĉreord . As a result, in an accumulator to produce the final row of C. The most the column range of each row in Ĉreord is much smaller than popular GPU accumulator is a hash map with column indices that of C, allowing a dense accumulator to fit entirely in local as keys. Upon insertion, new and existing values with duplicate memory and avoiding fallback to global-memory accumulation. keys are summed. Once all nonzeros are processed, the hash map Finally, gMAGNUS is input- and system-aware, automatically is sorted to produce the output row. The standard approach keeps selecting the number of chunks and hierarchy levels from the accumulators in local memory, where irregular accesses are much input matrix dimensions and available local memory. We implemented gMAGNUS in SYCL and CUDA, and present faster than global memory. This suffices for highly structured or small sparse matrices with sparse intermediate products. experimental results on both Intel Ponte Vecchio and NVIDIA Banded matrices exemplify this case due to their nearest- H200 GPUs. Five SpGEMM algorithms are used as baselines: neighbor structure, which produces small, localized intermediate Intel MKL [23], Kokkos [24], [25], NVIDIA cuSPARSE [26], products that can often be accumulated by a local-memory TileSpGEMM [27], and OpSparse [19]. We consider two datasets: dense accumulator. Random power-law matrices are far more the SuiteSparse matrix collection [28] and recursive power-law Abstract—We present gMAGNUS, a novel algorithm for sparse matrix-matrix multiplication (SpGEMM) of irregular matrices on GPUs. Such matrices often contain many heavy rows, those with large intermediate products that force local memory accumulators to spill to global memory. gMAGNUS addresses this by computing an intra-row reordering of intermediate products, subdividing heavy rows into independent chunks that can be accumulated completely in local memory. This reordering uses novel outer product and hierarchical multisplit operations. The algorithm is input- and system-aware, automatically determining the number of chunks and multisplit levels based on the input matrix dimensions and local memory size. Experimental results on two extensive datasets show that gMAGNUS achieves a geometric-mean speedup of 1.81–7.62× over five leading algorithms (including MKL and cuSPARSE) on Intel Ponte Vecchio and NVIDIA H200. Additionally, the core kernels of gMAGNUS are evaluated, achieving near-peak performance compared to their theoretical upper bound. Index Terms—SpGEMM, GPUs, Gustavson, Outer product

matrices (RMats) generated by PaRMAT [29]. Across all GPUs and datasets, gMAGNUS is faster than all baselines in terms of the geometric-mean speedup, with higher speedups for massive matrices that have large dimensions and intermediate products, such as large-scale RMat matrices. Additionally, we compare the core kernels of gMAGNUS (outer product and multisplit) to their “speed of light” performance, i.e., their theoretical peak performance. The results show that our core kernels achieve near-peak performance. The outer product often achieves 60-80% of the peak for large test matrices. The multisplit kernel, which has been studied for more regular applications, attains 40–50% of the peak, consistent with more extensive studies [22]. II. BACKGROUND A. Notation and Definitions Let A be a sparse matrix, with A.n, A.m, and A.nnz denoting its numbers of rows, columns, and nonzeros, respectively. We use direct indexing: A[i,:], A[:,j], and A[i,j] to denote row i, column j, and element (i,j), respectively. Sparse matrices are stored in the widely used compressed sparse row (CSR) format, where A.col , A.val , and A.rowPtr store the column indices, values, and row pointers, and have sizes A.nnz , A.nnz , and A.n+1, respectively. When indexing into a CSR matrix, row pointers determine the accessed range, e.g., A[i, :] denotes elements [A.rowPtr [i] : A.rowPtr [i+1]) of A.col and A.val . Vectors are denoted in lowercase (e.g., x, v), with x.n denoting the size. Subarrays and submatrices are indexed as [start : end ), so x[start : end ) denotes the consecutive elements {x[start],...,x[end −1]}. A masked subarray x[v] denotes the elements {x[v[0]],...,x[v[v.n− 1]]}. Similarly, A[v,:] and A[:,v] denote the submatrices formed by selecting the rows or columns indexed by v. On GPUs, we follow SYCL terminology: work group (thread block), subgroup, and local memory (shared memory). A device queue denotes a SYCL queue (CUDA stream). We use parfor to denote a GPU-style strided parallel loop over i ∈ [start,end ). B. Gustavson’s Method Gustavson’s method [30] is the basis for most GPU SpGEMM algorithms, and can be written as: A.m−1 X C[i, :] = A[i,j]×B[j, :]. (1) j=0

Algorithm 1: Gustavson’s Method Input: A, B, C.rowPtr Output: C.col , C.val 1 for i ∈ [0 : A.n) do 2 for j ∈ [A.rowPtr [i] : A.rowPtr [i+1]) do 3 for k ∈ [B.rowPtr [A.col[j]] : B.rowPtr [A.col[j]+1]) do 4 accumBuff [B.col[k]] += A.val[j]×B.val[k] 5 end 6 end 7 C[i,:] ← SortAndWrite(accumBuff ) 8 end

Algorithm 2: Outer Product Input: A, B, C.rowPtr Output: C.col, C.val 1 parfor i ∈ [0 : A.m) do 2 for j ∈ [ACSC .colP tr[i] : ACSC .colP tr[i+1]) do 3 for k ∈ [B.rowPtr [i], B.rowP tr[i+1]) do 4 dest ← AtomicAdd(Ĉ.rowPtr [A.row [j]+1], 1) 5 Ĉ.col[dest] ← B.col[k] 6 Ĉ.val[dest] ← A.val[j]×B.val[k] 7 end 8 end 9 end 10 parfor i ∈ [0, C.n) do 11 C[i,:] ← Accum(Ĉ[i,:]) 12 end

critical to efficiently handle the irregular accesses that arise during accumulation. The most common accumulator designs are dense arrays, hash maps, and sort-based algorithms. For “well-behaved” matrices, such as banded matrices, structural locality yields small, localized intermediate products, allowing a shortened dense accumulator whose size is proportional to the bandwidth rather than C.m. For more general matrices, hash maps and sort-based accumulators are more robust, since a dense accumulator requires storage proportional to C.m, which is often prohibitively expensive. However, for highly irregular matrices, such as random power-law matrices, many rows produce intermediate products that exceed local memory capacity, forcing all three approaches to use global memory and requiring more complex algorithms. We refer to such rows as heavy rows, and to the remaining rows as light rows. C. ESC and The Outer Product In some variants of Gustavson’s method, such as expand–sort– contract (ESC) [31], intermediate products are explicitly stored. A first pass expands these products to form the intermediate matrix Ĉ, and a subsequent accumulation pass merges rows of Ĉ to produce C. Alternative approaches adopt an outer-product formulation [21], [32]–[36], defined as: A.m−1 X C= A[:,i]⊗B[i,:], (2)

That is, the nonzeros in row i of A scale the corresponding rows of B, and the resulting scaled rows are summed to produce row i of C. Because the rows of C are independent, Gustavson’s method naturally parallelizes over rows. Algorithm 1 shows a CSR version of the numeric phase of Gustavson’s method. The numeric phase populates C.col and C.val , while the symbolic i=0 phase, performed beforehand, computes C.rowPtr and C.nnz where each rank-1 update (or slice) is given by the outer product (see [1] for more on size prediction). To combine the intermediate of column i of A and row i of B. This formulation mitigates product of a row, that is, the set of scaled rows A[i,j]×B[j, :], a key drawback of Gustavson’s method, namely uncoalesced an accumulator accumBuff sums the contributions on the fly. accesses and limited reuse in B, at the cost of explicitly After accumulation, a pass over accumBuff sorts and writes constructing Ĉ. Algorithm 2 shows a simple parallel version the row to C; in some implementations, the entire matrix C of the outer product. The first pass performs the construction is sorted afterwards. In typical parallel GPU implementations, of Ĉ using a parallel loop over the slices. The two inner loops the outer loop is partitioned across work groups and the inner illustrate the reuse of B, where row i of B is reused for each loops across threads within a work group. nonzero in column i of A. Since each slice modifies rows of Ĉ A key consideration in this on-the-fly accumulation is whether corresponding to the row indices in column i of A, an atomic the accumulator accumBuff fits in local memory, which is fetch-and-add allows for multiple slices to write to the same row.

While outer product improves reuse and mitigates the uncoalesced reads in Gustavson’s method, it introduces new challenges on GPUs. First, load balancing remains critical and requires careful distribution of work across all loop levels. Second, atomics can be avoided by expanding Ĉ in coordinate format, but the resulting triplets must be sorted during the final accumulation pass, which is prohibitive for large intermediate products. Third, achieving coalesced writes to Ĉ further requires careful partitioning of loops so that consecutive threads write to consecutive addresses. Finally, outer product primarily addresses inefficiencies in reading the input matrices, but the challenge of accumulating heavy rows still remains.

𝐴

𝐵

0 1 2 3 4 5 6 7

0 1 2

0 1 1

1 0 2 0

4 0

5 7

Outer product 1

Chunk boundary

መ 𝑐𝑜𝑙 0 5 0 7 1 3 5 0 5 0 𝐶. መ 𝑣𝑎𝑙 𝐴11𝐵10 𝐴11𝐵15 𝐴12𝐵20 𝐴12𝐵27 𝐴15𝐵51 𝐴15𝐵53 𝐴15𝐵55 𝐴21𝐵10 𝐴21𝐵15 𝐴24𝐵40 𝐶. Reorder

D. Related Work We focus on recent GPU SpGEMM algorithms most relevant to gMAGNUS. A comprehensive survey is available in [1]. Most prior work has focused on optimizing local-memory accesses and improving accumulator load balancing within Gustavson’s method [18]–[20], [37], [38]. A common strategy is to use one or more accumulator types and bin rows by size. Smaller work groups process smaller rows, while rows that exceed the local-memory capacity of the largest work group either require multiple passes [18], [20] or fall back to a global-memory accumulator [19], [24], [37]. Hybrid accumulators are also common, for example using hash maps for sparser rows or tiles and dense accumulators otherwise [20], [27], or using different accumulator types in the symbolic and numeric phases [38]. The ESC-based methods [31], [39]–[43] are generally less competitive than Gustavson-based approaches because they explicitly generate and sort the intermediate product in global memory, which is expensive. Recent ESC variants reduce the cost of global sorting by sorting smaller subarrays in local memory [41], but this requires additional passes over global memory to generate and merge those subarrays. Other work explores storage formats beyond CSR [27], [44], of which TileSpGEMM [27] is a leading approach. Inspired by tiled dense GEMM, TileSpGEMM converts the input matrices into 16×16 tiles, storing only nonzero tiles, which improves locality in both the input matrices and the accumulator. However, this introduces a substantial CSR-to-tile setup cost, and the output matrix C must later be converted from tiled form back to CSR.

0 1 3 1 3 1 0 0 1 𝐶መreord . 𝑐𝑜𝑙 0 𝐶መreord . 𝑣𝑎𝑙 𝐴11𝐵10 𝐴12𝐵20 𝐴15𝐵51 𝐴15𝐵53 𝐴11𝐵15 𝐴12𝐵27 𝐴15𝐵55 𝐴21𝐵10 𝐴24𝐵40 𝐴21𝐵15 Chunk 1

Chunk 0

Chunk 2

Chunk 3

Accumulate 𝐶. 𝑐𝑜𝑙

𝐴11 𝐵10 𝐴11 𝐵15 𝐴21 𝐵10 𝐶. 𝑣𝑎𝑙 + 𝐴15𝐵51 𝐴15𝐵53 + 𝐴12𝐵27 + 𝐴21𝐵15 𝐴12 𝐵20 𝐴15 𝐵55 𝐴24 𝐵40

Fig. 1. Example of processing of heavy rows by gMAGNUS, where two chunks and one level are used to compute rows 1 and 2 of C. The dashed lines represent the row pointers, whereas the column and value arrays are explicitly written.

for dense accumulators to fit in local memory (subsection III-D shows how the number of chunks is determined at runtime). Because it explicitly stores intermediate products, gMAGNUS is similar in spirit to ESC-based algorithms [31]. However, by expanding only the intermediate products of heavy rows, it requires substantially less memory than prior ESC approaches. Our optimized multisplit kernel performs a range-based multisplit operation, permuting a set of key-value pairs into chunks (also referred to as buckets or bins), where elements belonging to the same chunk are stored contiguously. The set of key-value pairs are rows of Ĉ and the chunks are rows of Ĉreord , e.g., if Ĉ stores two heavy rows that must be split into two chunks each, Ĉreord will have 4 rows, as shown in Figure 1. While Ĉreord has twice as many rows as Ĉ, it has half as many columns. The reduced column range enables local memory-only dense accumulation, while the increased number of rows exposes additional parallelism by effectively splitting each row of Ĉ across multiple work groups. Within a chunk, elements are unordered and processed in later III. MAGNUS stages by the accumulator. To construct the rows of Ĉreord , column A. Overview indices are shifted into their local chunk range. in Figure 1, gMAGNUS addresses heavy rows in large irregular matrices by chunks 1 and 3 shift indices from [4,8) to [0,4). During the final performing an intra-row reordering of the intermediate matrix write to C, these indices are shifted back to their original range. Ĉ (henceforth, Ĉ denotes the intermediate matrix generated Algorithm 3 shows the end-to-end gMAGNUS algorithm, only by the heavy rows). It splits rows into independent chunks where all inputs and outputs reside in device memory. The to produce a second intermediate matrix, Ĉreord , whose structure notation ⟨y1 , y2 , ... ⟩ ← Function⟨Q⟩(x1 , x2 , ... ) denotes a both enables local-memory accumulation and exposes additional function that launches one or more kernels on device queue Q. intra-row parallelism. This reordering is implemented through an Light and heavy rows are processed concurrently on separate outer product and a novel hierarchical multisplit operation, both asynchronous device queues. In the setup phase, a preprocessing optimized for load balancing and coalesced memory accesses. step queries the device and computes the heavy-row threshold τ , The final accumulation step uses a local-memory hybrid strategy: the number of chunks per row per level vchunks , and the number of hash map accumulation for light rows of C and Ĉreord , and levels vchunks .n (as detailed later in subsection III-D). τ is derived dense accumulation for heavy rows of Ĉreord . Our approach from the local memory size in bytes, sLM , which determines the is input- and system-aware, choosing the number of chunks maximum accumulator capacity. This O(1) computation overlaps per row based on C.m and the local-memory capacity. These with InterOffsets(), which computes per-row intermediate parameters determine the minimum number of chunks needed

Algorithm 3: gMAGNUS Input: A, B Output: C 1 /* Setup phase */ 2 ⟨τ, vchunks ⟩ ← GetMagnusParams(C.m, sLM ) 3 Ĉ.rowPtr ← InterOffsets⟨Qmain ⟩(A, B) 4 ⟨xheavy , xlight ⟩ ← ThreshPartition⟨Qmain ⟩(Ĉ.rowP tr, τ ) 5 Synchronize() 6 /* Outer product */ 7 ⟨ÃCSC ,B̃⟩ ← MaskedCsr2Csc⟨Qheavy ⟩(A[xheavy ,:],B) 8 Ĉ ← OuterProduct⟨Qheavy ⟩(ÃCSC , B̃) 9 /* Reorder */ 10 for i ∈ [0 : vchunks .n) do 11 Ĉreord .rowP tr ← Histogram⟨Qheavy ⟩(Ĉ, vchunks [i]) 12 if i < vchunks .n−1 then 13 ⟨x̂heavy , x̂light ⟩ ← ThreshPartition⟨Qheavy ⟩(Ĉreord .rowP tr, τ ) 14 end 15 Ĉreord ← Multisplit⟨Qheavy ⟩(Ĉ,Ĉreord , vchunks [i]) 16 if x̂heavy .n == 0 then 17 break 18 end 19 Ĉ ← Ĉreord [x̂heavy ,:] 20 end 21 C.rowPtr [xheavy ] ← Symbolic⟨Qheavy ⟩(Ĉreord ) 22 C.rowPtr [xlight ] ← Symbolic⟨Qlight ⟩(A[xlight ,:],B) 23 Synchronize() 24 C.rowPtr ← PrefixSum⟨Qmain ⟩(C.rowPtr ) 25 Synchronize() 26 C[xheavy ,:] ← Numeric⟨Qheavy ⟩(Ĉreord ) 27 C[xlight ,:] ← Numeric⟨Qlight ⟩(A[xlight ,:],B) 28 Synchronize()

product sizes and generates Ĉ.rowPtr using a device-wide prefix sum. The intermediate product sizes together with τ are used to partition rows into light and heavy groups using vendor primitives, such as cub::DevicePartition::Flagged() in CUDA. The remaining steps execute on the two concurrent queues for light and heavy rows. To process heavy rows, we first perform the outer product (lines 7-8), which requires matrix A in CSC format. Hence, we carry out a masked CSR-to-CSC conversion using a key-value sorting approach [45]. Additionally, we form the masked matrix B̃, which contains the rows of B corresponding to the nonzero columns of ÃCSC . The outer product then populates the intermediate matrix Ĉ. Our hierarchical multisplit (shown at lines 10–20) consists of multiple operations. Histogram() counts the number of elements per chunk and performs a device-wide prefix sum to compute Ĉreord .rowPtr . Multisplit() then permutes the entries of Ĉ to populate Ĉreord . Then we again use τ to split the remaining chunks into light and heavy groups, since only heavy chunks need further reordering, reducing data volume for subsequent multisplit operations. If there are only light chunks, we break from the loop (line 16). At the last level (i.e., when i == vchunks .n − 1), Ĉreord .m is small enough for the dense accumulator to fit in local memory, meeting the gMAGNUS design goal. The remaining steps of gMAGNUS follow the standard SpGEMM workflow, with accumulation in the symbolic and numeric stages (see subsection III-E). For algorithmic readability, we show explicit gather operations

Flattened space partitions

Work group 0

Work group 1

𝑡0 𝑡1 𝑡2 𝑡3 𝑡0 𝑡1 𝑡2 𝑡3 𝑡0 𝑡1 𝑡2 𝑡0 𝑡1 𝑡2 𝑡3 𝑡0 𝑡1 𝑡2 𝑡3 𝑡0 𝑡1 slice 0

slice 1

slice 2

slice 3

𝐴ሚ CSC . 𝑟𝑜𝑤

෨ 𝑐𝑜𝑙 𝐵.

𝑡0 𝑡1 𝑡2

𝑡2 𝑡3 𝑡0

𝑡0

𝑡0 𝑡1

0 1 3 𝑡3 𝑡0 𝑡1

1 2 𝑡1 𝑡2

𝑡1 𝑡2

3 𝑡3

0 𝑡1

1 3 𝑡3 𝑡0

መ 𝑐𝑜𝑙 𝐶.

Fig. 2. Example partitioning of ÃCSC and B̃ for the outer product kernel. The top sequence of elements shows the partitioning of the flattened index space, which maps to the partitioning of ÃCSC and B̃ shown below it. The blue dashed boundary lines represent row pointers. The elements read by each thread are also shown, as calculated by GetEntry() and GetExcess() in Algorithm 4.

like Ĉ ← Ĉreord [x̂heavy ,:]. However, in the implementation, this operation is fused with our GPU kernels. For example, this gather is actually performed within Multisplit() during the final write to global memory, ensuring heavy chunks are ordered first so that subsequent multisplits operate on contiguous elements for optimal load balancing. Additionally, the loop over levels uses double buffering: one buffer holds the input Ĉ, produced either by the outer product or by the reordering step from the previous iteration, while the other holds the output Ĉreord . After the loop completes, the buffer containing the final input is deallocated. B. Outer Product Kernel The outer product kernel (OuterProduct() in Algorithm 3) balances memory traffic across work groups and maximizes coalesced accesses. We use a static partitioning computed locally without any inter-thread coordination. Each work group reads the same number of ÃCSC –B̃ element pairs, multiplies them, and writes the results to Ĉ. We achieve this with a static partitioning of a flattened index space, as shown in the top array of Figure 2, where two work groups are assigned approximately the same number of slice (input) and intermediate (output) elements. This flattening unrolls the slices and requires mapping indices in the flattened space back to elements in ÃCSC and B̃. Both work groups read consecutive elements of ÃCSC and B̃ and write consecutive elements of B̃ to Ĉ. The mapping from the flattened space to ÃCSC and B̃ must handle two important edge cases. First, elements of ÃCSC and B̃ read by different work groups may overlap, for example, see slice 1 (column 1 and row 1 of ÃCSC and B̃, respectively). Second, as a result of overlap, work groups may read partial rows of B̃, as shown in the example where work group 0 reads two of the three elements on its last pass over row 1 of B̃, and work group 1 reads the remaining last element on its first pass over the same row of B̃. Algorithm 4 shows our outer product kernel, which has three main steps: (1) local load balancing, (2) updating Ĉ.rowPtr , and (3) populating Ĉ. In the algorithm, bold symbols are global and symbols with an ℓ subscript are local memory arrays. The map-

ping xmap is derived from the light-heavy row partitioning and maps the heavy rows, which are nonconsecutive in the original ordering, to consecutive rows in Ĉ. The first launch parameter is the number of work groups, and the second is the work group size. The slice offsets are computed within OuterProduct() as the prefix sum of sliceOffsets[i + 1] = (ÃCSC .colPtr [i + 1] − ÃCSC .colPtr [i]) × (B.rowPtr [i + 1] − B.rowPtr [i]) (the number of elements per slice) for i ∈ {0,1,...,A.m−1}.

flattened space to an element in ÃCSC as:

i−sliceOffsets[slice] , B̃.rowPtr [slice +1]− B̃.rowPtr [slice] (3) where slice = sliceRange.start and slice = sliceRange.end for the start and end elements, respectively, where ... in the Algorithm 4: OuterProductKernel function call denotes that the remaining arguments consist of Input: ÃCSC , B̃, sliceOffsets, xmap , ngroupElems the global variables and partition information. The first term Output: Ĉ represents the starting column of ÃCSC , and the second term Launch: ⌈ nĈ.nnz ⌉, nworkGroup groupElems gives the starting element within the column. The numerator 1 /* Local load balancing */ 2 groupRange ← GroupRange(ngroupElems ) in the second term denotes the mapping of the flattened index 3 sliceRange ← Search(sliceOffsets, groupRange) i to the local index within the current slice. Dividing by the 4 aRange ← ARange(groupRange,...) B̃ row size yields the element in the current column of ÃCSC . 5 /* Update Ĉ.rowP tr */ 6 slice ← sliceRange.start In step 2 (lines 7–15), updates to Ĉ.rowPtr determine the 7 parfor i ∈ [aRange.start : aRange.end) do write locations for populating Ĉ. The objective is to write 8 while i ≥ ÃCSC .colP tr[slice+1] do rows (or partial rows) of B̃ to contiguous locations in Ĉ. Each 9 slice ← slice+1 work group iterates over elements of ÃCSC in its partition 10 end 11 rowNnzB ← B̃.rowP tr[slice+1]− B̃.rowP tr[slice] and performs atomic updates to Ĉ.rowPtr using row sizes 12 excess ← GetExcess(i, slice, ...) of B̃. Because AtomicAdd() returns the previous value 13 rowNnzB ← rowNnzB −(excess.first +excess.last) of the target address, we store these return values in local 14 offsetsLocal ℓ [i−aRange.start] ← memory as the eventual write locations for populating Ĉ. AtomicAdd(Ĉ.rowP tr[xmap [ÃCSC .row[i]]+ 1], rowNnzB ) GetExcess() detects partial B̃ rows and adjusts the atomic 15 end updates accordingly. It returns the number of excess elements 16 SyncThreads() in the current B̃ row, that is, the elements that belong to one 17 /* Populate Ĉ.col and Ĉ.val */ 18 slice ← sliceRange.start or more other work groups, using a calculation similar to 19 parfor i ∈ [groupRange.start : groupRange.end) do GetAEntry(). Note that we update Ĉ.rowP tr, which was 20 while i ≥ sliceOffsets[slice +1] do computed earlier in the setup phase. Our implementation is 21 slice ← slice+1 22 end in-place: the setup phase modifies elements [2 : Ĉ.n+2), and the 23 /* Compute read locations of ÃCSC and B̃ */ outer product modifies the shifted elements [1 : Ĉ.n+1), yielding 24 ⟨entryA,entryB ⟩ ← GetABEntries(i, slice, ...) the correct Ĉ.rowP tr at completion of the outer product. 25 /* Retrieve column and value data from ÃCSC and B̃ */ In step 3 (lines 19–35), we populate Ĉ by writing consecutive 26 colB ← B̃.col[B̃.rowP tr[slice]+entryB] elements of B̃ to consecutive elements in Ĉ. These accesses are 27 valB ← B̃.val[B̃.rowP tr[slice]+entryB] 28 valA ← ÃCSC .val[ÃCSC .colP tr[slice]+entryA] not fully coalesced because the rows of B̃ are not necessarily 29 /* Compute write location in Ĉ */ multiples of the cache line (Intel) or sector length (NVIDIA). In 30 excess ← GetExcess(i,slice,...) practice, heavy rows usually consist of dense B̃ rows, enabling 31 dest ← offsetsLocal ℓ [ÃCSC .colP tr[slice]− many coalesced writes and good performance for gMAGNUS. aRange.start+entryA]+entryB −excess.f irst The function GetABEntries() calls GetAEntry() as 32 /* Write to Ĉ */ 33 Ĉ.col[dest] ← colB well as GetBEntry(), which is identical to GetAEntry() 34 Ĉ.val[dest] ← valA×valB except it uses a modulo instead of a division in the second 35 end term of Equation 3. When writing rows of B̃, we again adjust Local load balancing (lines 2–4) maps work group partitions for partial rows. Note that in both steps 2 and 3, a while loop to their starting elements in the slice offsets, ÃCSC , and B̃, adjusts the current slice, causing subgroup/warp divergence. where work group i processes elements [i ∗ ngroupElems , (i + However, because the outer product is applied to heavy rows, 1)∗ngroupElems ) in the flattened space (subsection III-D shows each work group usually touches only a few slices, keeping how ngroupElems is determined). GroupRange() returns the the number of while-loop iterations low. In summary, our static partitioning yields coalesced reads with flattened index range, and Search() (a modified binary search) returns the largest slice offset less than i ∗ ngroupElems . high reuse for input matrices. Writes to Ĉ are often highly coaEach thread performs its own O(log(A.m)) binary search1 . lesced as well since heavy rows tend to involve denser rows of B. ARange() returns the range of elements in ÃCSC for a work Unlike previous work that use dynamic global load balancing [34], group, which calls GetAEntry(groupRange.start, ... ) our local-only partitioning is lightweight, determined only by (not shown in the algorithm). GetAEntry(i,...) maps i in register-level calculations. Additionally, since our algorithm is only applied to heavy rows, which means that ÃCSC ≪ Ĉ.nnz, both the CSR-to-CSC conversion cost and the cost of the atomic operations in the first step of the outer product are often small 1While a group- or subgroup-wide parallel search could reduce this cost, our approach performed well in practice. compared to populating Ĉ, as shown in Section IV. ÃCSC .colPtr [slice]+

C. The Multisplit Kernel

Algorithm 5: MultisplitKernel

Input: Ĉ, xmap ,nchunks , ngroupElems The construction of Ĉreord involves three kernels: a histogram Output: Ĉ l reord m kernel followed by a device-wide prefix sum computes Launch: nĈ.nnz , nworkGroup Ĉreord .rowP tr, and a multisplit kernel [22] populates Ĉreord . workGroup Local load balancing */ We describe only the multisplit kernel, as the histogram kernel 12 /* groupRange ← GroupRange(ngroupElems ) is a subset of its steps, namely local histogramming and updates 3 rowRange ← Search(Ĉ.rowP tr, groupRange) to Ĉreord .rowP tr. Our design goal for multisplit is to divert the 4 /* Local histogram */ highly irregular permute phase to local memory, maximizing 5 row ← 0, offsetsLocal ℓ ← 0 6 parfor i ∈ [groupRange.start : groupRange.end) do coalesced accesses to global memory. Algorithm 5 shows 7 while i ≥ Ĉ.rowP tr[rowRange.start +row+1] do our multisplit kernel. As in the outer product, xmap and a 8 row ← row+1 end static partitioning of the flattened index space of Ĉ are used. 9 chunk ← GetChunk(Ĉ.col[i],row,...) Compared to the outer product, the partitioning is simplified 10 11 AtomicAdd(offsetsLocal ℓ [chunk+2], 1) since we only need to index into Ĉ instead of both ÃCSC and 12 end B̃, removing the need for GetEntry() and GetExcess(). 13 SyncThreads() Within a work group, multisplit has the following steps: 14 /* Updates to Ĉreord .chunkP tr */ 15 parfor i ∈ [0 : ngroupRows ×nchunks ) do local histogram, global update to Ĉreord .rowP tr, local prefix 16 chunkElems ← offsetsLocal ℓ [chunk+2] sum, local permute, and global write to Ĉreord . There are 17 chunk ← xmap [rowRange.start ×nchunks +i]+1 offsetsGlobal ℓ [i] ← four local memory arrays: offsetsLocal ℓ , offsetsGlobal ℓ , 18 colBuff ℓ , and valBuff ℓ . colBuff ℓ and valBuff ℓ store the 19 end AtomicAdd(Ĉreord .chunkPtr [chunk ], chunkElems) local chunks, i.e., reordered column indices and values 20 SyncThreads() of Ĉreord . offsetsLocal ℓ stores the local chunk offsets, 21 /* Local offsets via prefix sum */ PrefixSumInPlace(offsetsLocal ℓ [2 : ngroupRows ×nchunks +2)) where the elements of local chunk i are stored at positions 22 23 SyncThreads() [offsetsLocal ℓ [i] : offsetsLocal ℓ [i]+1) of colBuff and valBuff . 24 /* Permute in local memory */ offsetsGlobal ℓ stores the global write positions: elements at 25 row ← 0 positions [offsetsLocal ℓ [i] : offsetsLocal ℓ [i] + 1) of colBuff 26 parfor i ∈ [groupRange.start : groupRange.end) do 27 while i ≥ Ĉ.rowP tr[rowRange.start +row+1] do and valBuff are written to Ĉreord .col and Ĉreord .val starting at 28 row ← row+1 29 end position offsetsGlobal ℓ [i], respectively. chunk ← GetChunk(Ĉ.col[i],row,...) The size of colBuff ℓ and valBuff ℓ is ngroupElems . The size 30 dest ←AtomicAdd(offsetsLocal ℓ [chunk+1], 1) of offsetsGlobal ℓ and offsetsLocal ℓ is ngroupRows ×nchunks and 31 32 colBuff ℓ [dest] ← Ĉ.col[i] ngroupRows ×n l m chunks +1, respectively. The quantity ngroupRows = 33 valBuff ℓ [dest] ← Ĉ.val[i] ngroupElems −1 +1 denotes the maximum number of rows covered 34 end τ 35 SyncThreads() by any work group. The ngroupRows ×nchunks term accounts for 36 /* Coalesced writes to global memory */ row-boundary crossings, where we need to store a distinct set 37 row ← 0 of nchunks offsets per row. The additional 2 enables multiple 38 parfor i ∈ [groupRange.start : groupRange.end) do while i ≥ Ĉ.rowP tr[rowRange.start +row+1] do in-place histogramming passes over offsetsLocal ℓ , first to count 39 40 row ← row+1 local chunks and then to update offsets during the local permute. 41 end chunk ← In the local histogram step (lines 5-12), the number of elements 42 GetChunk(colBuff ℓ [i−groupRange.start],row,...) in each local chunk is counted. Column indices of Ĉ are 43 dest ← offsetsGlobal ℓ [chunk]+i−offsetsLocal ℓ [chunk] mapped to chunks which returns chunk ← 44 Ĉreord .col[dest] ← colBuff ℓ [i]−chunk× Ĉreord .m j using GetChunk(), k 45 Ĉreord .val[dest] ← valBuff ℓ [i] row×nchunks + Ĉ.col[i] , where row is the local row index in Ĉreord .m 46 end the current partition. The local histogram is then used to update Ĉreord .rowP tr (lines 15–19). As in the outer product kernel, we our design goal of coalesced reads and writes, with irregular employ an in-place update scheme for Ĉreord .rowP tr, where elenonconsecutive writes diverted to local memory during the local ments [2,Ĉreord .n+2) were updated in the histogram kernel prior permute phase. to the multisplit kernel. Next, a work group-wide parallel prefix Note that, as in the outer product, we do not achieve perfectly sum (line 22) of offsetsLocal ℓ computes the local chunk offsets coalesced writes since the local chunk counts are not guaranteed using vendor primitives, e.g., cub::DeviceScan on NVIDIA to align with the cache line (or sector) size. This behavior is GPUs. These offsets are used to then perform the local permute to expected given the irregularity of SpGEMM, yet the approach get the local chunks (lines 26-34). The reordered column indices still yields substantially more coalesced writes than performing and values are stored in colBuff ℓ and valBuff ℓ at locations deterthe reordering directly in global memory. Our experiments show mined by atomically updating offsetsLocal ℓ [i]. Finally, the local performance comparable to other multisplit implementations chunks are written to global memory, populating Ĉreord (lines 38- for more regular applications. Additionally, although not shown 46). To maximize coalesced writes, consecutive elements within a in Algorithm 5 for readability, we employ subgroup-private local chunk are written to consecutive elements in global memory: (warp-private) histograms, a common optimization that reduces elements [offsetsLocal ℓ [i],offsetsLocal ℓ [i+1]) of colBuff ℓ and atomic contention [46]. In this scheme, offsetsLocal has size ℓ valBuff ℓ are written to Ĉreord .col and Ĉreord .val, respectively, n groupRows × nchunks × nsubGroups + 1, where subgroup i updates starting at the locations stored in offsetsGlobal ℓ [i]. This achieves

elements [ngroupRows ×nchunks ×i,ngroupRows ×nchunks ×(i+1)). Our algorithm has several additional benefits. First, it performs few atomic operations relative to intermediate-element accesses: O(Ĉreord .n) atomics versus O(Ĉreord .nnz) reads and writes, where typically Ĉreord .nnz ≫ Ĉreord .n for heavy rows. Atomic contention is also low, since threads within a work group update distinct elements of Ĉreord .rowP tr, so at most ngroupRows threads access the same element (ngroupRows = 2 in our experiments). Second, we reduce register pressure to maximize occupancy. In particular, we avoid nested loops, which are often unrolled by the compiler and increase register usage. For example, in the final loop, even though the local intermediate elements are already ordered by chunk, we remap each element to its chunk number instead of using nested loops over chunks and elements. D. Parameters Selection In the initial host-side preprocessing step (GetMagnusParams() in Algorithm 3), the gMAGNUS parameters are computed using inputs that are easily queried at runtime. These inputs are C.m (the number of columns of C), and the size in bytes of the data types used to store the CSR arrays. We define sy as the number of bytes of the data type of variable y. All quantities in this section are integers, so division is assumed to use floor rounding unless a ceiling is explicitly specified. The computation of τ , the large chunk threshold (i.e., the maximum capacity of the numeric-phase local memory hash map), is given by: sLM τ= 1 , (4) α ×sC.col ×sC.val where α is the load factor used to reduce collisions. Our implementation of hash accumulators uses α = 2. We also have the maximum dense accumulator size for numeric phase as ndenseNumeric = sLM /(sC.val + sbitMap ), which we use to calculate nchunksRequired = C.m/ndenseNumeric , the minimum number of chunks per row required for local memory-only dense accumulation. Rounding to the nearest power of two allows us to use bit shift operations instead of division when mapping column indices to chunks, as explained in subsection III-C. Before using nchunksRequired to compute vchunks (an array of size nlevels containing the number of chunks per row per level), we need to calculate nchunksMax , the maximum number of chunks that our multisplit kernel can process in a single pass (a single level). Because the number of required chunks grows with C.m, nchunksMax represents the threshold at which offsetsLocal ℓ and offsetsGlobal ℓ no longer fit in local memory, beyond which our multilevel algorithm is needed. We compute nchunksMax by solving sLM = shisto ×(ngroupRows ×nsubGroups ×nchunksMax +1) + sĈ.rowP tr ×ngroupRows ×nchunksMax + (5) (sC.col +sC.val )×nminGroupElems for nchunksMax , where the right-hand side denotes the total size in bytes of all local memory arrays. The three terms represent the sizes in bytes of offsetsLocal ℓ , offsetsGlobal ℓ , and colBuff ℓ and valBuff ℓ . Multiplication by nsubGroups , the number of subgroups per work group, reflects histogram duplication across subgroups. nminGroupElems denotes the minimum number of elements per work group in the flattened partition. We set nminGroupElems = sLM 2×(sC.col +sC.val ) , i.e., half of the available local memory is allocated to the locally permuted intermediate elements. In practice, this choice balances local memory usage between the

locally permuted elements and the histogram arrays, which we found yields high local memory utilization. Solving gives us sLM −nminGroupElems ×(sC.col +sC.val )−shisto . nchunksMax = ngroupRows × nsubGroups ×shisto +sĈ.rowP tr (6) With nchunksRequired and nchunksMax , we can calculate the number of levels. If nchunksRequired < nchunksMax , we have one level with nchunksRequired chunks per row. Otherwise, we compute levels the number of levels by solving nnchunksMax = nchunksRequired . Intuitively, this corresponds to the number of times we need to divide C.m by nchunksMax to obtain our target local memory-only dense accumulator size. This represents the hierarchical component of gMAGNUS, giving us log2 (nchunksRequired ) nlevels = . (7) log2 (nchunksMax ) Since this gives us the minimum number of levels, one level may have fewer chunks than nchunksMax . We set this as the last level: vchunks [i] = nchunksMax for i ∈ [0, nlevels − 1) and levels −1 vchunks [nlevels − 1] = nchunksRequired − nnchunksMax . Using more chunks in earlier levels acts as a filter: subdividing the row into smaller parts produces fewer heavy chunks, increasing the chance of later levels reordering fewer elements. An additional optimization is that we floor nchunksRequired and nchunksMax to the nearest power of two. This lets us map column indices to chunks using lower-latency bit shifts instead of integer division. The final key parameter is the number of elements per work group in the flattened space, ngroupElems , which is kerneldependent. In general, we choose this quantity as large as possible to maximize either local memory utilization in the multisplit kernel or input reuse in the outer product kernel. For the multisplit kernel, ngroupElems is the number of elements that fit in the available shared memory after accounting for the local histograms, rounded down to the nearest multiple of the work-group size. For the outer product kernel, we use ngroupElems = s sLM , Ĉ.rowP tr which is the maximum size of offsetsLocal ℓ . E. Accumulation gMAGNUS is agnostic to the specific local-memory accumulator: any optimized local-memory accumulator can be used for light rows and chunks, and any optimized local-memory dense accumulator can be used for heavy chunks. While we use simplified implementations of both, integrating more optimized library accumulators into gMAGNUS could further improve performance and is left for future work. Our implementation uses a variant of the common binning approach [19], in which lighter rows are assigned smaller work-group sizes according to their intermediate product sizes. This requires a setup phase that uses a kernel similar to multisplit to map rows and chunks to bins, where each bin corresponds to a work-group size. We then launch one asynchronous kernel per bin, with one work group per row or chunk within the bin. For light rows, we use a standard hash map with modulo hashing and linear probing, followed by a bitonic sort to produce the final row of C (or subrow, for light chunks). For heavy chunks in the numeric phase, we use dense accumulation to merge elements, followed by a work-group-wide prefix sum over the bitmap to generate the sorted final chunk. The index arrays returned by FlaggedPartition() provide the mapping back to the original row and chunk

orderings (the inverse maps of those in Algorithm 4 and Algorithm 5), which we use when writing the final rows of C. IV. E XPERIMENTAL R ESULTS A. Experimental Setup We evaluate our SYCL and CUDA implementations against five state-of-the-art SpGEMM implementations spanning vendor libraries and prior academic work: Intel MKL [23], Kokkos [24], [25], NVIDIA cuSPARSE [26], TileSpGEMM [27], and OpSparse [19]. All methods are compiled with CUDA 13.1.1 or oneAPI 2025.3.0. TileSpGEMM and OpSparse are selected as academic baselines because they are recent, algorithmically distinct, and have been evaluated against a broad set of prior approaches. For cuSPARSE, we test all three available algorithms, including five separate runs of algorithm 3, with parameters 0.1, 0.2, 0.3, 0.4, and 0.5. This parameter controls the fraction of intermediate products processed at once as a way to reduce memory consumption. We report the minimum execution time across all seven cuSPARSE configurations, where algorithm 1 is usually the fastest followed by 2 and then 3. For all algorithms, we report the total end-to-end time, which is the sum of any pre-processing and/or setup, compute (symbolic and numeric), and post-processing phases. For all experiments, we perform one warmup run followed by 20 timed runs and report the mean time. Some baselines failed to run to completion on certain matrices due to segmentation faults, illegal memory accesses, out-of-memory errors, or timeouts (using a cutoff of 20 minutes). Table I shows our test data center GPUs: Intel Ponte Vecchio (PVC) 1100 and NVIDIA H200. The SYCL implementations (MKL and Kokkos) are evaluated on the Intel PVC GPU, and the CUDA implementations are evaluated on H200 (cuSPARSE, Kokkos, TileSpGEMM, and OpSparse). TABLE I S PECIFICATIONS OF THE TEST GPU S . Architecture Memory Peak memory bandwidth Max local memory per Xe-core / SM L1 / shared cache per Xe-core / SM L2 cache (global)

Intel PVC 1100 48 GB ∼1.2 TB/s 128 KB 512 KB 408 MB

NVIDIA H200 141 GB ∼4.8 TB/s 228 KB 256 KB 50 MB

TABLE II P ROPERTIES OF THE REPRESENTATIVE S UITE S PARSE MATRICES . Matrix para-9 Stanford Berkeley soc-Slashdot0902 HTC 336 4438 bloweya a0nsdsil c-57 brainpc2 TSOPF FS b39 c7 hangGlider 5 in-2004 net150 eu-2005 vsp south31 slptsk vsp model1 crew1 cr42 south31 c-big pkustk12 mult dcop 02 wb-edu rajat28

A.n 155,924 683,446 82,168 226,340 30,004 80,016 37,833 27,607 28,216 16,011 1,382,908 43,520 862,664 39,668 45,101 345,241 94,653 25,187 9,845,725 87,190

A.nnz 5,416,358 7,583,376 948,464 904,522 150,009 355,034 405,197 179,395 730,080 155,246 16,917,053 3,121,200 19,235,140 379,828 379,952 2,341,011 7,512,317 193,276 57,156,537 607,235

A.nnz A.n

34.7 11.1 11.5 4.0 5.0 4.4 10.7 6.5 25.9 9.7 12.2 71.7 22.3 9.6 8.4 6.8 79.4 7.7 5.8 7.0

A2 .nnz A.n

A2 .nnz 78,088,896 78,130,972 81,362,487 83,300,312 100,360,010 175,955,042 177,766,591 190,743,619 199,299,980 202,569,269 213,255,458 238,012,852 284,177,131 389,540,274 394,768,783 447,991,461 474,804,911 518,559,249 630,077,764 898,546,696

500.8 114.3 990.2 368.0 3344.9 2199.0 4698.7 6909.2 7063.4 12651.9 154.2 5469.0 329.4 9820.0 8753.0 1297.6 5016.3 20588.4 64.0 10305.6

TABLE III P ROPERTIES OF THE REPRESENTATIVE RM AT MATRICES . F OR THIS SET, B.nnz/B.n IS FIXED AT 8. Matrix Scale 15 16 17 18 19 20

A.nnz/A.n 4 C.nnz/C.n 477.6 633.5 835.0 1112.7 1471.9 1944.1

8 Ĉ.nnz/Ĉ.n 819.9 1068.2 1382.7 1802.7 2342.3 3043.1

C.nnz/C.n 782.3 1048.4 1397.9 1875.6 2496.0 N/A

16 Ĉ.nnz/Ĉ.n 1548.9 2021.2 2628.5 3432.7 4472.1 N/A

C.nnz/C.n 1233.3 1679.4 2267.3 3075.3 N/A N/A

32 Ĉ.nnz/Ĉ.n 2841.0 3746.7 4911.0 6427.6 N/A N/A

C.nnz/C.n 1895.6 2597.7 3574.5 4896.9 N/A N/A

Ĉ.nnz/Ĉ.n 5149.9 6779.1 8979.8 11805.3 N/A N/A

C and reduces the largest matrix scale that can be evaluated. When choosing which input density to vary more aggressively, we increase A.nnz/A.n rather than B.nnz/B.n. This increases the number of distinct B rows accessed rather than the length of individual B rows, introducing additional indirection and more irregular memory accesses. For each scale, we evaluate all combinations of A and B for a total of 384 matrices. These RMat matrices are challenging because their intermediate products vary widely in size, the distribution of column indices spans the full column range of C, and their irregular structure limits the effectiveness of accumulators that exploit structural matrix locality. Table III shows a representative set of RMat matrices. For this set, B.nnz/B.n = 8 and A.nnz/A.n ∈ {4,8,16,32}. Due to the nonuniform structure that arises from the Graph500 parameters, the ratio Ĉ.nnz/C.nnz increases as A.nnz/A.n increases. This means that as A.nnz/A.n increases, the intermediate matrix grows faster than the final matrix, increasing the amount of accumulation work required per nonzero in the final matrix.

We evaluate gMAGNUS on two matrix data sets: the SuiteSparse matrix collection [28] and recursive model B. SpGEMM Evaluation power-law matrices (RMats) [47]. For SuiteSparse, we compute Table IV shows the geometric-mean (geomean) speedup of A2 , which is standard practice in SpGEMM evaluation. gMAGNUS over the five baselines across all datasets and Following prior work [27], [38], we include all matrices (330 in GPUs. The results demonstrate that gMAGNUS is faster than total) that require at least 100 million floating-point operations all baselines, with a geomean speedup ranging from 1.19× over and run to completion for gMAGNUS and at least one baseline OpSparse to 5.52× over TileSpGEMM. Failed baseline runs are on at least one GPU. We also highlight a subset of 20 matrices excluded from the geomean calculation. On H200, gMAGNUS (shown in Table II) with the largest C.nnz that satisfy two and cuSPARSE ran to completion for all 330 matrices, whereas additional criteria: they run to completion on both GPUs, and Kokkos, TileSpGEMM, and OpSparse completed 311, 271, at least 10% of rows are classified as heavy on both GPUs. and 295 matrices, respectively. On PVC, gMAGNUS and MKL For the RMat matrices, we use PaRMAT [29] with the standard ran to completion for 312 matrices (the remaining 18 excluded Graph500 parameters (a = 0.57, b = c = 0.19) to generate various because C did not fit in device memory), whereas Kokkos pairs of A and B. The matrix scale ranges from 15 to 20, completed 295 matrices. corresponding to matrices with 2scale rows and columns. The Figure 3 shows the gMAGNUS speedup (ratio of the baseline density of A varies as A.nnz/A.n ∈ {4,8,12,...,64}, while the time to gMAGNUS time) in log scale for the 20 representative density of B varies as B.nnz/B.n ∈ {4,8,12,16}. Increasing the matrices from Table II. gMAGNUS is fastest for 34 out of 40 density of either input matrix raises the memory cost of storing instances (20 matrices and 2 GPUs), with a geomean speedup of

gMAGNUS Speedup

↑

MKL

65.6

Kokkos

0.5

↑

165.3

cuSPARSE OpSparse

95.6

TileSpGEMM Kokkos

0.5 2 9 38 ley 90 ra44 rke pa ot0 6_ Be hd 33 d_ las C_ S for T n c H so Sta

blo

ya a0

ds ns

0.5 7 c-5

5 c2 c7 er_ 9_ inp lid b3 gG S_ an _F h F OP

bra TS

in-

t15 ne

-2 eu

p_ vs

tsk

slp

3 uth

cre

c-b

_ el1

d mo

tk1

_ ult

-ed

at2

raj

Fig. 3. Speedup in log scale for the 20 representative SuiteSparse matrices on Intel PVC (top) and NVIDIA H200 (bottom). The geomean is shown in the last group of bars. The ×-shaped markers denote failed runs. TABLE IV G EOMETRIC - MEAN SPEEDUP OF gMAGNUS OVER FIVE BASELINES .

SuiteSparse RMat Total

MKL PVC 1.28 3.33 1.98

Kokkos PVC H200 4.92 2.66 26.78 8.53 6.63

cuSPARSE H200 1.74 1.91 1.81

TileSpGEMM H200 5.52 11.43 7.62

OpSparse H200 1.19 7.94 3.29

TABLE V gMAGNUS RUNTIME PARAMETERS FOR THE 20 REPRESENTATIVE S UITE S PARSE MATRICES . Heavy IP %: FRACTION OF TOTAL INTERMEDIATE PRODUCT ELEMENTS FROM HEAVY ROWS . Heavy Row %: FRACTION OF HEAVY ROWS . Num. Chunks: NUMBER OF CHUNKS PER ROW.

para-9 Stanford Berkeley soc-Slashdot0902 HTC 336 4438 bloweya a0nsdsil c-57 brainpc2 TSOPF FS b39 c7 hangGlider 5 in-2004 net150 eu-2005 vsp south31 * vsp model1 * c-big pkustk12 mult dcop 02 wb-edu rajat28

Ĉ Size (GB) PVC H200 0.45 0.44 0.88 0.74 0.63 0.47 2.13 2.13 0.80 0.80 1.61 0.40 5.51 5.45 3.81 3.81 15.93 15.93 1.63 1.63 4.53 4.43 3.65 3.65 2.18 1.06 5.38 5.33 5.15 5.12 3.50 3.12 21.41 19.67 4.16 4.16 2.03 1.66 7.19 6.88

Heavy IP % PVC H200 0.21 0.21 0.50 0.42 0.67 0.51 0.98 0.98 0.99 0.99 0.99 0.25 0.99 0.98 0.99 0.99 0.99 0.99 0.99 0.99 0.33 0.32 0.98 0.98 0.32 0.16 0.99 0.98 0.96 0.96 0.92 0.82 0.99 0.92 0.99 0.99 0.16 0.13 0.99 0.94

Heavy Row % PVC H200 0.0457 0.0445 0.0038 0.0023 0.0894 0.0470 0.0550 0.0550 0.3334 0.3334 0.4375 0.0627 0.4372 0.4053 0.4999 0.4999 0.4986 0.4986 0.8888 0.8888 0.0068 0.0051 0.6273 0.6273 0.0381 0.0069 0.5802 0.5533 0.4430 0.4294 0.0766 0.0567 0.9937 0.9062 0.9042 0.9042 0.0013 0.0004 0.4052 0.3358

Num. Chunks PVC H200 64: 32, 2 32 256: 32, 8 128: 64, 2 32 16 64: 32, 2 32 8 4 32 16 16 8 8 4 8 4 4 2 512: 32, 16 256: 64, 4 16 8 256: 32, 8 128: 64, 2 16 8 16 8 128: 32, 4 64 32 16 8 4 4096: 32, 32, 4 2048: 64, 32 32 16

2.81×, 13.98×, 2.32×, 2.74×, and 14.51× over MKL, Kokkos, cuSPARSE, TileSpGEMM, and OpSparse, respectively. Besides TileSpGEMM, the geomean speedup across the representative set increased compared to the the full 330 matrices, highlighting the benefit of gMAGNUS on large matrices with many heavy rows. The speedup calculation excludes timeouts, for which TileSpGEMM had the highest rate, failing often for larger matrices. The vendor libraries were the most competitive, likely due to adaptive strategies that provide robust performance across a wide range of applications. OpSparse and Kokkos were the least competitive, demonstrating the drawbacks of TABLE VI global hash-map accumulators. This is most apparent for highly gMAGNUS RUNTIME PARAMETERS FOR A REPRESENTATIVE SET OF RM AT MATRICES , WITH A.nnz/A.n = 4 AND B.nnz/B.n = 8. Heavy IP %: irregular applications, such as Stanford Berkeley (web graph) FRACTION OF TOTAL INTERMEDIATE PRODUCT ELEMENTS FROM HEAVY ROWS . and the vsp matrices (random unweighted graphs), and for Heavy Row %: FRACTION OF HEAVY ROWS . Num. Chunks: NUMBER OF CHUNKS PER ROW. applications where the input nonzeros are distributed broadly across the column range, such as hangGlider 5. Matrix Ĉ Size (GB) Heavy IP % Heavy Row % Num. Chunks Scale PVC H200 PVC H200 PVC H200 PVC H200 Table V shows runtime parameters of gMAGNUS. The first 15 0.14 0.10 0.64 0.47 0.0487 0.0227 8 4 16 0.41 0.34 0.73 0.60 0.0599 0.0345 16 8 column shows the total size, in GB, of the CSR arrays that 17 1.17 0.99 0.80 0.69 0.0718 0.0424 32 16 18 3.24 2.91 0.86 0.77 0.0826 0.0538 64: 32, 2 32 compose Ĉ, which ranges from 0.47 to 21.41. The H200 values 19 8.87 8.10 0.90 0.82 0.0941 0.0613 128: 32, 4 64 are always lower because its larger local memory results in fewer 20 N/A 22.60 N/A 0.89 N/A 0.0745 N/A 128: 64, 2 heavy rows. Columns 2-3 show the fraction of intermediate product elements from heavy rows and the fraction of heavy rows, in Table IV. Table VI shows the runtime parameters used demonstrating that a small number of rows can produce large by gMAGNUS for A.nnz/A.n = 4 and B.nnz/B.n = 8. As intermediate products. For example, in Stanford Berkeley, less expected, the required number of chunks is tied solely to matrix than 1% of heavy rows account for half of the total intermediate scale, doubling with each scale increment. When the maximum product elements. The last column shows the required number number of chunks is reached (scale 18 for PVC, scale 20 for of chunks, notated as required number of chunks: level 0 number H200), multiple levels are employed. Additionally, the fraction of chunks, level 1 number of chunks, .... The required number of of heavy rows and heavy intermediate products increases with chunks varies significantly due to the wide range of matrix sizes, matrix scale, even though the densities of A and B are fixed, with 7 of these matrices requiring more than one level and one which further underscores the challenges of scaling SpGEMM matrix requiring 3 levels (wb-edu on PVC). In general, PVC to massive RMats. requires twice as many chunks due to its smaller local memory. Figure 4 shows the gMAGNUS speedup versus matrix For the RMat matrices, gMAGNUS demonstrates higher scale for A.nnz/A.n = {4,8,16,32}, which shows important speedups than for SuiteSparse, with a geomean speedup ranging performance trends as both matrix scale and A.nnz/A.n from 1.91× over cuSPARSE to 26.78× over Kokkos, as shown increase. On PVC, Kokkos times out for larger matrix scales

A.nnz/A.n = 4 64

A.nnz/A.n = 8 64

↑

A.nnz/A.n = 16

↑

A.nnz/A.n = 32

64 MKL

gMAGNUS Speedup

↑

124.6

66.3

0.5

88.4

0.5 215

216

217

218

219

220

1 0.5

0.5 215

216

217

218

219

220

215

216

217

218

219

215

220

0.5

217

218

Number of Rows

219

220

215

216

217

218

219

220

Number of Rows

217

218

219

220

cuSPARSE OpSparse

TileSpGEMM Kokkos

217

219

1 0.5 216

216

215

Kokkos

215

216

217

218

Number of Rows

219

220

215

216

218

220

Number of Rows

Fig. 4. Speedup in log scale vs number of rows of C for the representative matrices from the RMat matrix set on Intel PVC (top row) and NVIDIA H200 (bottom row). From left to right, the figures correspond to increasing A.nnz/A.n (nonzeros per row of A), where × markers denotes failed runs.

For the histogram kernel, we achieve approximately 40–60% of peak performance. For multisplit, which contributes the largest share of overall SpGEMM time among the three kernels (as shown in Figure 5), we achieve approximately 40–50% of peak performance. These results are consistent with [22], an in-depth study of the multisplit kernel on random input streams, which reported 30–60% of peak performance. In that study, the lower end of this range was observed for larger numbers of chunks, up to a maximum of 256. In contrast, gMAGNUS supports an unbounded number of chunks while still sustaining C. Performance of gMAGNUS Core Kernels near-optimal performance. For instance, wb-edu requires 2048 We evaluate the core kernels of gMAGNUS against their total chunks across two levels, yet still achieves 40% and 50% “speed of light”, i.e., their theoretical peak performance. The of peak performance for histogram and multisplit, respectively. theoretical peak time of a kernel is its theoretical minimum Figure 8 and Figure 7 show the same metrics for RMats data volume divided by the peak memory bandwidth, where with A.nnz/A.n = 4 and B.nnz/B.n = 8, providing insight the peak bandwidth is measured using simple memory copy into performance-critical kernels as we increase the matrix benchmarks (e.g., bandwidthTest included in CUDA). For dimensions. As the matrix scale increases, the percentage of example, the theoretical minimum data volume of the histogram peak performance converges to approximately 75%, 60%, and kernel (which computes Ĉreord .rowP tr) is: 50% for the outer product, histogram, and multisplit, respectively. Ĉ.nnz×sĈ.col +(Ĉ.n+1)×sĈ.rowP tr + Ĉreord .n×sĈreord .rowP tr . The phase breakdown shows that more time is spent in the (8) matrix core kernels as we increase the matrix scale due to the The first term represents reading Ĉ, the second term represents higher numbers of heavy rows. Additionally, the fraction of time reading Ĉ.rowP tr, and the third term represents writing to spent in the setup phase diminishes as matrix scale increases, Ĉreord .rowP tr. We only present results for H200, since PVC demonstrating that our setup overhead becomes negligible as the shows the same trends. heavy row intermediate product sizes increase. This convergence Figure 5 shows the breakdown of gMAGNUS into its demonstrates that our core kernels scale excellently for our performance-critical phases, and Figure 6 compares the core most irregular datasets, resulting in the excellent scaling of kernels against their speed-of-light performance. Because the gMAGNUS compared to other SpGEMM algorithms. fraction of heavy intermediate products varies across matrices, the V. C ONCLUSION time spent in each SpGEMM phase also varies widely. For example, para-9, soc-Slashdot0902, a0nsdsil, eu-2005, wb-edu, and ra- This paper presents gMAGNUS, a novel GPU algorithm for jat28 spend most of their time in the symbolic and numeric phases SpGEMM designed for massive irregular matrices. Such matrices because they have the smallest fraction of heavy intermediate often contain many heavy rows, whose intermediate products products, meaning that most intermediate products belong to light exceed local memory capacity and force conventional localrows. For the remaining matrices, heavy-row processing accounts memory accumulators to fall back to global-memory solutions. for a significant portion of runtime and often exceeds half of the The key idea behind gMAGNUS is to compute an intra-row total execution time. The outer product kernel achieves more than reordering of the intermediate products using an optimized outer 50% of theoretical peak performance for most matrices, and most product and a hierarchical multisplit, enabling accumulation of the largest matrices reach 60–80% of peak. The cases below entirely in local memory. We evaluate SYCL and CUDA 50% typically arise when dense columns of ACSC are multiplied implementations of gMAGNUS on Intel Ponte Vecchio and by highly sparse rows of B̃. This structure reduces coalesced NVIDIA H200, respectively, using matrices from the SuiteSparse writes to Ĉ and increases atomic updates to Ĉ.rowP tr, since collection and recursive power-law graphs. The results show that the number of atomic operations is proportional to ACSC .nnz. gMAGNUS outperforms five widely used baselines, including and is orders of magnitude slower when it completes. While MKL is competitive with gMAGNUS at smaller matrix scales, it becomes orders of magnitude slower beyond scale 16. On H200, all baselines except cuSPARSE are orders of magnitude slower than gMAGNUS, with TileSpGEMM and OpSparse failing for larger scales. cuSPARSE remains competitive for smaller sizes and lower values of A.nnz/A.n, but diverges from gMAGNUS as these parameters increase. Overall, gMAGNUS exhibits the best scaling as the dimensions and density scale up.

% of Total Time

1.2

Inter. Prod. Size

CSR-to-CSC

Outer Prod

Histo

Multisplit

Symbolic

Numeric

1 0.8 0.6 0.4 0.2 0

Fig. 5. Performance breakdown of gMAGNUS phases for the representative SuiteSparse matrices on H200. % of Theoretical Peak

Outer Prod

Histo

Multisplit

0.8 0.6 0.4 0.2 0 2 38 ley 90 44 rke ot0 6_ Be hd 33 las C_ S T c H so

ra-

rd_

nfo Sta

blo

il ds

7 c-5

ns a0

bra

c2 inp FS

_ PF

O TS

7 r_5 _c 39 ide Gl _b ng ha

0 20

in-

t1 ne

-20

lpt

uth

o _s

rew

el1

d mo

c-b

tk1

lt_

-ed

at2

raj

Fig. 6. Comparison of the core kernels in gMAGNUS to their theoretical upper bound for the representative SuiteSparse matrices on H200.

% of Total Time

1.5 Inter. Prod. Size Histo Numeric

CSR-to-CSC Multisplit

Outer Prod Symbolic

0.5

0 215

216

217 218 Number of Rows

219

220

Fig. 7. Performance breakdown of gMAGNUS phases for an RMat with A.nnz/A.n = 4 on H200. 1

% of Theoretical Peak

Outer Prod

Histo

Multisplit

0.8 0.6 0.4 0.2 0 215

216

217 218 Number of Rows

219

220

Fig. 8. Comparison of the core kernels in gMAGNUS to their theoretical upper bound for an RMat with A.nnz/A.n = 4 on H200.

cuSPARSE and MKL, with the largest gains on massive matrices with many heavy rows. In addition, its core kernels (outer product and multisplit) achieve near-optimal performance relative to the theoretical “speed of light” upper bound. R EFERENCES [1] J. Gao, W. Ji, F. Chang, S. Han, B. Wei, Z. Liu, and Y. Wang, “A systematic survey of general sparse matrix-matrix multiplication,” ACM Computing Surveys, vol. 55, no. 12, March 2023. [Online]. Available: https://doi.org/10.1145/3571157 [2] G. Guidi, O. Selvitopi, M. Ellis, L. Oliker, K. A. Yelick, and A. Buluç, “Parallel string graph construction and transitive reduction for de novo genome assembly,” in 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). Los Alamitos, CA, USA: IEEE Computer Society, May 2021, pp. 517–526. [Online]. Available: https://doi.ieeecomputersociety.org/10.1109/IPDPS49936.2021.00060 [3] K. Yelick, A. Buluç, M. Awan, A. Azad, B. Brock, R. Egan, S. Ekanayake, M. Ellis, E. Georganas, G. Guidi, S. Hofmeyr, O. Selvitopi, C. Teodoropol, and L. Oliker, “The parallelism motifs of genomic data analysis,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 378, no. 2166, p. 20190394, 01 2020.

[Online]. Available: https://doi.org/10.1098/rsta.2019.0394 [4] G. Guidi, M. Ellis, D. Rokhsar, K. Yelick, and A. Buluç, BELLA: Berkeley Efficient Long-Read to Long-Read Aligner and Overlapper, pp. 123–134. [Online]. Available: https://epubs.siam.org/doi/abs/10.1137/1.9781611976830.12 [5] X. Feng, Y. Xie, M. Song, W. Yu, and J. Tang, “Fast randomized PCA for sparse data,” in Asian conference on machine learning. PMLR, 2018, pp. 710–725. [Online]. Available: https://proceedings.mlr.press/v95/feng18a.html [6] O. Selvitopi, M. T. Hussain, A. Azad, and A. Buluç, “Optimizing high performance Markov clustering for pre-exascale architectures,” in 2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS), 2020, pp. 116–126. [Online]. Available: htpps://www.doi.org/10.1109/IPDPS47924.2020.00022 [7] E. Qin, A. Samajdar, H. Kwon, V. Nadella, D. Srinivasan, Sudarshan an Das, B. Kaul, and T. Krishna, “SIGMA: A sparse and irregular GEMM accelerator with flexible interconnects for DNN training,” in 2020 IEEE International Symposium on High Performance Computer Architecture (HPCA), 2020, pp. 58–70. [Online]. Available: htpps://www.doi.org/10.1109/HPCA47549.2020.00015 [8] A. Azad, G. A. Pavlopoulos, C. A. Ouzounis, N. C. Kyrpides, and A. Buluç, “HipMCL: a high-performance parallel implementation of the Markov clustering algorithm for large-scale networks,” Nucleic Acids Res., vol. 46, no. 6, pp. e33–e33, January 2018. [Online]. Available: htpps://www.doi.org/10.1093/nar/gkx1313 [9] T. Hoefler, D. Alistarh, T. Ben-Nun, N. Dryden, and A. Peste, “Sparsity in deep learning: pruning and growth for efficient inference and training in neural networks,” J. Mach. Learn. Res., vol. 22, no. 1, Jan. 2021. [Online]. Available: https://dl.acm.org/doi/abs/10.5555/3546258.3546499 [10] R. Li, B. Sjögreen, and U. M. Yang, “A new class of amg interpolation methods based on matrix-matrix multiplications,” SIAM Journal on Scientific Computing, vol. 43, no. 5, pp. S540–S564, 2021. [Online]. Available: https://doi.org/10.1137/20M134931X [11] R. Falgout, “An introduction to algebraic multigrid,” Computing in Science and Engineering, vol. 8, no. 6, pp. 24–33, 2006. [Online]. Available: https://doi.org/10.1109/MCSE.2006.105 [12] J. R. Gilbert, S. Reinhardt, and V. B. Shah, “High-performance graph algorithms from parallel sparse matrices,” in Applied Parallel Computing: State of the Art in Scientific Computing. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 260–269. [Online]. Available: https://doi.org/10.1007/978-3-540-75755-9 32 [13] V. Gleyzer, A. J. Soszynski, and E. K. Kao, “Leveraging linear algebra to count and enumerate simple subgraphs,” in 2020 IEEE High Performance Extreme Computing Conference (HPEC), 2020, pp. 1–8. [Online]. Available: https://doi.org/10.1109/HPEC43674.2020.9286191 [14] M. M. Wolf, J. W. Berry, and D. T. Stark, “A task-based linear algebra building blocks approach for scalable graph analytics,” in 2015 IEEE High Performance Extreme Computing Conference (HPEC), 2015, pp. 1–6. [Online]. Available: https://doi.org/10.1109/HPEC.2015.7322450 [15] A. Azad, A. Buluç, and J. Gilbert, “Parallel triangle counting and enumeration using matrix algebra,” in 2015 IEEE International Parallel and Distributed Processing Symposium (IPDPS) Workshop, 2015, pp. 804–811. [Online]. Available: https://doi.org/10.1109/IPDPSW.2015.75

[16] H. Kaplan, M. Sharir, and E. Verbin, “Colored intersection searching via sparse rectangular matrix multiplication,” in Proceedings of the Twenty-Second Annual Symposium on Computational Geometry (SCG). New York, NY, USA: Association for Computing Machinery, 2006, pp. 52–60. [Online]. Available: https://doi.org/10.1145/1137856.1137866 [17] J. Li, F. Wang, T. Araki, and J. Qiu, “Generalized sparse matrix-matrix multiplication for vector engines and graph applications,” in 2019 IEEE/ACM Workshop on Memory Centric High Performance Computing (MCHPC), 2019, pp. 33–42. [Online]. Available: https://doi.org/10.1109/MCHPC49590.2019.00012 [18] P. N. Q. Anh, R. Fan, and Y. Wen, “Balanced hashing and efficient GPU sparse general matrix-matrix multiplication,” in Proceedings of the 2016 International Conference on Supercomputing, ser. ICS16. New York, NY, USA: Association for Computing Machinery, 2016. [Online]. Available: https://doi.org/10.1145/2925426.2926273 [19] Z. Du, Y. Guan, T. Guan, D. Niu, L. Huang, H. Zheng, and Y. Xie, “OpSparse: A highly optimized framework for sparse general matrix multiplication on GPUs,” IEEE Access, vol. 10, pp. 85 960–85 974, 2022. [Online]. Available: https://doi.org/10.1109/ACCESS.2022.3196940 [20] M. Parger, M. Winter, D. Mlakar, and M. Steinberger, “spECK: accelerating GPU sparse matrix-matrix multiplication through lightweight analysis,” in Proceedings of the 25th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (PPoPP). New York, NY, USA: Association for Computing Machinery, 2020, pp. 362–375. [Online]. Available: https://doi.org/10.1145/3332466.3374521 [21] J. Wolfson-Pou, J. Laukemann, and F. Petrini, “Magnus: Generating data locality to accelerate sparse matrix-matrix multiplication on cpus,” in Proceedings of the 39th ACM International Conference on Supercomputing, ser. ICS ’25. New York, NY, USA: Association for Computing Machinery, 2025, p. 442–457. [Online]. Available: https://doi.org/10.1145/3721145.3725773 [22] S. Ashkiani, A. Davidson, U. Meyer, and J. D. Owens, “GPU multisplit: An extended study of a parallel algorithm,” ACM Transactions on Parallel Computing, vol. 4, no. 1, Aug. 2017. [Online]. Available: https://doi.org/10.1145/3108139 [23] Intel Corporation, Intel oneAPI Math Kernel Library Developer Reference for C, 2025. [24] M. Deveci, C. Trott, and S. Rajamanickam, “Multithreaded sparse matrix-matrix multiplication for many-core and GPU architectures,” Parallel Comput., vol. 78, no. C, p. 33–46, Oct. 2018. [Online]. Available: https://doi.org/10.1016/j.parco.2018.06.009 [25] S. Rajamanickam, S. Acer, L. Berger-Vergiat, V. Dang, N. Ellingwood, E. Harvey, B. Kelley, C. R. Trott, J. Wilke, and I. Yamazaki, “Kokkos kernels: Performance portable sparse/dense linear algebra and graph kernels,” 2021. [Online]. Available: https://arxiv.org/abs/2103.11991 [26] NVIDIA Corporation, cuSPARSE Library, 2025. [27] Y. Niu, Z. Lu, H. Ji, S. Song, Z. Jin, and W. Liu, “TileSpGEMM: a tiled algorithm for parallel sparse general matrix-matrix multiplication on GPUs,” in Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (PPoPP). New York, NY, USA: Association for Computing Machinery, 2022, pp. 90–106. [Online]. Available: https://doi.org/10.1145/3503221.3508431 [28] T. A. Davis and Y. Hu, “The university of florida sparse matrix collection,” ACM Transactions on Mathematical Software, vol. 38, no. 1, December 2011. [Online]. Available: https://doi.org/10.1145/2049662.2049663 [29] F. Khorasani, R. Gupta, and L. N. Bhuyan, “Scalable SIMD-efficient graph processing on GPUs,” in 2015 International Conference on Parallel Architecture and Compilation (PACT), 2015, pp. 39–50. [Online]. Available: https://doi.org/10.1109/PACT.2015.15 [30] F. G. Gustavson, “Two fast algorithms for sparse matrices: Multiplication and permuted transposition,” ACM Transactions on Mathematical Software, vol. 4, no. 3, pp. 250–269, September 1978. [Online]. Available: https://doi.org/10.1145/355791.355796 [31] S. Dalton, L. Olson, and N. Bell, “Optimizing sparse matrix-matrix multiplication for the GPU,” ACM Transactions on Mathematical Software, vol. 41, no. 4, October 2015. [Online]. Available: https://doi.org/10.1145/2699470 [32] Z. Gu, J. Moreira, D. Edelsohn, and A. Azad, “Bandwidth optimized parallel algorithms for sparse matrix-matrix multiplication using propagation blocking,” in Proceedings of the 32nd ACM Symposium on Parallelism in Algorithms and Architectures, ser. SPAA ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 293–303. [Online]. Available: https://doi.org/10.1145/3350755.3400216 [33] S. Pal, J. Beaumont, D. hyeon Park, A. Amarnath, S. Feng, C. Chakrabarti, H.-S. Kim, D. Blaauw, T. N. Mudge, and R. G. Dreslinski, “OuterSPACE: An outer product based sparse matrix multiplication accelerator,” 2018 IEEE International Symposium on High Performance Computer Architecture (HPCA), pp. 724–736, 2018. [Online]. Available: https://api.semanticscholar.org/CorpusID:4571588

[34] J. Lee, S. Kang, Y. Yu, Y.-Y. Jo, S.-W. Kim, and Y. Park, “Optimization of GPU-based sparse matrix multiplication for large sparse networks,” in 2020 IEEE 36th International Conference on Data Engineering (ICDE), 2020, pp. 925–936. [Online]. Available: https://doi.org/10.1109/ICDE48307.2020.00085 [35] Z. Zhang, H. Wang, S. Han, and W. J. Dally, “Sparch: Efficient architecture for sparse matrix multiplication,” 2020 IEEE International Symposium on High Performance Computer Architecture (HPCA), pp. 261–274, 2020. [Online]. Available: https://api.semanticscholar.org/CorpusID:211205022 [36] Y.-Y. Jo, S.-W. Kim, and D.-H. Bae, “Efficient sparse matrix multiplication on GPU for large social network analysis,” in Proceedings of the 24th ACM International on Conference on Information and Knowledge Management, ser. CIKM ’15. New York, NY, USA: Association for Computing Machinery, 2015, p. 1261–1270. [Online]. Available: https://doi.org/10.1145/2806416.2806445 [37] Y. Nagasaka, A. Nukada, and S. Matsuoka, “High-performance and memorysaving sparse general matrix-matrix multiplication for NVIDIA Pascal GPU,” in 2017 46th International Conference on Parallel Processing (ICPP), 2017, pp. 101–110. [Online]. Available: https://doi.org/10.1109/ICPP.2017.19 [38] M. Wu, H. Luo, F. Li, Y. Zhang, Z. Tang, K. Li, J. Zhang, and C. Liu, “HSMU-SpGEMM: Achieving high shared memory utilization for parallel sparse general matrix-matrix multiplication on modern GPUs,” in 2025 IEEE International Symposium on High Performance Computer Architecture (HPCA), 2025, pp. 1452–1466. [Online]. Available: https://doi.org/10.1109/HPCA61900.2025.00109 [39] S. Dalton, S. Baxter, D. Merrill, L. Olson, and M. Garland, “Optimizing sparse matrix operations on GPUs using merge path,” in 2015 IEEE International Parallel and Distributed Processing Symposium, 2015, pp. 407–416. [Online]. Available: https://doi.org/10.1109/IPDPS.2015.98 [40] S. Dalton, N. Bell, L. Olson, and M. Garland, “Cusp: Generic parallel algorithms for sparse matrix and graph computations,” 2026, version 0.6.0. [Online]. Available: https://github.com/cusplibrary/cusplibrary [41] M. Winter, D. Mlakar, R. Zayer, H.-P. Seidel, and M. Steinberger, “Adaptive sparse matrix-matrix multiplication on the GPU,” in Proceedings of the 24th Symposium on Principles and Practice of Parallel Programming (PPoPP). New York, NY, USA: Association for Computing Machinery, 2019, pp. 68–81. [Online]. Available: https://doi.org/10.1145/3293883.3295701 [42] W. Liu and B. Vinter, “A framework for general sparse matrix-matrix multiplication on GPUs and heterogeneous processors,” J. Parallel Distrib. Comput., vol. 85, no. C, p. 47–61, Nov. 2015. [Online]. Available: https://doi.org/10.1016/j.jpdc.2015.06.010 [43] F. Gremse, A. Höfter, L. O. Schwen, F. Kiessling, and U. Naumann, “GPU-accelerated sparse matrix-matrix multiplication by iterative row merging,” SIAM Journal on Scientific Computing, vol. 37, no. 1, pp. C54–C71, 2015. [Online]. Available: https://doi.org/10.1137/130948811 [44] O. Zachariadis, N. Satpute, J. Gómez-Luna, and J. Olivares, “Accelerating sparse matrix–matrix multiplication with GPU tensor cores,” Computers & Electrical Engineering, vol. 88, p. 106848, 2020. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045790620307011 [45] N. Bell, S. Dalton, and L. N. Olson, “Exposing fine-grained parallelism in algebraic multigrid methods,” SIAM Journal on Scientific Computing, vol. 34, no. 4, pp. C123–C152, 2012. [Online]. Available: https://doi.org/10.1137/110838844 [46] C. Nugteren, G.-J. van den Braak, H. Corporaal, and B. Mesman, “High performance predictable histogramming on gpus: exploring and evaluating algorithm trade-offs,” in Proceedings of the Fourth Workshop on General Purpose Processing on Graphics Processing Units, ser. GPGPU-4. New York, NY, USA: Association for Computing Machinery, 2011. [Online]. Available: https://doi.org/10.1145/1964179.1964181 [47] D. Chakrabarti, Y. Zhan, and C. Faloutsos, “R-MAT: A recursive model for graph mining,” in Proceedings of the 2004 SIAM International Conference on Data Mining (SDM). SIAM, 2004, pp. 442–446. [Online]. Available: https://doi.org/10.1137/1.9781611972740.43

Record · ID 405622 · SHA-256 7c5a9f05cc331aa8
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.