Piece of Algebra

The most useful C-programming snippets on the Internet
to help you find your X

There are certain moments in programming when you stumble on some problem that has definitely been solved previously. You don't need to reinvent the wheel! There are around lots of libraries, like CBLAS and LAPACK that offer a great deal of well-tested and well-documented functions for all sorts of algebraic operations. There are however cases when you do need to re-implement matrix-vector multiplication and other trivial operations. Such a case is when you're writing code for Arduino or some other low-level system. In this post I'm disclosing some of the finest code snippets known to humanity - feel free to copy and paste them in your code, change the variable names and pretend it's your own exploit and impress your friends with your skills.  :~P

Matrices as vectors in C:

Matrices in C are most times represented as arrays (i.e., vector-like entities) - after all, RAM is a linear thing.  It is pretty much straightforward to pack the values of a matrix column-wise or row-wise and map any matrix to a vector. Knowing the dimensions of the initial matrix, it is then easy to access its elements and perform various algebraic operations as we'll see in what follows. The same holds for 3-dimensional and higher dimensional matrices. In what follows, we assume that the matrices are mapped onto arrays column-wise.

Note 1: Here I'm making use of certain types such as size_t and ptrdiff_t.  You have to make sure that there are defined on your system. If not, consider replacing them by unsigned long or unsigned short or unsigned int.

Note 2: Here I'm using the notation double a[] because I work with fixed-size arrays. However, exactly the same code will work for dynamically allocated arrays of type double*.

Matrix-Vector Multiplication:

This functions performs the multiplication alpha*a*x+b and updates (adds on) the array y, i.e., it performs the update y += alpha*a*x+b. If b is NULL, it will not be taken into account. It is assumed that all arrays are of compatible dimensions.

void zmvmult(

      const double a[], 

      const double x[], 

      const double alpha, 

      const double b[], 

      const size_t nRows, 

      const size_t nCols, 

      double y[]

     ) 

{
    ptrdiff_t i, j;
    for (i = 0; i < nRows; i++) {
        for (j = 0; j < nCols; j++) 

           y[i] += alpha * a[i + j * nRows] * x[j];
        if (b != NULL) y[i] += b[i];
    }
}

Multiply the transpose of a matrix with a vector:

void zmvmult_trans (

      const double a[], 

      const double x[], 

      const double alpha, 

      const double b[],

      const size_t nRows, 

      const size_t nCols, 

      double y[]

     ) 

{
    ptrdiff_t i, j;
    for (j = 0; j < nCols; j++) {
        for (i = 0; i < nRows; i++) 

          y[j] += alpha * a[i + j * nRows] * x[i];
        if (b != 0) y[j] += b[j];
    }
}

Solve a linear system whose matrix is lower-triangular:

Consider a system in the form Ly=x, where L is a lower triangular matrix. In this case, the array low_triang contains only the lower-triangular elements of the matrix to avoid redundancy (we don't need to store all those zeros). We assume that the elements of L are packed row-wise, i.e. low_triang = {L(0,0), L(1,0), L(1,1), L(2,0), ...}. It should be simple to modify the code for column-wise packed arrays.

void zlowslv (

      const double low_triang[], 

      double x[], 

      const size_t n

     ) 

{
    ptrdiff_t i, j, offset;
    double sum = 0.0;
    x[0] = x[0] / low_triang[0];
    for (i = 1; i < n; i++) {
        offset = ((i + 1) * i) / 2;
        sum = x[i];
        for (j = 0; j < i; j++) {
            sum -= (x[j] * low_triang[offset + j]);
        }
        x[i] = sum / low_triang[offset + i];
    }
}

this function solves the system Ly=x and updates the array x with the solution of the system using forward substitution. Notice also that the variable x is not declared to be constant.

Solve a linear system whose matrix is upper-triangular:

Let us now see how we can write a function to solve a system of the form L'y=x. where L is a lower-triangular matrix and L' is its transpose. The following function solves this system using backward substitution and updates x with the solution.

void zlowtrslv(

      const double low_triang[], 

      double x[], 

      const size_t n

     ) 

{
    ptrdiff_t i, j, s;
    const size_t size_low = (n * (n + 1)) / 2;
    double sum = 0.0;
    x[n - 1] = x[n - 1] / low_triang[size_low - 1];
    for (i = 1; i < n; i++) {
        sum = x[n - i - 1];
        s = size_low - i - 1;
        for (j = 0; j < i; j++) {
            sum -= (low_triang[s] * x[n - 1 - j]);
            s -= (n - j - 1);
        }
        x[n - i - 1] = sum / low_triang[s];
    }
}