/****************************************************************************
* Program 3.1: Gauss Elimination Method for systems of linear equations.	*
*				Coded by Dilip Datta										*
*																			*
*****************************************************************************/

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

/********************************************************************/

int main()
{
	int i, j, n;
	double **a, *b, *x;

	void gauss_elimination(int n, double **a, double *b, double *x);

	printf("\nNumber of equations: ");										/* Step-1 */
	scanf("%d", &n);

	x = (double *)malloc( (n+1)*sizeof(double) );
	b = (double *)malloc( (n+1)*sizeof(double) );
	a = (double **)malloc( (n+1)*sizeof(double *) );
	for(i = 1; i <= n; i++)
		a[i] = (double *)malloc( (n+1)*sizeof(double) );

	for(i = 1; i <= n; i++)
	{
		for(j = 1; j <= n; j++)
		{
			printf("a[%d][%d] = ", i, j);
			scanf("%lf", &a[i][j]);
		}
		printf("b[%d] = ", i);
		scanf("%lf", &b[i]);
	}

	gauss_elimination(n, a, b, x);

	printf("\n\n");
	printf("Solution of the given system of linear equations\n");
	printf("------------------------------------------------\n");
	printf("x = (");
	for(i = 1; i <= n-1; i++)	printf("%lf, ", x[i]);
	printf("%lf)\n\n", x[n]);

	return(0);
}
/********************************************************************/

void gauss_elimination(int n, double **a, double *b, double *x)
/* Sub-routine for Guass elimination method for systems of linear equation */
{
	int i, j, k;
	int p;
	double factor, sum;
	double big, dummy;

	// Forward elimination with partial pivoting
	for(k = 1; k <= n-1; k++)
	{
		// Partial (column) pivoting
		p = k;
		big = fabs(a[k][k]);
		for(i = k+1; i <= n; i++)
		{
			if(big < fabs(a[i][k]))
			{
				big = fabs(a[i][k]);
				p = i;
			}
		}

		if(p != k)
		{
			for(j = k; j <= n; j++)
			{
				dummy = a[p][j];
				a[p][j] = a[k][j];
				a[k][j] = dummy;
			}

			dummy = b[p];
			b[p] = b[k];
			b[k] = dummy;
		}

		// Forward elimination
		for(i = k+1; i <= n; i++)
		{
			factor = a[i][k]/a[k][k];
			for(j = k+1; j <= n; j++)	a[i][j] -= a[k][j]*factor;

			b[i] -= b[k]*factor; 
		}
	}

	// Back substitution
	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;
}
/********************************************************************/