testing gaussian elimination to solve system of equations.

gaussian.c

#include <stdio.h>
#include <math.h>

void gaussian(size_t const n, double A[n][n], double v[const restrict]) {
	for( size_t k = 0; k < n-1; k++ ) {
		/// Partial pivot
		double cur_max = fabs(A[k][k]);
		size_t m = k;
		for( size_t i = k+1; i < n; i++ ) {
			double const potential_max = fabs(A[i][k]);
			if( cur_max < potential_max ) {
				cur_max = potential_max;
				m = i;
			}
		}
		if( m != k ) {
			double const b_temp = v[k];
			v[k] = v[m];
			v[m] = b_temp;
			for( size_t j = k; j < n; j++ ) {
				double const A_temp = A[k][j];
				A[k][j] = A[m][j];
				A[m][j] = A_temp;
			}
		}
		/// Forward elimination
		for( size_t i = k+1; i < n; i++ ) {
			double const factor = A[i][k] / A[k][k];
			for( size_t j = k+1; j < n; j++ ) {
				A[i][j] = (A[i][j] - (factor * A[k][j]));
			}
			v[i] = (v[i] - (factor * v[k]));
		}
	}
	/// Back substitution
	for( size_t i = n-1; i < n; i-- ) {
		v[i] = v[i];
		for( size_t j = i+1; j < n; j++ ) {
			v[i] = (v[i] - (A[i][j] * v[j]));
		}
		v[i] = (v[i] / A[i][i]);
	}
}

int main() {
	enum {
		M_SIZE = 3,
	};
	/// [1  4 6, 2 3 4, 5 2 1]
	double matrix[M_SIZE][M_SIZE] = {
	  {-6,1,1},
	  {1,-2,2},
	  {1,1,-1},
	};
	double v[M_SIZE] = {11, 4, 1};
	
	printf("Matrix :: \n");
	for( size_t row=0; row < M_SIZE; row++ ) {
		for( size_t col=0; col < M_SIZE; col++ ) {
			printf("| %+f |", matrix[row][col]);
		}
		printf("| %+f\n", v[row]);
	}
	
	gaussian(M_SIZE, matrix, v);
	printf("\nSolution :: \n[x]: %+f\n[y]: %+f\n[z]: %+f\n", v[0], v[1], v[2]);
}