100 DAYS OF GPU Course preview
Day 37Module 4
in-technical-review

Sparse matrices in CUDA: COO, CSR and SpMV

Here are the first 32 rows of two sparse matrices, written as row lengths.

even    16 16 16 16 16 16 16 16 ... 16 16 16 16       512 nonzeros
skewed   3  3  3  3  3  3  3  3 ...  3  3  3 419      512 nonzeros

Same dimensions, same number of nonzeros, the same rule for columns and the same rule for values. Put one thread on each row, which is the decomposition many programmers write first, and those 32 rows form one warp. On the top matrix, that warp runs 16 loop iterations.

On the bottom matrix it runs 419 iterations because a warp stops only when its slowest lane stops. The other 31 lanes finish after three iterations. The same kernel and arithmetic now issue 26 times as many lane slots.

This is the first lesson where the shape of the data, and not the code, decides the performance. By the end you can look at a row-length histogram and say which of two SpMV kernels your matrix wants.

Three arrays, and what dropping one of them buys

Store a 524,288 by 524,288 matrix densely and it is 2^40 bytes, one tebibyte, whatever is in it. With 8,388,608 nonzeros in it, the interesting part is 64 mebibytes. Every sparse format is a way of writing down only that part, plus enough bookkeeping to find it again.

COO, coordinate format, is the obvious one: three arrays, each of length nnz, holding a row index, a column index and a value per nonzero. cuSPARSE describes exactly that, "the pointers to the row indices array of length nnz", and the same for columns and values (https://docs.nvidia.com/cuda/cusparse/index.html section 3.3.2, checked 2026-08-30). You can build it in any order, append to it, and read it without knowing anything about the matrix.

CSR, compressed sparse row, keeps the columns and values and throws the row-index array away. In its place goes one offset per row: "the row offsets array of length number of rows + 1 that represents the starting position of each row in the columns and values arrays" (same page, section 3.3.3). Row r owns everything in [rowPtr[r], rowPtr[r + 1]).

That drops 32 mebibytes of row indices and adds 2 mebibytes of offsets. More important, a row is now a contiguous slice one thread or one warp can own outright. In COO nothing tells you where a row ends, so two threads can hold nonzeros belonging to the same output element and the accumulate has to be atomic:

y[rowIdx[e]] += val[e] * x[colIdx[e]];            // a race, silently wrong
atomicAdd(&y[rowIdx[e]], val[e] * x[colIdx[e]]);  // correct, and it queues

Day 26 measured what that queue costs when many threads target the same address, and the heavy row below aims 419 of them at one float. CSR removes the atomic by giving each output a single owner, which is why every SpMV kernel on this page is a CSR kernel.

Diagram: one 32-row group, 512 nonzeros, three ways. Three horizontal bands. Each band is 32 lane boxes stacked vertically, with a bar to the right of each lane whose length is the number of loop iterations that lane runs. Filled bar means the lane had a nonzero that iteration, hollow means it was along for the ride. Band 1, one thread per row on the even matrix: 32 bars of length 16, all filled. Caption "512 lane-slots issued, 512 of them do work." Band 2, one warp per row on the skewed matrix: 32 separate warps handle the 32 rows, 31 of them taking a single 32-wide pass with 3 lanes filled and one taking fourteen passes. Caption "1,440 lane-slots issued, 512 do work." Band 3, one thread per row on the skewed matrix: 31 bars of length 3 and one bar of length 419, with the 31 short bars followed by hollow boxes out to 419. Caption "13,408 lane-slots issued, 512 do work." Alt text: "For one 32-row group with 512 nonzeros, the three mappings issue 512, 1,440, and 13,408 lane-slots, yielding 100, 35.6, and 3.8 percent utilization."

The row is not the unit of work you think it is

One thread per row is the decomposition every CPU programmer brings. It is correct, it needs no atomics and no reduction, and each thread writes exactly one output. On a machine with 32 independent cores it would also be the right answer.

The 32 lanes of a warp are not 32 independent workers. They share one instruction stream, so a per-lane loop bound is not a shortcut for 31 of them, it is a wait. That is divergence again, which day 22 measured on a branch; here it arrives as a trip count and the cost is set by the longest row in the warp, not the mean.

The addresses are the second problem, even when all rows have the same length. At step k, lane L reads val[rowPtr[r0 + L] + k]. On a matrix with 16 nonzeros per row, those addresses are 16 floats apart.

The warp's 32 lanes land in 32 different 32-byte sectors, so one request pulls 32 sectors. Day 11 measured that stride as the lowest global memory throughput in its test.

Bell and Garland named both kernels and the trade between them in 2008, and nothing since has changed the shape of the answer: give the row to a warp instead, and each lane takes every 32nd nonzero, so the 32 addresses are coalesced again (Efficient Sparse Matrix-Vector Multiplication on CUDA, NVIDIA Technical Report NVR-2008-004, https://research.nvidia.com/publication/2008-12_efficient-sparse-matrix-vector-multiplication-cuda , checked 2026-08-30). The imbalance does not vanish, it moves from between the lanes of one warp to between warps, where the scheduler can hide it behind other resident work (latency hiding). The cost is that a row of 3 nonzeros still occupies a full 32-wide pass and a five-step shuffle reduction.

Holding the code still and changing the data

Full program in code/day37-spmv/spmv.cu. The kernels are short. The discipline around them is the lesson.

Two matrices, one difference. Both are 524,288 square, both hold 8,388,608 nonzeros, both use the same rule for column indices and the same rule for values. Only the row lengths change: 16 everywhere, or 31 rows of 3 and one row of 419 in each aligned group of 32 rows. A group of 32 rows is exactly one warp of the thread-per-row kernel, so the imbalance lands inside a warp rather than being averaged away across the grid.

The heavy row length is derived, not typed. It is 32 * 16 - 31 * 3, and a static_assert fails the build if the two matrices ever stop holding the same number of nonzeros. A timing table that compares matrices with different nonzero counts is the bug day 11 opens with, one level up.

Lane-slots are computed before anything runs. A lane-slot is one lane of one warp on one loop iteration, whether or not that lane has a nonzero. It is host arithmetic over the row lengths, printed above the timings, and it is the prediction the clock gets judged against.

A floor, measured in the same process. A third kernel walks the nonzeros in a grid-stride loop with no row structure at all, so its 32 lanes read 32 consecutive values. Its output is not a matrix-vector product, which is why it is a probe, but its time is the ceiling neither SpMV kernel can beat. SpMV does two flops per eight bytes of matrix, so it is memory bound on any card and that ceiling is a bandwidth number.

The conversion that turns COO into CSR is three passes and no sorting:

    csr->rowPtr.assign(kRows + 1, 0);
    for (size_t e = 0; e < cooRow.size(); ++e) {
        ++csr->rowPtr[static_cast<size_t>(cooRow[e]) + 1];
    }
    for (size_t r = 0; r < kRows; ++r) {
        csr->rowPtr[r + 1] += csr->rowPtr[r];
    }
    std::vector<int> cursor(csr->rowPtr.begin(), csr->rowPtr.end() - 1);
    csr->colIdx.assign(cooCol.size(), 0);
    csr->val.assign(cooVal.size(), 0.0f);
    for (size_t e = 0; e < cooRow.size(); ++e) {
        const size_t row = static_cast<size_t>(cooRow[e]);
        const size_t slot = static_cast<size_t>(cursor[row]++);
        csr->colIdx[slot] = cooCol[e];
        csr->val[slot] = cooVal[e];
    }

Count, prefix sum, scatter. The middle loop is a serial scan, which is the same primitive day 31 puts on the GPU.

Then the kernel everybody writes first:

__global__ void spmvCsrThreadPerRow(const int* __restrict__ rowPtr,
                                    const int* __restrict__ colIdx,
                                    const float* __restrict__ val,
                                    const float* __restrict__ x,
                                    float* __restrict__ y, size_t nRows) {
    const size_t row =
        blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    if (row < nRows) {
        float sum = 0.0f;
        for (int j = rowPtr[row]; j < rowPtr[row + 1]; ++j) {
            sum += val[j] * x[colIdx[j]];
        }
        y[row] = sum;
    }
}

And the same product with the row given to a warp:

__global__ void spmvCsrWarpPerRow(const int* __restrict__ rowPtr,
                                  const int* __restrict__ colIdx,
                                  const float* __restrict__ val,
                                  const float* __restrict__ x,
                                  float* __restrict__ y, size_t nRows) {
    const size_t t = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    const size_t row = t / kWarpSize;
    const unsigned int lane = threadIdx.x % kWarpSize;
    if (row < nRows) {
        const int end = rowPtr[row + 1];
        float sum = 0.0f;
        for (int j = rowPtr[row] + static_cast<int>(lane); j < end;
             j += kWarpSize) {
            sum += val[j] * x[colIdx[j]];
        }
        for (int offset = kWarpSize / 2; offset > 0; offset /= 2) {
            sum += __shfl_down_sync(0xffffffffu, sum, offset);
        }
        if (lane == 0) {
            y[row] = sum;
        }
    }
}

The guard is warp-uniform on purpose. A warp's 32 threads hold 32 consecutive values of t, so t / 32 is the same for all of them. Every lane named in the shuffle mask reaches the shuffle.

A block size that is not a whole number of warps, or a bound computed per lane, breaks that rule. Day 23 covers the result.

Note. This matrix is banded and synthetic. Its columns run consecutively from the diagonal, so the gather into x is coalesced for both kernels, and x is 2 mebibytes against the T4's 4 mebibytes of L2. A real matrix scatters those reads across a vector that does not fit, and that cost is missing here. Holding the column pattern fixed is what makes the row length the only variable, so the ratio between the two matrices is the number to carry away, not the level.

Results

Re-verified on a Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88) on 2026-09-02. Correctness, lane-slot arithmetic and the input-dependent winner reproduced. Absolute SpMV timings rose by roughly 1.7x to 2.4x, and the skewed warp-per-row advantage grew from 1.07x to 1.51x.

The original CUDA 12.6 table below and the CUDA 13 transcript both remain in evidence.

GPU: Tesla T4 (compute capability 7.5, 40 SMs)

Part 1: two matrices, 524288 rows, 8388608 nonzeros each
matrix     nonzeros  min row  max row  mean row
-------- ---------- -------- -------- ---------
even        8388608       16       16     16.00
skewed      8388608        3      419     16.00

Lane-slots issued, and the share with a nonzero to do
matrix   kernel             lane-slots   useful
-------- ---------------- ------------ --------
even     thread per row        8388608   100.0%
even     warp per row         16777216    50.0%
skewed   thread per row      219676672     3.8%
skewed   warp per row         23592960    35.6%

Part 2: y = A x, 10 timed runs after 3 warm-ups, 16777216 flops per run
matrix   kernel                   time (ms)     GFLOP/s   x floor
-------- ---------------------- ----------- ----------- ---------
even     stream probe (floor)         0.265        63.4      1.00
even     thread per row               0.303        55.4      1.14
even     warp per row                 0.604        27.8      2.28
skewed   stream probe (floor)         0.285        58.8      1.00
skewed   thread per row               0.849        19.8      2.98
skewed   warp per row                 0.793        21.1      2.78

All six rows above ran the same 16777216 flops over the same 8388608
values and 8388608 column indices. The two matrices hold the same
nonzeros in a different set of rows, and nothing else differs.

The data shape decides the winner, and the lane-slot table shows why before any timing does. Both matrices have 8,388,608 nonzeros and a mean row length of 16. On the even matrix, thread-per-row issues 8,388,608 lane slots and every one has work: 100 percent useful.

On the skewed matrix the same kernel issues 219,676,672 lane slots and 3.8 percent of them do anything. A warp is 32 lanes stepping together, so a warp holding one row of length 419 and 31 rows of length 3 runs for 419 steps with 31 lanes idle for 416 of them. The work did not grow. The waiting did.

Warp-per-row inverts the trade. It is worse on the even matrix, 50 percent useful against 100, because a 32-lane warp on a 16-element row wastes half itself every time. On the skewed matrix it is 35.6 percent useful against 3.8, because the long rows now have a whole warp to chew through them and the short rows finish quickly.

Lane-slot counts do not map directly to time. Under CUDA 12.6, thread-per-row won the even matrix 2.0x and warp-per-row won the skewed matrix by 7 percent. Under CUDA 13.0, the same kernels won by 1.52x (0.696 against 1.059 ms) and 1.51x (1.365 against 2.066 ms).

The lane-slot arithmetic did not change. It predicts a 9.3x advantage for the skewed matrix, much larger than either measured result. Both kernels stay several times above the stream floor because the gather of x also limits them.

The lane-slot table predicts which kernel wins. The gather helps set the margin.

This is the first lesson in the course where the right kernel depends on the input rather than the algorithm. You cannot pick between these two from the source code. You have to look at the row-length distribution, which is why the program prints min, max and mean before it prints a single time.

Run it yourself

Use a CUDA GPU with compute capability 7.5 or later. The README gives the build line:

nvcc -std=c++17 -O3 -arch=sm_75 -o spmv spmv.cu

Compiler Explorer cannot host this one. Building both matrices on the host uses about 235 mebibytes, and two CPU references each walk 8,388,608 nonzeros before the first kernel launch.

That work exceeds the embed's 20-second limit. The kernels use no shared memory, and the device side uses 136 mebibytes.

Exercise

Sort the rows of the skewed matrix by length, rebuild CSR in that order, and run the thread-per-row kernel on the permuted matrix. Predict the new lane-slot count before you build it, then report the measured ratio against the unsorted skewed run and say in one sentence why the kernel got faster without changing.

Time: 30 to 40 minutes. Submit: the predicted lane-slot count, the measured ratio, and the sentence.

Check: the harness un-permutes your output and compares all 524,288 rows against the same CPU reference, so a permutation that drops or duplicates a row fails instead of winning. It prints the first bad row, what it got, what it wanted, and the length of the row that owns it. It then runs the program's own lane-slot arithmetic on your sorted row lengths and prints it beside the measured ratio, pass or fail, because those two numbers are the answer and the pass is only the gate.

Hint 1

You cannot make a warp faster than its longest row. You can choose which rows share a warp.

Hint 2

Count the lane-slots after sorting. Both row lengths occur a multiple of 32 times here, so every group of 32 rows now holds 32 rows of one length. What does that do to the sum of 32 * max(row length) over the groups?

Solution

8,388,608 lane-slots, which is exactly the nonzero count, so the kernel wastes none. Sorting puts the 16,384 heavy rows into 512 warps of their own, each costing 32 * 419. It puts the 507,904 light rows into 15,872 warps, each costing 32 * 3.

Those groups cost 6,864,896 plus 1,523,712 lane slots. The unsorted matrix cost 219,676,672 for the same work, so the sorted layout issues 26 times fewer lane slots without changing the kernel.

The measured ratio will be smaller than 26, for the slots-are-not-seconds reason the Results section measured: both kernels run 2.8 to 3.0 times off the stream floor on the skewed matrix, bound by the gather rather than by issue slots, so removing issue slots cannot produce the full 26x gain.

Sorting costs one reordering of the matrix's rows up front and an un-permute of y after every multiply. You can avoid the un-permute if the rest of the solver keeps the sorted order.

ELL, hybrid, and adaptive CSR formats use the same principle. They rearrange the data so that rows sharing a warp have similar lengths.

Pitfalls

Your two kernels are timed on two different matrices. Two SpMV numbers are comparable only when the nonzero count, the values and the column pattern are identical. Derive the second matrix from the first and assert the nonzero count at compile time, the way the program here does; a kernel that looks twice as fast because it was handed half the work is the bug day 11 opens with.

Your warp-per-row kernel returns garbage for some rows. __shfl_down_sync(0xffffffff, ...) promises that all 32 lanes reach it. If the row index is not warp-uniform, because the block size is not a whole number of warps or the guard is computed per lane, some lanes have already exited and the mask names threads that will never arrive.

Keep the guard on a value that is constant across the warp. Day 23 has the rules.

You decide one of these kernels is the better one. Neither is. The winner is a property of your matrix's row-length distribution, and a kernel that wins on a finite-element mesh loses on a web graph. Print the histogram of row lengths before you print a benchmark.

Your offsets are int and your matrix has more than two billion nonzeros. CSR's row offsets and column indices are 32-bit in almost every published kernel, including the one on this page, because that is the common cuSPARSE path. Past 2^31 nonzeros the offsets wrap silently and the answer is wrong in the middle of the matrix, not at the end. cuSPARSE has a 64-bit index type for this; a hand-written kernel needs the change made by hand.

You profile as a normal user and get nothing. ncu reports ERR_NVGPUCTRPERM, in full: The user running <tool_name/application_name> does not have permission to access NVIDIA GPU Performance Counters or the Hardware Event System on the target device (https://developer.nvidia.com/nvidia-development-tools-solutions-err_nvgpuctrperm-permission-issue-performance-counters , checked 2026-08-30). On a machine where you have root, sudo ncu works, because the stock driver default RmProfilingAdminOnly is 1. Colab and Kaggle give you neither.

Go deeper

Next

Day 38 runs a breadth-first search as a sequence of these products, where the frontier changes shape every level, so the row-length distribution you just measured becomes a property of the iteration rather than of the matrix. Day 40's PageRank is this SpMV in a loop with a convergence reduction on the end, and day 82 puts both of today's kernels against cusparseSpMV on this same skewed matrix.