/* Computer Soft/c2-1.c    Gauss Elimination */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define TRUE 1
/*             a[i][j] : matrix element, a(i,j)
               n : order of matrix
               eps : machine epsilon
               det : determinant                                */
void main()
{
int i, j, _i, _r;
static n = 3;
static float a_init[10][11] = {{1, 2, 3,  6},
                               {2, 2, 3,  7},
                               {3, 3, 3,  9}};
static double a[10][11];
void gauss();
/*static int _aini = 1; */

   printf( "\nComputer Soft/C2-1    Gauss Elimination \n\n" );
   printf( "Augmented matrix\n" );
   for( i = 1; i <= n; i++ ){   
      for( j = 1; j <= n+1; j++ ) {   
	 a[i][j]=a_init[i-1][j-1];  printf(  "  %13.5e", a[i][j] );
      }
      printf(  "\n" );
   }    
   gauss( n, a );
   printf( " Solution\n" );
   printf( "-----------------------------------------\n" );
   printf( "       i        x(i)\n" );
   printf( "-----------------------------------------\n" );
   for( i = 1; i <= n; i++ ) printf( "   %5d  %16.6e\n", i, a[i][n+1] );
   printf( "-----------------------------------------\n\n" );
   exit(0);
}   
         
void gauss(n, a)
int n;   double a[][11];
{
int i, j, jc, jr, k, kc, nv, pv;
double det, eps, ep1, eps2, r, temp, tm, va;
   eps = 1.0; ep1 = 1.0 ;           /* eps = Machine epsilon */
   while( ep1  > 0 ){ 
      eps = eps/2.0; ep1 = eps*0.98 + 1; ep1 = ep1 - 1;
   }    
   eps = eps*2;       eps2 = eps*2;
   printf( "                Machine epsilon=%g \n", eps );
   det = 1;                      /* Initialization of determinant */
   for( i = 1; i <= (n - 1); i++ ){   
      pv = i;
      for( j = i + 1; j <= n; j++ ){   
         if( fabs( a[pv][i] ) < fabs( a[j][i] ) )   pv = j;
      }
      if( pv != i ){   
         for( jc = 1; jc <= (n + 1); jc++ ){   
            tm = a[i][jc];  a[i][jc] = a[pv][jc];     a[pv][jc] = tm;
         }
         det = -det;
      }
      if( a[i][i] == 0 ){       /* Singular matrix */
         printf( "Matrix is singular.\n" );    exit(0);
      }
      for( jr = i + 1; jr <= n; jr++ ){ /* Elimination of below-diagonal. */
         if( a[jr][i] != 0 ){   
            r = a[jr][i]/a[i][i];
            for( kc = i + 1; kc <= (n + 1); kc++ ){   
               temp = a[jr][kc];
               a[jr][kc] = a[jr][kc] - r*a[i][kc];
               if( fabs( a[jr][kc] ) < eps2*temp ) a[jr][kc] = 0.0;
/*                           If the result of subtraction is smaller than
 *                           2 times machine epsilon times the original
 *                           value, it is set to zero. */
         
            }
         }
      }
   }    
   for( i = 1; i <= n; i++ ) {   
      det = det*a[i][i];           /*  Determinant is calculated. */
   }    
   if( det == 0 ){
      printf( "Matrix is singular.\n" );    exit(0);
   }
   else{                         /* Backward substitution starts. */  
      a[n][n+1] = a[n][n+1]/a[n][n];
      for( nv = n - 1; nv >= 1; nv-- ){   
         va = a[nv][n+1];   
         for( k = nv + 1; k <= n; k++ ) {va = va - a[nv][k]*a[k][n+1];}
         a[nv][n+1] = va/a[nv][nv];
      }
      printf( "                Determinant = %g \n", det );
      return;
   }    
}