Gaussian Elimination on C

Gaussian elimination is one of popular methods of solving of linear systems. Let me present it’s implementation on C.

Input and output data

This implementation accepts input data as one dimension array of doubles with size N x (N + 1), where N is a dimension of the matrix. Rows of extended matrix are placed in this array consequently, and every N-th item of array is a free member. As a result of execution of Gaussian elimination this array will be transformed in the same place to form [1, 0, … , 0, x(0), 0, 1, … 0, x(1), … 0, 0, … 1, x(N-1)], where x(0), x(1), … x(N-1) are roots of this system, if system can be solved. If system can not be solved, algorithm should determine this situation and return gracefully.

gaussian-elimination.h

#pragma once

/* Index of matrix cell with row “r” and column “c” in array with matrix’s size “s” */
#define I(s,r,c) ((r)*((s) + 1) + (c))

/**
* Perform Gaussian elimination
*
* Arguments:
* size – size of matrix
* matrix – matrix of a system
*
* Result: 0 – system can not be solved, 1 – system is solved
*/
const int gaussianElimination(const int size, double *matrix);

gaussian-elimination.c

#include "gaussian-elimination.h"

/* Swap two rows of matrix */
void swap(const int size, double *matrix, int r1, int r2) {
if (r1 == r2) {
return;
}
for (int c = 0; c < size + 1; c++) {
const double t = matrix[I(size, r1, c)];
matrix[I(size, r1, c)] = matrix[I(size, r2, c)];
matrix[I(size, r2, c)] = t;
}
}

/* Divide matrix’s row to non-zero double */
void divide(const int size, double *matrix, int r, double k) {
if (1 == k) {
return;
}
for (int c = 0; c < size + 1; c++) {
matrix[I(size, r, c)] /= k;
}
}

/* Subtract from matrix’s row with index “rTo” matrix’s row with index “rFrom”, multiplied to “k”. If row without free members contains only zeros, system can not be solved. */
int subtractMultiplied(const int size, double *matrix, int rTo, int rFrom, double k) {
if (0 == k) {
return 1;
}
int hasNonZero = 0;
for (int c = 0; c < size + 1; c++) {
matrix[I(size, rTo, c)] -= k * matrix[I(size, rFrom, c)];
if (0 != matrix[I(size, rTo, c)]) {
hasNonZero = 1;
}
}
return hasNonZero;
}

/* Find row with index greater than “r”, when cell in column “r” is not a zero, and swap it with the row with index “r”.
If such row can not be found, it means, that system can not be solved. */
const int findFirstNotZero(const int size, double *matrix, const int r) {
for (int i = r + 1; i < size; i++) {
if (0 != matrix[I(size, i, r)]) {
swap(size, matrix, i, r);
return 1;
}
}
return 0;
}

/* Transform matrix with purpose to make main diagonal does not contain zeros. If this transformation impossible, than system can not be solved */
const int sort(const int size, double *matrix) {
for (int r = 0; r < size; r++) {
if (0 == matrix[I(size, r, r)]) {
if (!findFirstNotZero(size, matrix, r)) {
return 0;
}
}
}
return 1;
}

/* Gaussian elimination */
const int gaussianElimination(const int size, double *matrix) {
/* Preliminary matrix transformation: main diagonal should not contain zeros */
if (!sort(size, matrix)) {
return 0;
}

/* “Move down”: main diagonal should contain only 1 and all cells below it should be 0 */
for (int r = 0; r < size; r++) {
const double k = matrix[I(size, r, r)];
divide(size, matrix, r, k);
for (int i = r + 1; i < size; i++) {
const double n = matrix[I(size, i, r)];
if (!subtractMultiplied(size, matrix, i, r, n)) {
return 0;
}
}
}

/* “Move up”: all cells above main diagonal should be 0. Column of free members should contain roots of this system. */
for( int r=size-1; r>0; r– ) {
for( int i=r-1; i>=0; i– ) {
const double n = matrix[I(size, i, r )];
if( !subtractMultiplied( size, matrix, i, r, n )) {
return 0;
}
}
}
return 1;
}

 How to print all calculated roots?

It’s simple:
...
for( int i=0; iprintf("X(%d)=%f\n", i, matrix[I(size,i,size)]);
}
...

Source https://greywraithblog.wordpress.com/2017/05/28/gaussian-elimination-on-c/