Start New Discussion within our Software Development Community

Gaussian elimination is designed to solve systems of linear algebraic equations,
[A][X] =

This program reads input from a file which is passed as an argument when calling the program. Say the matrix is in file name matrix.txt so we call the program in this way:
-------------------------
gauss matrix.txt
------------------------

In the matrix.txt file the inputs should be in this way:
1. first the dimension of the matrix n, this program works ONLY with square matrices.
2. The coeffiecients of the variables of the equations
3. The constants of the equations, respectively.

For example a sample input file would be like this
------------------------------------------
3

3 -0.1 -0.2
0.1 7 -0.3
0.3 -0.2 10

7.85
-19.3
71.4
------------------------------------------
Entries just need to be white-space separated.

About the algorithm:
This program includes modules for the three primary operations of the Gauss elimination algorithm: forward elimination, back substitution, and partial pivoting.

It implements scaled partial pivoting to avoid division by zero, and during pivoting it also checks if any diagonal entry is zero, thus detecting a singular system.

The equations are not scaled, but scaled values are used to determine the best partial pivoting possible.

Matrices are allocated in such a way so that the index starts from 1, not from zero. So the first element of a matrix A is A[1][1],
not A[0][0].

The co-efficient matrix must be a square one. The constant matrix has been implemented as a column. Although it could have been implemented as an n x 1 matrix, I chose not to.

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

typedef double ** Matrix;
typedef double * Row;
typedef double * Col;
typedef double Elem;

Matrix allocate_matrix(int n);
Col allocate_col(int n);
Row allocate_row(int n);
void free_matrix(Matrix M, int n);

void pivot_partial(Matrix A, Col S,Col B, int n);
void forward_elimination(Matrix A,Col B,int n);
Col back_substitution(Matrix A, Col B, int n);
Col scale_factor(Matrix A,int n);
void gauss(Matrix A, Col B, int n);

void swap_rows(Row *r1, Row*r2);
void print_matrix(Matrix M, int n, char * name);
void print_col(Col C, int n, char *name);
void print_row(Row R, int n, char *name);



int main(int argc, char *argv[])
{
	FILE *ifp;
	int n,i,j;
	Matrix A;
	Col B;
	if(argc < 2)
	{
		printf("\nInput filename not passed \n");
		exit(1);
	}
		
	ifp = fopen(argv[1],"r");
	
	if(ifp == NULL)
	{
		printf("\nCould not open file %s\n",argv[1]);
		exit(1);
	}
	fscanf(ifp,"%i",&n);
	printf("\nDimension = %i\n",n);
	
	A = allocate_matrix(n);
	for( i = 1; i <= n; ++i)
		for(j = 1; j <= n; ++j)
			fscanf(ifp,"%lf", &A[i][j]);
	
	B = allocate_col(n);	
	
	for(j = 1; j <= n; ++j)
		fscanf(ifp,"%lf",&B[j]);
	fclose(ifp);	
	
	
	print_matrix(A,n,"A");
	print_col(B,n,"B");
	
	gauss(A,B,n);	
	
	free_matrix(A,n);
	free(B + 1);
	
	printf("\nDone");
	getchar();
	return 0;
}




void print_matrix(Matrix M, int n, char * name)
{
	int i,j;
	printf("\n[%s] = ",name);
	printf("\n\n");
	for(i = 1; i <= n; i++)
	{
		for(j = 1; j <= n; ++j)
			printf("%6lG ",M[i][j]);
		printf("\n");
	}
}

void print_col(Col C, int n, char * name)
{
	int j;
	printf("\n[%s] = ",name);
	printf("\n\n");
	for(j = 1; j <= n; ++j)
		printf("%6lg\n",C[j]);
			
}
void print_row(Row R, int n, char * name)
{
	int i;
	printf("\n[%s] = ",name);
	for(i = 1; i <= n; ++i)
		printf("%6lg ",R[i]);
	printf("\n");
}

Matrix allocate_matrix(int n)
{
	Matrix A;
	int i,j;
	A = malloc(n * sizeof(Row));
	if(!A)
	{
		printf("\nError : Could not allocate memory for matrix\n");
		exit(1);
	}		
	--A;
		
	for(i = 1; i <= n; ++i)
	{
		A[i] = malloc(n * sizeof(Elem));
		if(!A[i])
		{
			printf("\nError : Could not allocate memory for matrix\n");
			exit(1);
		}
		--A[i];
	}
	return A;
}

void free_matrix(Matrix M, int n)
{
	int i;
	for(i = 1; i <= n; ++i)
			free(M[i] + 1);
	free(M + 1);
}

Col allocate_col(int n)
{
	Col B;
	
	B = malloc(n * sizeof(Elem));
	
	if(!B)
	{
		printf("\nError : could not allocate memory\n");
		exit(1);
	}
	--B;
	return B;
}

Row allocate_row(int n)
{
	Row B;
	B = malloc(n * sizeof(Elem));
	if(!B)
	{
		printf("\nError : could not allocate memory\n");
		exit(1);
	}
	--B;
	return B;
}


Col scale_factor(Matrix A, int n)
{
	int i,j;
	Col S ;
	S = allocate_col(n);
	
	for(i = 1; i <= n; ++i)
	{
		S[i] = A[i][1];
		for(j = 2; j <= n; ++j)
		{
			if(S[i] < fabs(A[i][j]))
				S[i] = fabs(A[i][j]);
		}
		
	}
	return S;
		
}

void pivot_partial(Matrix A, Col S,Col B, int n)
{
	int i,j;
	Elem temp;
	for(j = 1; j <= n; ++j)
	{
		for(i = j + 1; i <= n; ++i)
		{
			if(S[i] == 0)
			{
				if(B[i] == 0)
					printf("\nSystem doesnt have a unique solution");
				else 
					printf("\nSystem is inconsistent");
				exit(1);
			}
			if(fabs(A[i][j]/S[i]) > fabs(A[j][j]/S[j]))
			{
				swap_rows(&A[i],&A[j]);
				temp = B[i];
				B[i] = B[j];
				B[j] = temp;
			}
		}
		
		if(A[j][j] == 0)
		{
			printf("\nSingular System Detected\n");
			exit(1);
		}
	}
	
}

void swap_rows(Row *r1, Row*r2)
{
	Row temp;
	temp = *r1;
	*r1 = *r2;
	*r2 = temp;	
}

void forward_elimination(Matrix A,Col B,int n)
{
	int i,j,k;
	double m;
	
	for(k = 1; k <= n-1; ++k)
	{
		for(i = k + 1; i <= n; ++i)
		{
			m =  A[i][k] / A[k][k];
			for(j = k + 1; j <= n; ++j)
			{
				A[i][j] -= m * A[k][j];
				if(i == j && A[i][j] == 0)
				{
					printf("\nSingular system detected");
					exit(1);
				}					
			}
			B[i] -= m * B[k];
			
		}
	}
			
}

Col back_substitution(Matrix A, Col B, int n)
{
	int i,j;
	Elem sum;
	Col X = allocate_col(n);
	X[n] = B[n]/A[n][n];
	for(i = n - 1; i >= 1; --i)
	{
		sum = 0;
		for(j = i + 1; j <= n; ++j)
			sum += A[i][j] * X[j];
		X[i] = (B[i] - sum) / A[i][i];
	}
	return X;
}

void gauss(Matrix A, Col B, int n)
{
	int i,j;
	Col S, X;
	S = scale_factor(A,n);
	pivot_partial(A,S,B,n);
	forward_elimination(A,B,n);
	X = back_substitution(A,B,n);
	print_col(X,n,"X");
	
	free(S + 1);
	free(X + 1);
}

I have tried to apply the same implementation on Matrix generated based on Finite Field
I made sure that I have changed all the operations of (add, sub, div and mul) to be FF operations.. However I did not get and Identity Matrix in the end ? Can you please help me