/*
  Gaussian Elimination

  Copyright (c) Paul Bourke, 1997 
  http://astronomy.swin.edu.au/pbourke/analysis/gausselim/

  Modified for use in Mandarin Tone Recognizance software
  by Jim Gilsinan IV
  http://www.fas.harvard.edu/~gilsinan/

  Senior Thesis
  Harvard College, 2001

  Solve a system of n equations in n unknowns using Gaussian Elimination
  Solve an equation in matrix form Ax = b
  The 2D array a is the matrix A with an additional column b.
  This is often written (A:b)
  
  A0,0    A1,0    A2,0    ....  An-1,0     b0
  A0,1    A1,1    A2,1    ....  An-1,1     b1
  A0,2    A1,2    A2,2    ....  An-1,2     b2
  :       :       :             :          :
  :       :       :             :          :
  A0,n-1  A1,n-1  A2,n-1  ....  An-1,n-1   bn-1
  
  The result is returned in x, otherwise the function returns 1
  if the system of equations is singular.
*/

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

#define EPS 0.00001

int GSolve(double **a, int n, double *x)
{
   int i, j, k, maxrow;
   double tmp;

   for (i = 0; i < n; i++) {

      /* Find the row with the largest first value */
      maxrow = i;
      for (j = i + 1; j < n; j++)
         if (fabsf(a[i][j]) > fabsf(a[i][maxrow]))
            maxrow = j;

      /* Swap the maxrow and ith row */
      for (k=i;k<n+1;k++) {
         tmp = a[k][i];
         a[k][i] = a[k][maxrow];
         a[k][maxrow] = tmp;
      }

      /* Singular matrix? */
      if (fabsf(a[i][i]) < EPS)
         return 1;

      /* Eliminate the ith element of the jth row */
      for (j=i+1;j<n;j++) {
         for (k=n;k>=i;k--) {
            a[k][j] -= a[k][i] * a[i][j] / a[i][i];
         }
      }
   }

   /* Do the back substitution */
   for (j=n-1;j>=0;j--) {
      tmp = 0;
      for (k=j+1;k<n;k++)
         tmp += a[k][j] * x[k];
      x[j] = (a[n][j] - tmp) / a[j][j];
   }

   return 0;
}