/**********************************************************************/
/* Function to perform Gaussian elimination with partial pivoting     */
/* Translated from Fortran code, pages 29-30 of "Numerical Analysis"  */
/* Johnson and Riess, Addison-Wesley Publishing Company		      */
/**********************************************************************/
gauss(A,B,X,n,ierror,rnorm)
float A[NO_POINTS+4][NO_POINTS+4],B[NO_POINTS+4],X[NO_POINTS+4],*rnorm;
int *ierror,n;
{
float AUG[NO_POINTS+4][NO_POINTS+3],pivot,temp,Q,residual,Rmag,RSQ;
int nm1,np1,i,j,k,ipivot,ip1;
char s[50];

nm1 = n-1;
np1 = n+1;

/* Set up the augmented matrix */
for(i=0;i<n;i++)
	{
	for(j=0;j<n;j++)
		AUG[i][j] = A[i][j];
	AUG[i][n] = B[i];
	}
if(DEBUG) print_matrix(AUG,n,np1,"Before Triangulation");

/* The outer loop uses elementary row operations to transform the
   augmented matrix to echelon form, i.e. triangular form.	*/

for(i=0;i<nm1;i++)
	{
	/* Search for the largest entry in column i, rows i to n.
	   ipivot is the row index of the largest entry.   */
	pivot = 0.0;
	for(j=i;j<n;j++)
		{
		temp = fabs(AUG[j][i]);
		if(pivot < temp)
			{
			pivot = temp;
			ipivot = j;
			}
		}
	if(fabs(pivot) < VERY_SMALL)
		{ 
		(*ierror) = 2;
		printf("Error: input matrix is singular\n");
		return;
		}
	/* Interchange row i and row ipivot, if necessary */
	if(ipivot != i)
	for(k=i;k<np1;k++)
		{
		temp = AUG[i][k];
		AUG[i][k] = AUG[ipivot][k];
		AUG[ipivot][k] = temp;
		}
	/* Zero entries in the columns of AUG */
	ip1 = i+1;
	for(k=ip1;k<n;k++)
		{
		Q = -AUG[k][i]/AUG[i][i];
		AUG[k][i] = 0.0;
		for(j=ip1;j<np1;j++)
			{
			AUG[k][j] = Q*AUG[i][j]+AUG[k][j];
			}
		}
	if(DEBUG)
	{
	sprintf(s,"After Processing Row %d",i);
	print_matrix(AUG,n,np1,s);
	}
	}
	if(fabs(AUG[nm1][nm1]) < VERY_SMALL) 
		{ 
		(*ierror) = 2;
		printf("Error: input matrix is singular\n");
		return;
		}


if(DEBUG) print_matrix(AUG,n,np1,"After Triangulation");

/* Backsolve to obtain the solution to MX=B. */
X[nm1] =  AUG[nm1][n]/AUG[nm1][nm1];
for(k=1;k<=nm1;k++)
	{
	Q = 0.0;
	for(j=0;j<k;j++)
		{
		Q = Q+AUG[nm1-k][nm1-j]*X[nm1-j];
		}
	X[nm1-k] = (AUG[nm1-k][n]-Q)/AUG[nm1-k][nm1-k];
	}
/* Calculate the norm of the residual error, B-AX */
RSQ = 0.0;
for(i=0;i<n;i++)
	{
	Q = 0.0;
	for(j=0;j<n;j++)
		Q = Q+A[i][j]*X[j];
	residual = B[i] - Q;
	Rmag = fabs(residual);
	if(DEBUG) printf("i=%d residual=%f\n",i,Rmag);
	RSQ = RSQ+Rmag*Rmag;
	}
(*rnorm) = sqrt(RSQ);
(*ierror) = 1;
}