Sparse matrices as linked lists

Four row heads, each with a chain of nodes sorted by column. Watch seven updates splice nodes in, change one, and unlink one, with nothing shifted.

⋯6struct entry {7    int col;8    double val;9    struct entry *next;10};11 12static struct entry *rows[N];13 14void set(int i, int j, double x)15{16    struct entry **link = &rows[i];17    while (*link != NULL && (*link)->col < j)18        link = &(*link)->next;19    if (*link != NULL && (*link)->col == j) {20        if (x == 0.0) {21            struct entry *dead = *link;22            *link = dead->next;23            free(dead);24        } else {25            (*link)->val = x;26        }27        return;28    }29    if (x == 0.0)30        return;31    struct entry *e = malloc(sizeof *e);32    if (e == NULL)33        exit(1);34    e->col = j;35    e->val = x;36    e->next = *link;37    *link = e;38}⋯48int main(void)49{50    set(0, 2, 5);51    set(0, 0, 1);52    set(1, 3, 2);53    set(0, 1, 4);54    set(0, 2, 0);55    set(1, 3, 6);56    set(3, 0, 7);57    for (int i = 0; i < N; i++) {

A chain per row

CSR packs the entries tightly, but a new entry in the middle means shifting every entry after it. When the pattern of non-zeros keeps changing, chains do better. Each row gets a head, a pointer to the first node of a linked chain, and each node holds one entry:

struct entry {
    int col;
    double val;
    struct entry *next;
};

static struct entry *rows[N];

A node is 24 bytes, and malloc gives it a 32-byte block. The heads are statics, so they start as NULL: an empty matrix. Each chain is kept sorted by column. A search can then stop at the first column larger than the one it wants, and the chain prints in column order.

The code and its cost

set: three cases

set(i, j, x) walks row i’s chain with a pointer to a link, as deleting from a list does: link points at the head or at the next field of a node, so changing *link changes whichever pointer leads to the next node.

struct entry **link = &rows[i];
while (*link != NULL && (*link)->col < j)
    link = &(*link)->next;

The walk stops at the first node whose column is j or more, or at the end. Then there are three cases.

The entry exists, and *link is the column-j node. If x is 0, the node is unlinked and freed. Setting an entry to zero does not store a zero; it removes the node:

if (x == 0.0) {
    struct entry *dead = *link;
    *link = dead->next;
    free(dead);
} else {
    (*link)->val = x;
}

Otherwise the value changes in place.

The entry is missing, and x is 0. The matrix already has 0 there, so nothing happens.

The entry is missing, and x is not 0. A new node is spliced in where the walk stopped, before the first larger column:

struct entry *e = malloc(sizeof *e);
e->col = j;
e->val = x;
e->next = *link;
*link = e;

Because link may point at the head or at a node’s next, the same two lines insert a new first entry and a middle one. No other node moves.

The run

The program makes seven calls on an empty 4 × 4 matrix:

The last node gets the block that set(0, 2, 0) gave back: on this platform, glibc’s malloc hands out a recently freed block of the same size first. C does not promise this; another allocator may return any free block. The program prints the chains and then the matrix through get:

row 0: 0:1 1:4
row 1: 3:6
row 2:
row 3: 0:7
1 4 0 0
0 0 0 6
0 0 0 0
7 0 0 0

At the end, main frees each row’s nodes and then sets rows[i] = NULL. Without that line the heads of rows 0, 1 and 3 would still hold the addresses of freed blocks: dangling pointers, which any later get or set would follow into memory the program no longer owns.

get

get walks the chain and stops as soon as the columns pass j:

for (const struct entry *e = rows[i]; e != NULL && e->col <= j; e = e->next)
    if (e->col == j)
        return e->val;
return 0.0;

Against CSR

Chains never shift: an insertion or a removal changes one or two pointers once the place is found. But finding the place means following pointers from node to node, which is slower than reading neighbouring cells of an array, and each entry takes a 32-byte block, where CSR takes 12 bytes (a double and an int). Reaching an entry costs a walk along its row, as long as the row’s chain.

A fuller version gives each node a second pointer and links the columns too, with an array of column heads, so that a column can be walked as easily as a row.

Running it yourself

gcc -std=c11 -Wall -Wextra -pedantic chains.c -o chains
./chains

Worked example

Appending. On the final state, set(0, 3, 9) walks past columns 0 and 1 and reaches the end of row 0’s chain, so its node goes after the column-1 node.

Memory. Ten non-zeros in an 8 × 8 matrix take 10 blocks of 32 bytes, 320 bytes, plus 8 heads of 8 bytes, 64 bytes: 384 bytes. CSR takes 29 numbers, 10 doubles and 19 ints, 156 bytes.

Nothing to do. set(2, 1, 0) on the empty row 2 finds no entry and x is 0, so it returns without allocating.

Two checks. Why is each chain sorted? (A search can stop at the first larger column, and printing goes in column order.) Which block does (3, 0) get? (0x1000, freed when (0, 2) was set to zero, under glibc’s malloc.)

Where it’s used

Elimination. Solving a sparse system by elimination creates new non-zeros, called fill-in, as it goes; chains take them without shifting.

Assembly. Simulations often build a sparse matrix entry by entry, in no particular order, in a linked form, and convert it to CSR once it is complete, for fast products.

Reference card

Next, polynomials: storing a polynomial keeps its terms the same two ways, in an array or in a chain.

Platform: C11, gcc on x86-64 Linux: int 4 bytes, double 8, pointers 8, little-endian. Addresses in the pictures are short stand-ins, by one rule for the whole course; sizes are the real ones.