DZone Snippets is a public source code repository. Easily build up your personal collection of code snippets, categorize them with tags / keywords, and share them with the world

Adrian has posted 1 posts at DZone. View Full User Profile

Matrix operations: Gaussian elimination, Upper triangulization and Back-substitution example in C

11.03.2012
| 6659 views |
  • submit to reddit
/** Gaussian elimination / Upper triangulization / Back-substitution example by Adrian Boeing
 adrianboeing.blogspot.com
 www.adrianboeing.com
 */

#include <stdio.h>

void PrintMatrix(double **mat, int m, int n) {
    for (int j=0;j<n;j++) {
        for (int i=0;i<m;i++) {
            printf("%+5.2f ",mat[j][i]);
        }
        printf("\n");
    }
}

int main(int argc, char*argv[]) {
    int m,n; //m - col, n - rows
    int i;
    m=5;
    n=4;
    
    //allocate the matrix
    double **mat = new double* [n];
    for (i=0;i<n;i++)
        mat[i] = new double [m];
    
    //initialize matrix
    mat[0][0] = 1; mat[0][1] = 2; mat[0][2] = 1; mat[0][3] = 4; mat[0][4] = 13;
    mat[1][0] = 2; mat[1][1] = 0; mat[1][2] = 4; mat[1][3] = 3; mat[1][4] = 28;
    mat[2][0] = 4; mat[2][1] = 2; mat[2][2] = 2; mat[2][3] = 1; mat[2][4] = 20;
    mat[3][0] =-3; mat[3][1] = 1; mat[3][2] = 3; mat[3][3] = 2; mat[3][4] = 6;
    
    printf("Initial matrix\n");
    PrintMatrix(mat,m,n);

    printf("Gaussian elimination\n");
    //p is the pivot - which row we will use to eliminate
    for (int p=0;p<n-1;p++) { //pivot through all the rows
        for (int r=p+1; r < n; r++) { //for each row that isn't the pivot
            float multiple = mat[r][p] / mat[p][p]; //how many multiples of the pivot row do we need (to eliminate this row)?
            for (int c = 0; c<m; c++) { //for each element in this row
                mat[r][c] = mat[r][c] - mat[p][c]*multiple; //subtract the pivot row element (multiple times)
            }
        }
    }
    PrintMatrix(mat,m,n);
    
    printf("Upper triangulization\n");
    //p is the pivot - which row we will normalize
    for (int p=0;p<n;p++) { //pivot through all the rows
        for (int c=p+1;c<m;c++) { //normalize the pivot row, so that the pivot element can be set to 1.
            mat[p][c] /= mat[p][p]; //divide all elements in the row by the pivot element
        }
        mat[p][p] = 1; //set the pivot element to 1.
    }
    PrintMatrix(mat,m,n);
    
    printf("Back-substitution\n");
    for (int p=n-1;p>0;p--) { //pivot backwards through all the rows
        for (int r=p-1;r>=0;r--) { //for each row above the pivot
            float multiple = mat[r][p]; //how many multiples of the pivot row do we need (to subtract)?
            for (int c=p-1;c<m;c++) {
                mat[r][c] = mat[r][c] - mat[p][c]*multiple; //subtract the pivot row element (multiple times)
            }
        }
    }
    PrintMatrix(mat,m,n);
}

/*
 Another example:
 m = 6;
 n = 3;
 
 mat[0][0] = 2; mat[0][1] = 4; mat[0][2] =-2; mat[0][3] = 1; mat[0][4] = 0; mat[0][5] = 0;
 mat[1][0] = 4; mat[1][1] = 9; mat[1][2] =-3; mat[1][3] = 0; mat[1][4] = 1; mat[1][5] = 0;
 mat[2][0] =-2; mat[2][1] =-3; mat[2][2] = 7; mat[2][3] = 0; mat[2][4] = 0; mat[2][5] = 1;
 */

A simple example of matrix operations (Gaussian elimination, Upper triangulization and Back-substitution) which can be used for solving systems of linear equations, or matrix inversion.