Compressed Sparse Row Storage and the Cost of Insertion

A sparse matrix from a discretized partial differential equation has perhaps five to seven nonzeros per row, regardless of how fine the mesh becomes. Storing it as a dense two-dimensional array wastes memory in proportion to the square of the grid dimension, and worse, it makes the matrix-vector product read gigabytes of zeros. The compressed sparse row format is the standard response, and its layout is shaped almost entirely by one operation: multiplying the matrix by a vector.

That focus explains both what the format does well and the one thing it does badly. Assembly, the phase in which the matrix is filled in entry by entry, is the operation the format handles worst, which is why every serious sparse library assembles in a different representation and converts once.

Three arrays, no pointers

CSR stores a matrix as three flat arrays. values holds the nonzero entries in row-major order, col_ind holds the column index of each entry, and row_ptr holds n_rows + 1 integers such that the entries of row i occupy positions row_ptr[i] through row_ptr[i+1] - 1. The last entry of row_ptr equals the number of nonzeros, so the length of a row is a subtraction rather than a stored count.

The format is exactly the adjacency structure of a directed graph with edge weights, which is why the same layout appears in graph libraries under a different name. Row i lists its out-neighbors in col_ind, with the weights alongside.

#include <stdlib.h>

/* y = A*x for an n x n matrix in CSR form. */
void csr_matvec(int n, const int *row_ptr, const int *col_ind,
                const double *values, const double *x, double *y)
{
    for (int i = 0; i < n; i++) {
        double sum = 0.0;
        int start = row_ptr[i], end = row_ptr[i + 1];
        for (int k = start; k < end; k++) {
            sum += values[k] * x[col_ind[k]];   /* gather from x */
        }
        y[i] = sum;                             /* scatter-free: y is sequential */
    }
}

The byte cost is easy to state: eight bytes for a double, four for a 32-bit column index, so twelve bytes per nonzero, plus four bytes per row for row_ptr and eight per row for the output vector. A three-dimensional seven-point stencil on a grid with N unknowns has 7N nonzeros, which is 84N bytes of matrix, against 8N^2 bytes for the dense form. The sparse form wins by a factor of N/10, and the factor grows with the problem.

Column indices are usually sorted within a row. The format does not require it, but several algorithms rely on it: merging rows, computing the transpose, and the triangular solves in a direct factorization all become linear passes when the indices are ordered, and degrade to searches when they are not.

The product is limited by memory, not by arithmetic

A matrix-vector product performs two flops per nonzero, one multiply and one add. It moves twelve bytes per nonzero. The ratio is about one flop per six bytes, which on any modern processor is far below the ratio at which the floating-point units become the bottleneck; the kernel is memory bound by roughly an order of magnitude.

The accesses are also irregular. The value and column index for a given nonzero are contiguous, so reading them is a sequential stream, but the corresponding entry of x is at a position determined by the column index, which for an unstructured mesh is effectively random. Each such access may miss in the last-level cache, and the hardware cannot prefetch it because the address is not known until the index has been loaded.

Two structural choices follow. First, the rows must be processed in a way that exposes parallelism, since a single row’s work is a serial chain of dependent additions; splitting the row into chunks and combining partial sums, or processing independent rows concurrently, restores parallelism at the cost of extra additions. Second, the ordering of the unknowns changes the locality of the gather. A bandwidth-reducing permutation of the matrix, such as reverse Cuthill-McKee, makes the column indices of nearby rows overlap, so the entries of x that are needed stay in cache across several rows. The permutation does not change the number of nonzeros or the flop count; it changes how many of the gathers hit.

There is a third effect that is easy to miss: rows of very different lengths create load imbalance when the rows are distributed across threads. A row with a thousand nonzeros and a row with two are not comparable units of work, and a static partition by row index gives some threads far more work than others. Splitting long rows or distributing by nonzero count addresses it.

Assembly in two passes

Filling a CSR matrix entry by entry is expensive because the arrays are sized in advance and every insertion into the middle of a row shifts the tail. Building an n x n matrix with nnz nonzeros one entry at a time therefore costs O(nnz) per insertion in the worst case, and O(nnz^2) overall for an unlucky order.

The standard construction avoids insertion entirely by using the coordinate format as an intermediate. Each nonzero is recorded as a triple of row, column, and value in an append-only array. Once all triples are known, the conversion to CSR is a counting sort with two linear passes: count the entries per row, prefix-sum the counts into row_ptr, then place each triple at the next free position of its row.

/* Convert COO triplets to CSR. Returns 0 on success. */
int coo_to_csr(int n, int nnz, const int *rows, const int *cols,
               const double *vals, int *row_ptr, int *col_ind, double *values)
{
    for (int i = 0; i <= n; i++) row_ptr[i] = 0;
    for (int k = 0; k < nnz; k++) row_ptr[rows[k] + 1]++;
    for (int i = 0; i < n; i++) row_ptr[i + 1] += row_ptr[i];   /* prefix sum */

    int *next = malloc((size_t)n * sizeof *next);               /* cursor per row */
    if (!next) return -1;
    for (int i = 0; i < n; i++) next[i] = row_ptr[i];
    for (int k = 0; k < nnz; k++) {
        int p = next[rows[k]]++;
        col_ind[p] = cols[k];
        values[p] = vals[k];
    }
    free(next);
    return 0;
}

The conversion is O(nnz + n) and needs one temporary array of n integers. Duplicate entries, which arise when the same matrix position is assembled from more than one element in a finite element code, survive the conversion as separate entries; whether that is correct depends on the semantics the consumer expects. Some libraries sum duplicates during conversion, others leave them and let the product accumulate them naturally, and a few reject them.

The format’s structural weakness

Insertion into an existing CSR matrix is the operation the layout makes expensive, and it is not an accident of implementation. The arrays are contiguous and densely packed; making room in row i requires moving every entry after position row_ptr[i], which touches the entire tail of the matrix. Amortized growth strategies, such as leaving gaps at row boundaries, trade memory and a more complex index structure for cheaper insertion, and they slow down the product because the gaps break the sequential stream.

This is why sparse direct solvers do not build their factors incrementally in CSR. The symbolic phase determines the exact nonzero pattern of the factor before any arithmetic happens, so the arrays can be allocated once at their final size and filled in a fixed order. That pattern computation is the subject of Fill-In, Elimination Ordering, and Symbolic Factorization, and it is what makes the format viable for factors that are far denser than the original matrix.

Iterative methods have an easier time because they only need the product, and the matrix is fixed for the duration of the solve. The cost model of those methods is the subject of Krylov Subspace Iteration and Why GMRES Restarts, and it is expressed almost entirely in units of the product that this layout implements. A matrix that changes every iteration, as in a nonlinear solver that reassembles the Jacobian, pays the conversion cost each time. That cost is linear in the number of nonzeros, so it is usually small compared with the arithmetic, but it is not negligible when the product itself is only a few flops per entry.

Formats for regular structure

CSR treats every row independently, which is wasteful when the rows are all the same length and the pattern is known in advance. A five-point stencil on a rectangular grid has exactly five nonzeros per interior row, at offsets that depend only on the row’s position in the grid, and storing the five column indices for every row repeats information the index arithmetic already encodes.

The diagonal format stores a fixed set of diagonals as dense vectors, one per offset, with padding at the boundaries. It needs no column indices at all, so it costs eight bytes per stored entry instead of twelve, and the product becomes a sequence of vector operations over contiguous memory that vectorizes cleanly. The price is paid at the boundaries, where the padding entries are either zero or masked out, and in the loss of generality: a matrix whose pattern does not match the diagonal set cannot be stored this way.

Block formats sit between the two. The block sparse row format stores small dense r x c blocks, with the column index referring to block columns. A matrix arising from a system of coupled equations at each mesh point has a natural block structure, with a small dense block for each pair of neighboring points, and storing it as blocks reduces the index overhead by a factor of r*c and turns the product into small dense kernels. The blocks must be uniform in size, or the format degenerates into something close to CSR with padding.

Integer width and the limits that follow

A 32-bit signed column index caps the number of stored entries at 2^31 - 1, about 2.1 billion. That is not a large number for a three-dimensional problem: a 128^3 grid with seven nonzeros per row already exceeds 14 million, and a factor with substantial fill can reach the cap on a machine with enough memory to hold it. When the cap is approached, the index arrays must widen to 64 bits, which increases the storage per entry from twelve bytes to sixteen and slows the product by roughly the ratio of the bytes moved.

The choice is usually made at build time by a library-wide configuration, so a code linked against a 32-bit-index build cannot handle a matrix that exceeds the limit no matter how much memory is available. Mixed setups, where the row pointers are wide and the column indices narrow, exist but add a conversion step at every interface.

The other limit is the assumption that the matrix is truly sparse. A factorization of a sparse matrix can be much denser than the original, and a matrix that fits comfortably in CSR form before elimination may not fit after it. The memory for the factor is determined by the ordering rather than by the input, and estimating it before committing to the factorization is a separate problem that the symbolic phase exists to solve.