/** 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;
*/