/*  main is a driver for a general linear solver called gauss */

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

void trisolve(double [][100], int, double [], double []);

void rswap(double [][100], int , int , int );    

void rzero(double [][100], int , int , int , int );

int gauss(double [][100],double [],int ,double []);

main()
{
  /* DRIVER: get the vector, b, 
	     get the matrix, A,
	     solve A x = b,
	     display x or inform user of reason for lack of x */

        int j, k,		/* counters */
	    n, 			/* dimension of A */
	    info;		/* gauss indicator */
 
        double b[100],		/* the right hand side */
	      x[100];		/* the solution */

        double A[100][100];	/* the matrix */

        printf("\n");
        printf("How large is your problem? ");
        scanf("%d",&n);
        printf("\n");

        printf("Please enter the components of the right hand side\n");
        for (j=1;j<=n;j=j+1)
                scanf("%lf",&b[j]);

        printf("\n");
        printf("\n");
        for (j=1;j<=n;j=j+1){
            printf("Please enter row %d of the matrix\n",j);
            for (k=1;k<=n;k=k+1){
                scanf("%lf",&A[j][k]);
            }
        }
	printf("\n");
 
        info = gauss(A,b,n,x);

	if (info==0) {		/* gauss succeeded */

		printf("The solution vector is\n");
		printf("\n");
	        for (j=1;j<=n;j=j+1)
        	     printf("%lf\n ",x[j]);

        	}
	else {			/* gauss failed */
    
		printf("Sorry, no reliable solution\n");
		printf("A very small pivot was found in row %d\n",info);
        }

    printf("\n");

}

/* GAUSSIAN ELIMINATION : 

	reduce A x = b to triangular form and call trisolve.

   arguments:  A	the matrix 
	       b	the right hand side
	       n	the size of the matrix
	       x	the solution to Ax=b

   returned:

	if the reduction meets a small pivot then return the
	respective row number, else return 0
*/

int gauss(double A[][100],double b[],int n,double x[])

{
	int j, k,		/* counters */
	    index;		/* index of dominant row */

	double max;		/* fabs of dominant element */

	for (k=1;k<=n;k=k+1)		/* augment A with b */
	    A[k][n+1] = b[k];
     
	for (j=1;j<n;j=j+1) {  /* j counts columns */ 

	     /* compute index, the row number, greater than or
	        equal to j, in column j, where the element of
		greatest magnitude resides */
     
	     max = fabs(A[j][j]);   /* set max to fabs of diagonal element */
	     index = j;
     
	     for (k=j+1;k<=n;k=k+1) {  /* search down and reset max and
	  				  index when appropriate */

       		  if ( fabs(A[k][j]) > max ){ 
		           max = fabs(A[k][j]);
		           index = k;
 	   	  }
	   
             }

	     if (max < 1e-12) 		/* if max is small return the bad j */
		     return j;

      	     rswap(A,n+1,index,j);	/* swap row#j and row#index */

	     for (k=j+1;k<=n;k=1+k)	/* use max to make zeros in all */
   	         rzero(A,n+1,j,j,k);	/* places below the diagonal */

	}

	if ( fabs(A[n][n]) < 1e-12)	/* if final max is small return n */
	  return n;

	for (k=1;k<=n;k=k+1)		/* strip off the modified b */
	    b[k] = A[k][n+1];

	trisolve(A,n,b,x);		/* A is now triangular, use trisolve */

	return 0;			/* return all clear signal */

}


void trisolve(double U[][100], int n, double b[], double x[])

/* arguments:  U	the upper triangular matrix 
	       n	the size of the matrix
	       b	the right hand side
	       x	the solution to Ux=b
*/

{
	int j,		/* row counter */
	    k;		/* column counter */
	double tmp;	/* temporary storage, note, for fixed j
			   we want

			         k=n
			tmp = sum     U[j][k]*x[k]
				 k=j+1


			*/

	x[n] = b[n]/U[n][n]; 	/* solve at last row first */

	for (j=n-1;j>=1;j=j-1){		/* now back up through the rows */
            tmp = 0;
	    for (k=j+1;k<=n;k=k+1){
	    	tmp = tmp + U[j][k]*x[k];
	    }
            x[j] = (b[j] - tmp)/U[j][j];
	}
}

/* ROW SWAPPER 
swap two rows in a matrix  */
    
void rswap(double a[][100], int n, int r1, int r2)

/* arguments:  a	the matrix in need of a swapping
	       n	the size of the matrix
	       r1	one of the row numbers
	       r2	the other row number
*/

{
	int j;		/* row counter */
	double tmp;	/* temporary storage */

	for (j=1;j<=n;j=j+1){
	    tmp = a[r1][j];
	    a[r1][j] = a[r2][j];
            a[r2][j] = tmp;
	}
}

  
/* ZERO-MAKER
    eliminate the nonzero in column c1 row c2 by performing 
    the following elementary row operation:

	row r2 = row r2 + (mult * row r1)

	where mult = -a[r2][c1]/a[r1][c1]
*/
    
void rzero(double a[][100], int n, int c1, int r1, int r2)

/* arguments:  a	the matrix 
	       n	the size of the matrix
	       c1	the column in which elimination occurs
	       r1	the dominant row
	       r2	the submissive row
*/

{
	int j;		/* column counter */
        double mult;	/* multiplier */

	mult = -a[r2][c1]/a[r1][c1];

	for (j=1;j<=n;j=j+1)
            a[r2][j] = a[r2][j] + mult*a[r1][j];
}