FlashSparse: Efficient Sparse Matmul on GPUs

• FlashSparse is a Tensor-Core-oriented method that accelerates sparse matrix multiplication and SDDMM by reducing sparse granularity from 16×1 to 8×1.
• It employs a swap-and-transpose strategy to remap sparse operands, significantly cutting compute redundancy and enhancing memory coalescing.
• The approach utilizes a compact ME-BCRS storage format and a memory-efficient thread mapping to deliver substantial performance gains in GNN workloads and scientific computing.

FlashSparse is a Tensor-Core-oriented method for accelerating sparse matrix-matrix multiplication (SpMM) and sampled dense-dense matrix multiplication (SDDMM) on GPUs under highly irregular, unstructured sparsity. It is designed to close the gap between sparse workloads and Tensor Core Units (TCUs), with the central claim that prior TCU-based methods waste compute and bandwidth because their sparse packing granularity is too coarse. FlashSparse addresses this by combining a swap-and-transpose matrix multiplication strategy, a memory-efficient thread mapping strategy for coalesced data access, and a compact sparse storage format, ME-BCRS (Shi et al., 2024).

1. Problem domain and motivating inefficiency

FlashSparse targets two sparse operators that are described as important in scientific computing and deep learning. In the paper’s notation for SpMM, the operation is

$C = AB,$

with elementwise form

$C_{ij} = \sum_k A_{ik} B_{kj}.$

For SDDMM, the output is sparse and sampled from a dense-dense multiplication. The paper situates these operators in graph neural network workloads: in GCNs, neighborhood aggregation is an SpMM, while in attention-based GNNs such as AGNN and GAT, attention-score generation is often expressed through SDDMM and then consumed by SpMM (Shi et al., 2024).

The architectural problem is a mismatch between unstructured sparsity and the rigid tile shapes of Tensor Core MMA instructions. The paper lists representative instruction shapes such as WMMA TF32 $m16n16k8$, MMA TF32 $m16n8k4$ and $m16n8k8$, and MMA FP16 $m16n8k8$ and $m16n8k16$. Prior Tensor Core methods are described as packing sparse data into units aligned with the $m=16$ dimension, which yields $16\times 1$ nonzero-vector granularity. Because sparse matrices are highly irregular, many of these vectors contain only a few true nonzeros and many padded zeros. The paper defines the resulting excess arithmetic as computation redundancy: extra MMA invocations and multiply-adds on structural zeros that arise solely from Tensor Core packing (Shi et al., 2024).

This redundancy is quantified directly. On real graph datasets, the number of zeros inside such “nonzero vectors” is reported to exceed the true nonzeros by $5.6\times$ to $C_{ij} = \sum_k A_{ik} B_{kj}.$0. In one SpMM case with dense matrix column count $C_{ij} = \sum_k A_{ik} B_{kj}.$1, reducing the granularity from $C_{ij} = \sum_k A_{ik} B_{kj}.$2 to $C_{ij} = \sum_k A_{ik} B_{kj}.$3 lowers the number of MMA invocations by an average of 43%. The motivating claim is therefore not merely that sparse workloads are irregular, but that the prevailing sparse-to-Tensor-Core mapping itself is the dominant source of inefficiency (Shi et al., 2024).

2. Swap-and-transpose strategy and minimum sparse granularity

The central technical idea in FlashSparse is the identity

$C_{ij} = \sum_k A_{ik} B_{kj}.$4

The paper uses this apparently simple transformation to change which MMA dimension is associated with the sparse operand. In prior mappings, the sparse matrix is treated as the left operand and therefore constrained by the minimum $C_{ij} = \sum_k A_{ik} B_{kj}.$5-dimension of the MMA shape, which is 16. FlashSparse instead swaps operands and performs the computation on transposed fragments, so that the sparse operand is mapped onto the smaller $C_{ij} = \sum_k A_{ik} B_{kj}.$6-dimension, which is 8. This is what enables the reduction of sparse granularity from $C_{ij} = \sum_k A_{ik} B_{kj}.$7 to $C_{ij} = \sum_k A_{ik} B_{kj}.$8 while still using native Tensor Core instructions (Shi et al., 2024).

The paper describes the method as consisting of four coordinated steps: operands swapping, transposed access, transposed computation, and transposed output. For SpMM, conventional FP16 Tensor Core mapping uses sparse left-operand blocks of size $C_{ij} = \sum_k A_{ik} B_{kj}.$9, dense right-operand blocks of size $m16n16k8$0, and outputs of size $m16n16k8$1. Under FlashSparse, the sparse TC block becomes $m16n16k8$2 for FP16 or $m16n16k8$3 for TF32, while the dense TC block is $m16n16k8$4. The actual register-level MMA is then performed on $m16n16k8$5 and $m16n16k8$6, and the final result is written back through a transposed layout-compatible path (Shi et al., 2024).

The paper argues that this achieves the minimum sparse granularity available under the targeted instruction set, because the sparse packing is now aligned with the smaller MMA dimension. That reduction is directly tied to higher useful density inside packed sparse TC blocks. The illustrative example given in the paper shows a case that requires 4 MMAs with $m16n16k8$7 vectors but only 2 MMAs under FlashSparse’s $m16n16k8$8 packing. A plausible implication is that FlashSparse’s novelty lies less in a new sparse kernel primitive than in an algebraic reorientation of the multiplication that changes the hardware-imposed sparsity unit itself (Shi et al., 2024).

3. Data representation, memory access, and kernel organization

FlashSparse organizes a sparse matrix into vertical vectors and row-wise windows. A vector of size $m16n16k8$9 is retained if it contains at least one nonzero; such retained units are termed nonzero vectors. Within each window, every $m16n8k4$0 nonzero vectors are grouped into a sparse TC block, where $m16n8k4$1 depends on the instruction shape: for FP16, $m16n8k4$2; for TF32 under $m16n8k4$3, $m16n8k4$4. These packed structures are then stored in ME-BCRS, a custom sparse format comprising three arrays: RowPointers, ColumnIndices, and Values (Shi et al., 2024).

RowPointers stores the start index of each row window in ColumnIndices. ColumnIndices stores the column indices of the nonzero vectors. Values stores the packed sparse TC-block values in row-major order. The row-major choice is explicitly tied to the operand-access requirements induced by swap-and-transpose execution. The format differs from padding-based blocked schemes because FlashSparse does not physically store padded zero vectors in the last TC block of a window. Instead, the kernel reconstructs whether a thread corresponds to a valid vector by computing

$m16n8k4$5

and injecting zeros into registers when the thread would otherwise read beyond the valid vectors in the final block (Shi et al., 2024).

A second major systems contribution is the memory-efficient thread mapping strategy. Under a direct mapping, the dense TC block $m16n8k4$6 needed for SpMM is loaded in a way that causes groups of 8 threads to touch only 16 bytes per phase, even though NVIDIA memory transactions are at least 32 bytes. The paper states that this requires 16 memory transactions to load the dense TC block. FlashSparse instead shuffles columns in registers so that each thread loads a local $m16n8k4$7 region rather than two far-apart columns. In the FP16 example, this reduces the load to 8 transactions instead of 16, a 50% reduction, and the same mapping is reused for output because the register layouts of $m16n8k4$8 and $m16n8k4$9 are consistent (Shi et al., 2024).

The kernel organization is warp-centric. Each warp is assigned sparse TC blocks and the corresponding output tile region. The sparse block is loaded from ME-BCRS, the required rows of the dense matrix are gathered using the vector column indices, the fragments are transposed in registers to form the left and right MMA operands, and partial results are accumulated over successive sparse TC blocks from the same window. The implementation is therefore built around reducing both arithmetic redundancy and dense-row gather overhead, rather than treating storage savings alone as the main optimization objective (Shi et al., 2024).

4. Operator workflows for SpMM and SDDMM

For both operators, FlashSparse begins with a preprocessing stage that translates the sparse structure into TCU-oriented packed blocks. The paper states that this CSR-to-ME-BCRS translation is parallelized on the GPU using CUDA and is one-time when sparsity is static. In end-to-end GNN training with static graphs, the reported preprocessing overhead is less than 1% of total runtime, which frames the method as particularly suited to workloads with reusable sparse structure (Shi et al., 2024).

In SpMM, a warp first loads a sparse TC block $m16n8k8$0 from ME-BCRS. Under FlashSparse, $m16n8k8$1 is $m16n8k8$2 for FP16 or $m16n8k8$3 for TF32. Using the stored vector column indices, the warp gathers the corresponding rows from dense matrix $m16n8k8$4, producing a dense TC block of shape $m16n8k8$5. The sparse and dense blocks are then transposed in registers to become $m16n8k8$6 and $m16n8k8$7, and Tensor Core MMA accumulates

$m16n8k8$8

After accumulation over all relevant sparse TC blocks, the transposed output is written back through the coalesced mapping (Shi et al., 2024).

In SDDMM, the same swap-and-transpose principle is applied to dense inputs with sparse output. The paper notes that the data layouts are especially favorable here because one dense input is row-major and the other is column-major, matching the swapped transposed access pattern. The output sparse TC block generated by FlashSparse is $m16n8k8$9, whereas downstream SpMM expects sparse TC block shapes of $m16n8k8$0 for FP16 or $m16n8k8$1 for TF32. FlashSparse therefore splits the SDDMM output into sub-blocks before storage so that it can be consumed efficiently by subsequent SpMM kernels. Algorithm 1 in the paper gives explicit per-thread offset formulas for the FP16 $m16n8k8$2 and TF32 $m16n8k8$3 output sub-blocks, showing that the output transformation is treated as a first-class part of the kernel design rather than a generic postprocessing step (Shi et al., 2024).

This operator-specific differentiation is important. SpMM starts from sparse input and produces dense output, while SDDMM starts from dense inputs and produces sparse output. FlashSparse’s contribution is to place both within a single Tensor-Core-oriented framework based on $m16n8k8$4 granularity and transposed MMA execution, while adapting the storage and writeback path to the different operator semantics (Shi et al., 2024).

5. Experimental results and reported performance

The evaluation covers 515 sparse matrices: 500 representative sparse matrices from a SuiteSparse-selected corpus satisfying more than 10k rows, more than 10k columns, and more than 100k nonzeros, plus 15 real graph datasets including GitHub, Artist, Blog, Ell, Yelp, DD, Reddit, Amazon, Amazon0505, Comamazon, Yeast, OGBProducts, AmazonProducts, IGB-small, and IGB-medium. The experiments use NVIDIA H100 PCIe and NVIDIA GeForce RTX 4090 GPUs, with FlashSparse integrated into PyTorch. Baselines include CUDA-core sparse kernels such as RoDe, Sputnik, GNNAdvisor, GE-SpMM, and cuSPARSE, as well as Tensor Core methods DTC-SpMM and TC-GNN (Shi et al., 2024).

The reported speedups are summarized below.

Workload	Platform	Reported result
SpMM	RTX 4090	5.5× geomean over DTC-SpMM; 3.22× over RoDe
SpMM	H100	4.41× mean over DTC-SpMM; 2.23× mean over RoDe
SDDMM	H100	2.92× geomean over RoDe
SDDMM	RTX 4090	2.18× geomean over RoDe
GCN end-to-end	RTX 4090	1.57× geomean over latest DGL
AGNN end-to-end	RTX 4090	1.79× geomean over latest DGL

For SpMM on RTX 4090, FlashSparse reports a geometric mean speedup of 5.5× over DTC-SpMM and 3.22× over RoDe, with maximum speedups of 25.26× and 14.2×, respectively. On H100 for $m16n8k8$5, the corresponding mean speedups are 4.41× over DTC-SpMM and 2.23× over RoDe, with maxima of 16.03× and 6.9×. The paper also reports geometric mean throughput on RTX 4090 of 4888 GFLOPS for FP16 and 2697 GFLOPS for TF32, with peaks up to 26 TFLOPS FP16 and 16 TFLOPS TF32 (Shi et al., 2024).

Distributional results are similarly strong. On RTX 4090, FlashSparse is reported to beat DTC-SpMM on 98.8% of matrices by at least 2×, to beat RoDe by at least 2× on 80.16% of matrices, and to beat Sputnik by at least 2× on 94.59% of matrices. On H100, it beats DTC-SpMM by at least 2× on 99.01% of matrices and RoDe by at least 2× on 50.2% of matrices. For SpMM with $m16n8k8$6, FP16, the $m16n8k8$7 granularity reduces data-access cost by up to 49%, with an average reduction of 35% relative to $m16n8k8$8 packing (Shi et al., 2024).

For SDDMM, FlashSparse reports a geometric mean speedup of 2.92× over RoDe on H100 and 2.18× on RTX 4090, with maxima of 18.59× and 14.93×. For $m16n8k8$9, FlashSparse is at least 2× faster than RoDe on 68.34% of matrices on H100 and 49.79% on RTX 4090. The corresponding data-access reduction for FP16 SDDMM reaches 49% with an average of 28% relative to $m16n8k16$0 granularity (Shi et al., 2024).

The ablation studies attribute the gains to three distinct mechanisms. First, replacing $m16n8k16$1 with $m16n8k16$2 while holding other factors fixed yields 1.89× geomean speedup for SpMM and 2.61× for SDDMM on H100. Second, the coalesced thread mapping contributes 1.34× average speedup on H100 and 1.18× on RTX 4090, with maxima up to 2.0×. Third, ME-BCRS reduces memory footprint by 11.72% on average over SR-BCRS, with a maximum reduction of 50.0%, and 336/515 matrices showing more than 10% footprint reduction (Shi et al., 2024).

FlashSparse is also evaluated end-to-end in GNN training. On RTX 4090, relative to the latest DGL, the paper reports 1.57× geometric mean speedup for GCN and 1.79× for AGNN, with maximum end-to-end speedups of 1.8× and 2.83×. The reported end-to-end time includes format translation, forward pass, backward pass, and optimizer/model update. For GCN on several datasets, FlashSparse in FP16 and TF32 is reported to show accuracy comparable to FP32 frameworks, with no reported accuracy loss (Shi et al., 2024).

6. Scope, limitations, and relation to adjacent systems

FlashSparse is presented as preferable when the workload is SpMM or SDDMM, the sparse pattern is unstructured and highly sparse, the platform offers strong Tensor Core performance, and sparsity is static or slowly changing so preprocessing can be amortized. The method is explicitly implemented for TF32 and FP16 on Tensor Cores, and the paper frames graph neural networks, graph analytics, and some scientific sparse linear algebra as natural application domains. It is especially suited to cases where irregular sparsity makes $m16n8k16$3-row grouping inefficient and where dense feature dimensions are large enough for Tensor Core throughput to dominate (Shi et al., 2024).

The paper also makes its constraints clear. The approach is specialized to Tensor Core MMA shapes with unbalanced dimensions and relies on mapping the sparse operand to the smaller $m16n8k16$4 dimension. Its scope is limited to SpMM and SDDMM rather than arbitrary sparse tensor contractions. Although preprocessing is cheap in the reported static-sparsity regime, the design is less obviously advantageous when the sparse structure changes frequently. The reported comparisons are against FP32 CUDA-core baselines and TF32 or FP16 Tensor Core baselines; FlashSparse is not presented as a generic FP32 sparse kernel. The paper further notes that for very regular or hardware-native sparse regimes, such as block-structured sparsity or 2:4 sparsity, other approaches such as cuSPARSELt may be preferable (Shi et al., 2024).

Within the broader sparse-systems landscape, FlashSparse occupies a specific design point. It differs from SPLAT, which targets regular sparse multi-head self-attention on GPUs through the affine-compressed-sparse-row format and JIT code generation for R-SDDMM and R-SpMM under static regular masks (Gupta et al., 2024). It also differs from work on sparse transformer MLP execution such as “Sparser, Faster, Lighter Transformer LLMs,” which focuses on unstructured activation sparsity in feedforward layers through tile-wise ELLPACK packing and fused sparse kernels on modern NVIDIA GPUs (Cetin et al., 24 Mar 2026). A plausible implication is that FlashSparse should be understood as a sparse linear algebra method specialized for Tensor Core execution of SpMM and SDDMM, rather than as a general sparse-attention framework or a sparse-MLP runtime. The title may also invite confusion with flash-storage systems, but that is a distinct line of work; for example, BigSparse is a fully external graph analytics system built around SSD arrays and sort-reduce-based sequentialization of random vertex updates, not a Tensor Core sparse multiplication method (Jun et al., 2017).

The enduring significance of FlashSparse lies in the paper’s claim that the decisive optimization target is not only sparse storage or load balancing, but the granularity at which sparsity is exposed to Tensor Core MMA. By changing the algebraic orientation of the matrix product so that sparsity aligns with the smaller Tensor Core dimension, FlashSparse reframes unstructured sparse acceleration as a problem of operand placement within the hardware’s fixed fragment geometry (Shi et al., 2024).