Gaussian Elimination C Code

The internet is a wild place full of college homework. It appears that there is a lot of incorrect C code online for Gaussian Elimination with off by one indexing errors. This is an attempt to remedy that.

#include <stdio.h>

#define n 3
#define my_fabs(a) (((a) < 0) ? -(a) : (a))

void gaussian(float A[n][n], float v[n]) {
    float *b = v;
    float *x = v;
    for (int k = 0; k < n-1; k++) {
        // Partial pivot
        float cur_max = my_fabs(A[k][k]);
        int m = k;
        for (int i = k+1; i < n; i++){
            float potential_max = my_fabs(A[i][k]);
            if (potential_max > cur_max) {
                cur_max = potential_max;
                m = i;
            }
        }
        if (m != k) {
            float b_temp = b[k];
            b[k]  = b[m];
            b[m]  = b_temp;
            for (int j = k; j < n; j++) {
                float A_temp = A[k][j];
                A[k][j] = A[m][j];
                A[m][j] = A_temp;
            }
        }
        // Forward elimination
        for (int i = k+1; i < n; i++) {
            float factor = A[i][k]/A[k][k];
            for (int j = k+1; j < n; j++) {
                A[i][j] -= factor*A[k][j];
            }
            b[i] -= factor*b[k];
        }
    }
    // Back substitution
    for (int i = n-1; i >= 0; i--) {
        x[i] = b[i];
        for (int j = i+1; j < n; j++) {
            x[i] -= A[i][j]*x[j];
        }
        x[i] /= A[i][i];
    }
}

int main() {
    float A[n][n] = {{ 0,  1,  1},
                     { 1, -2,  2},
                     { 1,  1, -1}};
    float b[n] = {11, 4, 1};
    float *v = b;
    gaussian(A, v);
    float *x = v;
    printf("The solution is:\n");
    for (int i = 0; i < n; i++) {
        printf("x[%d] = %f\n", i, x[i]);
    }
    return 0;
}

andrew

See author's posts