/*
    gep.c
    Steve Cox, October 1, 2000
    A driver and routines for solving a general linear system
    The system is a dynamically allocated one dimensional
    array of pointers. The right hand side is overwritten
    with the solution. All arrays begin indexing from 0.
*/

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

void trisolve(double *[], int);
void rzero(double *[], int , int , int , int );
int gauss(double *[], int);
void matview(double *[], int);

main()
{

int j, k,		/* counters */
    n, 			/* dimension of A */
    info;		/* gauss indicator */

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

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

for (j=0;j<n;j=j+1)     /* allocate space */
    ptr[j] = (double *)malloc( (n+1)*sizeof(double) );

printf("\n Please enter your matrix (row wise):\n\n");
for (j=0;j<n;j=j+1){
for (k=0;k<n;k=k+1){
    scanf("%lf", ptr[j]+k);
}
}
printf("\nPlease enter the components of the right hand side\n");
for (j=0;j<n;j=j+1)
    scanf("%lf",ptr[j]+n);
 
printf("The system is\n");
matview(ptr,n);

info = gauss(ptr,n); 	/* solve the system and place the answer in
			   "column" n of ptr */

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

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

    }
    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 zero (or small) pivot then return the
	respective row number, else return 0
*/

int gauss(double *ptr[], int n)

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

double max;		/* fabs of dominant element */
double *tmp;

k = 0;
while(k<n-1){

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

    	  if ( fabs(ptr[j][k]) > max ){ 
             max = fabs(ptr[j][k]);
             r1 = j;
   	  }
 
      }

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

     tmp = ptr[k];		/* swap row#k and row#r1 */
     ptr[k] = ptr[r1];
     ptr[r1] = tmp;
     printf("after swapping\n");
     matview(ptr,n);
             
     for (j=k+1;j<n;j=j+1){	/* use max to make zeros in all */
 	 rzero(ptr,n+1,k,k,j);	/* places below the diagonal */
         printf("after rzero \n");
         matview(ptr,n);
     }

     k = k + 1;

}	/* close while */

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

trisolve(ptr,n);		/* a is now triangular, so use */
					/* trisolve */

return 0;			/* return all clear signal */

}


void trisolve(double *ptr[], int n)

/* 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


			*/

        ptr[n-1][n] = ptr[n-1][n]/ptr[n-1][n-1];

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

/* 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 *ptr[], 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 = -ptr[r2][c1]/ptr[r1][c1];

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

void matview(double *ptr[], int n)
{
int j,k;
for (j=0;j<n;j=j+1){
for (k=0;k<=n;k=k+1){
    printf("%lf ",ptr[j][k]);
    }
    printf("\n");
}
}

/*

Script started on Tue Oct 23 08:36:06 2001
harvey.caam.rice.edu% a.out

How large is your problem? 3


 Please enter your matrix (row wise):

2 -1 0
-1 2 -1
0 -1 2

Please enter the components of the right hand side
1 1 1
The system is
2.000000 -1.000000 0.000000 1.000000 
-1.000000 2.000000 -1.000000 1.000000 
0.000000 -1.000000 2.000000 1.000000 
after swapping
2.000000 -1.000000 0.000000 1.000000 
-1.000000 2.000000 -1.000000 1.000000 
0.000000 -1.000000 2.000000 1.000000 
after rzero 
2.000000 -1.000000 0.000000 1.000000 
0.000000 1.500000 -1.000000 1.500000 
0.000000 -1.000000 2.000000 1.000000 
after rzero 
2.000000 -1.000000 0.000000 1.000000 
0.000000 1.500000 -1.000000 1.500000 
0.000000 -1.000000 2.000000 1.000000 
after swapping
2.000000 -1.000000 0.000000 1.000000 
0.000000 1.500000 -1.000000 1.500000 
0.000000 -1.000000 2.000000 1.000000 
after rzero 
2.000000 -1.000000 0.000000 1.000000 
0.000000 1.500000 -1.000000 1.500000 
0.000000 0.000000 1.333333 2.000000 
The solution vector is

1.500000
 2.000000
 1.500000
 
harvey.caam.rice.edu% a.out

How large is your problem? 3


 Please enter your matrix (row wise):

1 1 1
2 2 2
3 3 3

Please enter the components of the right hand side
1 1 1
The system is
1.000000 1.000000 1.000000 1.000000 
2.000000 2.000000 2.000000 1.000000 
3.000000 3.000000 3.000000 1.000000 
after swapping
3.000000 3.000000 3.000000 1.000000 
2.000000 2.000000 2.000000 1.000000 
1.000000 1.000000 1.000000 1.000000 
after rzero 
3.000000 3.000000 3.000000 1.000000 
0.000000 0.000000 0.000000 0.333333 
1.000000 1.000000 1.000000 1.000000 
after rzero 
3.000000 3.000000 3.000000 1.000000 
0.000000 0.000000 0.000000 0.333333 
0.000000 0.000000 0.000000 0.666667 
Sorry, no reliable solution
A very small pivot was found in row 1

harvey.caam.rice.edu% ^D
script done on Tue Oct 23 08:36:52 2001

*/