Tridiagonal and banded matrices

A 5 × 5 matrix whose numbers sit on three central diagonals. Watch each diagonal move out into a strip of its own.

⋯5static double lo[N - 1], di[N], up[N - 1];6 7double get(int i, int j)8{9    if (i == j)10        return di[i];11    if (i == j + 1)12        return lo[j];13    if (j == i + 1)14        return up[i];15    return 0.0;16}⋯34int main(void)35{36    for (int i = 0; i < N; i++)37        di[i] = 2;38    for (int i = 0; i < N - 1; i++) {39        lo[i] = -1;40        up[i] = -1;41    }42    for (int i = 0; i < N; i++) {43        for (int j = 0; j < N; j++)44            printf("%2g ", get(i, j));45        printf("\n");46    }47    printf("stored %d of %d cells\n", 3 * N - 2, N * N);48    const double b[N] = { 0, 0, 0, 0, 6 };49    double x[N];50    solve(b, x);51    printf("x =");52    for (int i = 0; i < N; i++)53        printf(" %.6g", x[i]);54    printf("\n");55    return 0;56}

Three strips

A tridiagonal matrix is zero except on its main diagonal and the two diagonals beside it: A(i, j) = 0 whenever |i − j| > 1. The main diagonal has n entries and each side diagonal has n − 1, so the matrix holds at most

n + 2(n − 1) = 3n − 2

non-zero entries. For n = 5 that is 13 of 25. The program keeps each diagonal in an array of its own:

static double lo[N - 1], di[N], up[N - 1];

di[i] is A(i, i), lo[j] is A(j + 1, j), just below the diagonal, and up[i] is A(i, i + 1), just above it. The program fills di with 2 and lo and up with −1.

The code and its cost

Reading an entry

get asks which diagonal, if any, holds (i, j):

double get(int i, int j)
{
    if (i == j)
        return di[i];
    if (i == j + 1)
        return lo[j];
    if (j == i + 1)
        return up[i];
    return 0.0;
}

get(3, 2) is lo[2] = −1; get(0, 3) is on none of the three diagonals, so it is 0 with no memory read. The program prints all 25 entries through get:

 2 -1  0  0  0
-1  2 -1  0  0
 0 -1  2 -1  0
 0  0 -1  2 -1
 0  0  0 -1  2
stored 13 of 25 cells

Bands

A banded matrix has its non-zeros within w diagonals of the main one: A(i, j) = 0 whenever |i − j| > w. The number w is the half-width; a tridiagonal matrix has w = 1, a pentadiagonal one w = 2. The band holds 2w + 1 diagonals: the main one with n entries, and two of each length n − 1, n − 2, …, n − w. That makes

n + 2((n − 1) + (n − 2) + … + (n − w)) = (2w + 1)n − w(w + 1)

cells. Store the band as a (2w + 1) × n array, one row per diagonal, and A(i, j) is in row w + j − i and column j: each column of A keeps its column in the array. A few cells at the ends of the side rows go unused. (The program’s three strips are a variant: up[i] holds A(i, i + 1), indexed by row, so its cells sit one place to the left of column j.)

Solving with the strips

The point of the strips is not only the memory. A system Ax = b with a tridiagonal A is solved without ever touching the zeros. The program solves it for b = (0, 0, 0, 0, 6), which is A times (1, 2, 3, 4, 5):

void solve(const double *b, double *x)
{
    double c[N], d[N];
    c[0] = up[0] / di[0];
    d[0] = b[0] / di[0];
    for (int i = 1; i < N; i++) {
        double m = di[i] - lo[i - 1] * c[i - 1];
        if (i < N - 1)
            c[i] = up[i] / m;
        d[i] = (b[i] - lo[i - 1] * d[i - 1]) / m;
    }
    x[N - 1] = d[N - 1];
    for (int i = N - 2; i >= 0; i--)
        x[i] = d[i] - c[i] * x[i + 1];
}

This is the tridiagonal algorithm, also called the Thomas algorithm after the physicist who described it. The sweep down removes the lower diagonal one row at a time. Row i subtracts lo[i − 1] times the row above, which leaves m on the diagonal, and divides by m, so that each row reads x[i] + c[i] x[i + 1] = d[i]. The pass up starts at the last row, which now reads x[4] = d[4], and substitutes each answer into the row above. Replay it here:

⋯18void solve(const double *b, double *x)19{20    double c[N], d[N];21    c[0] = up[0] / di[0];22    d[0] = b[0] / di[0];23    for (int i = 1; i < N; i++) {24        double m = di[i] - lo[i - 1] * c[i - 1];25        if (i < N - 1)26            c[i] = up[i] / m;27        d[i] = (b[i] - lo[i - 1] * d[i - 1]) / m;28    }29    x[N - 1] = d[N - 1];30    for (int i = N - 2; i >= 0; i--)31        x[i] = d[i] - c[i] * x[i + 1];32}⋯48    const double b[N] = { 0, 0, 0, 0, 6 };49    double x[N];50    solve(b, x);51    printf("x =");52    for (int i = 0; i < N; i++)53        printf(" %.6g", x[i]);54    printf("\n");55    return 0;56}

 2 -1  0  0  0 
-1  2 -1  0  0 
 0 -1  2 -1  0 
 0  0 -1  2 -1 
 0  0  0 -1  2 
stored 13 of 25 cells

The values on the diagonal after elimination, m, are 2, 1.5, 1.333…, 1.25 and 1.2. c is −0.5, −0.667, −0.75 and −0.8, and d is 0, 0, 0, 0 and 5. The pass up gives x[4] = 5, then 4, 3, 2 and 1, and the program prints

x = 1 2 3 4 5

In floating point x[0] comes out as 0.9999999999999998, which %.6g prints as 1. The solver needs no row swaps here because each diagonal entry, 2, is at least as large as its two neighbours together, and larger in the first and last rows, so m never gets close to 0.

The cost

Storage is 3n − 2 cells instead of n². The sweep does 6 arithmetic operations a row (m takes two, c one, d three) and the pass up 2, so the solver takes about 8n operations. Counted the same way, each multiplication and each subtraction one operation, Gaussian elimination on a full n × n matrix takes about 2n³/3. Solving a linear system does not always cost about n³: with a band, it costs about n.

Running it yourself

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

Worked example

A large system. For n = 1,000 the three strips hold 3 × 1,000 − 2 = 2,998 cells against 1,000,000 for the full grid. The solver takes about 8,000 operations; elimination on the full matrix takes about 2 × 1,000³ / 3, roughly 667,000,000, counted the same way.

A band of half-width 2. For n = 5 and w = 2 the band holds (2 × 2 + 1) × 5 − 2 × 3 = 25 − 6 = 19 cells: the main diagonal’s 5, two diagonals of 4 and two of 3.

Two reads. get(3, 2) is lo[2] = −1, since row 3 is one below column 2. get(0, 3) is 0: column 3 is three away from row 0.

Two checks. Which array holds the entry in row 2, column 1? (lo, cell 1.) How many cells does a pentadiagonal (w = 2) 100 × 100 matrix need? (5 × 100 − 6 = 494.)

Where it’s used

Heat and vibration along a rod. Split a rod into points; each point’s temperature or displacement depends only on its two neighbours, so the equations form a tridiagonal system.

Cubic splines. Fitting a smooth curve through points with cubic pieces gives one equation per point, linking it to its neighbours: a tridiagonal system.

Implicit time steps. Simulations that step a diffusion equation implicitly solve a tridiagonal or banded system at every step, which takes about n operations each time. The same elimination, written for any matrix, is LU decomposition.

Reference card

Back to diagonal matrices, the band with w = 0. Next, Toeplitz matrices need only one number per diagonal.

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.