//Modified Code from C Numerical Methods Text

#include <stdio.h>
#include <math.h>
#define MAXSIZE 20

//function prototype
int gauss (double a[][MAXSIZE], double b[], int n, double *det);

int main(void)
{
    double a[MAXSIZE][MAXSIZE], b[MAXSIZE], det;
    int i, j, n, retval;
    
	printf("\n \t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS");
	printf("\n \t USING GAUSSIAN ELIMINATION \n");
	printf("\n This program uses Gaussian Elimination to solve the");
	printf("\n system Ax = B, where A is the matrix of known");
	printf("\n coefficients, B is the vector of known constants");
	printf("\n and x is the column matrix of the unknowns.");

	//get number of equations
	n = 0;
    while(n <= 0 || n > MAXSIZE)
	{
		printf("\n Number of equations: ");
		scanf ("%d", &n);
	}

	//read matrix A
	printf("\n Enter elements of matrix [A]\n");
	for (i = 0; i < n; i++)
		for (j = 0; j < n; j++)
		{
			printf(" A(%d,%d) = ", i + 1, j + 1);
			scanf("%lf", &a[i][j]);
		}
    //read {B} vector
    printf("\n Enter elements of [b] vector\n");
	for (i = 0; i < n; i++)
	{
		printf(" B(%d) = ", i + 1);
		scanf("%lf", &b[i]);
	}
    
	//call Gauss elimination function
	retval = gauss(a, b, n, &det);

	//print results
	if (retval == 0)
	{
		printf("\n\t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS\n");
		printf("\n\t The solution is");
		for (i = 0; i < n; i++)
			printf("\n \t x(%d) = %lf", i + 1, b[i]);
		printf("\n \t Determinant = %lf \n", det);
	}
	else
		printf("\n \t SINGULAR MATRIX \n");

    return 0;
 }

/* Solves the system of equations [A]{x} = {B} using       */
/* the Gaussian elimination method with partial pivoting.  */
/* Parameters:                                             */
/*      n         - number of equations                    */
/*      a[n][n]   - coefficient matrix                     */
/*      b[n]      - right-hand side vector                 */
/*      *det      - determinant of [A]                     */

int gauss (double a[][MAXSIZE], double b[], int n, double *det)
{
    double tol, temp, mult;
	int npivot, i, j, l, k, flag;

	//initialization
	*det = 1.0;
	tol = 1e-30;        //initial tolerance value
	npivot = 0;

	//forward elimination
	for (k = 0; k < n; k++)
	{
		//search for max coefficient in pivot row- a[k][k] pivot element
		for (i = k; i < n; i++)
		{
			if (fabs(a[i][k]) > fabs(a[k][k]))
			{
				//interchange row with maxium element with pivot row
				npivot++;
				for (l = 0; l < n; l++)
				{
				   temp = a[i][l];
				   a[i][l] = a[k][l];
				   a[k][l] = temp;
				}
				temp = b[i];
				b[i] = b[k];
				b[k] = temp;
			}
		}
		//test for singularity
		if (fabs(a[k][k]) < tol)
		{
           //matrix is singular- terminate
		   flag = 1;
		   return flag;
		}
        //compute determinant- the product of the pivot elements
		*det = *det * a[k][k];

		//eliminate the coefficients of X(I)
		for (i = k; i < n; i++)
		{
			mult = a[i][k] / a[k][k];
			b[i] = b[i] - b[k] * mult;  //compute constants
			for (j = k; j < n; j++)     //compute coefficients
				a[i][j] = a[i][j] - a[k][j] * mult;
		}
	}
	//adjust the sign of the determinant
	if(npivot % 2 == 1)
	  *det = *det * (-1.0);

	//backsubstitution
	b[n] = b[n] / a[n][n];
	for(i = n - 1; i > 1; i--)
	{
		for(j = n; j > i + 1; j--)
			b[i] = b[i] - a[i][j] * b[j];
		b[i] = b[i] / a[i - 1][i];
	}
	flag = 0;
	return flag;
}

//eliminate the coefficients of X(I)
		for (i = k + 1; i < n; i++)
		{
			mult = a[i][k] / a[k][k];
			b[i] = b[i] - b[k] * mult;  //compute constants
			for (j = k; j < n; j++)     //compute coefficients
				a[i][j] = a[i][j] - a[k][j] * mult;
		}
	}