/*
      	Title: This is a demo program for subroutine linsys.
 
      	It will solve a linear system a*x=b, given by the user.
		This program is from Section 8.2.
*/
 
#include <stdio.h>
#include <math.h>

#define  MAX_N	20


void linsys(float mat[][MAX_N+1], float b[MAX_N+1], int pivot[MAX_N+1], int n, int *ier);

main()
{

	float a[MAX_N+1][MAX_N+1], b[MAX_N+1];
	int pivot[MAX_N+1];
	int n, ier, i, j;


	while (1)  
	{ 		

		/*  Input order of linear system  */
		
		printf("\n\n Give the order of the linear system");
		printf("\n Give zero to stop : ");
		scanf("%d", &n);
		if (n==0)
			return(0);

		/*  Input linear system  */
      				
		printf("\n Give the coefficients of the linear system");
		printf("\n One equation at a time ");
		printf("\n Conclude each equation with its right-hand constant");
		
		for (i=1; i <= n; i++)  
		{
			printf("\n Give coefficients of equation : ");
			for(j=1; j <= n; j++)
				scanf("%f", &a[i][j]);
			scanf("%f", &b[i]);
		}

		/* Solve the linear system  */
      		
		linsys(a,b,pivot,n,&ier);

		/* Print the results  */
      		
		printf("\n\n\n n = %2d     ier = %3d\n\n  I        Solution\n",n, ier);
		for (i=1; i <= n; i++)
			printf("\n %3d\t %20.10e", i, b[i]);
	}
}

	

void linsys(float mat[][MAX_N+1], float b[MAX_N+1], int pivot[MAX_N+1], int n, int *ier)
{

/*	
	This routine solves a system of linear equations
        	a*x = b
	The method used is gaussian elimination with
	partial pivoting.
	
	Input:
	The coefficient matrix  a  is stored in the array mat.
	The right side constants are in the array  b.
	The order of the linear system is  n.
	The variable  md  is the number of rows that  mat
	is dimensioned as having in the calling program.
	The size of  maxpiv, given below, must be greater
	than  n.  If not, it is a fatal error.
	
	Output:
	The array  b  contains the solution  x.
	mat contains the upper triangular matrix  u
	obtained by elimination.  The row multipliers
	used in the elimination are stored in the
	lower triangular part of  mat.
	ier=0 means the matrix  a  was computationally
	nonsingular, and the gaussian elimination
	was completed satisfactorily.
	ier=1 means that the matrix  a  was
	computationally singular.
*/

int i, j, k;
float mult, temp, sum, absa, amax;


	/*  Begin elimination steps  */

	for (k=1; k <= n-1; k++)  
	{ 
		
		/*  Choose pivot row  */
 		pivot[k] = k;
 		amax = fabs(mat[k][k]);
  		for (i=k+1; i <= n; i++)  
		{
			absa = fabs(mat[i][k]);
			if (absa > amax)  {
				pivot[k] = i;
				amax = absa;
			}
		}
        
		if (amax == 0.0)  
		{ 	
			/*  Coefficient matrix is singular  */            		
			*ier = 1;
			return;
		}
        
		if (pivot[k] != k)  
		{ 	
			/*  Switch rows k and pivot[k]  */
			i = pivot[k];
			temp = b[k];
			b[k] = b[i];
			b[i] = temp;
			for (j=k; j <= n; j++)  
			{
				temp = mat[k][j];
				mat[k][j] = mat[i][j];
				mat[i][j] = temp;
			}
		}

		/*  Perform step #k of elimination */
		for (i=k+1; i <=n; i++)  
		{
			mult = mat[i][k]/mat[k][k];
			mat[i][k] = mult;
			b[i] = b[i] - mult*b[k];
			for (j=k+1; j <= n; j++)
				mat[i][j] = mat[i][j] - mult*mat[k][j];
		}
	}
      
	if (mat[n][n] == 0.0)  { 	/* Coefficient matrix is singular  */
		*ier = 1;
		return;
	}


	/*  Solve for solution x using back substitution  */
	for (i=n; i >= 1; i--)  
	{
		sum = 0.0;
		for (j=i+1; j <= n; j++)
			sum = sum + mat[i][j]*b[j];
		b[i] = (b[i] - sum)/mat[i][i];
	}
      	
	*ier = 0;
	return;
}