/* gaussElim.cpp is a driver program for Gaussian Elimination.
 * 
 * Input:  a series of linear equation coefficients
 * Output: the solution of the linear equation or a "singular
 *           system" message
 * Joel Adams, Fall 1996 at Calvin College.
 ***************************************************************/
 
#include <iostream.h>     // cout, cin, ...
#include <stdlib.h>       // exit()
#include "Matrix.h"       // Matrix class
 
void ReadEquations(Matrix & augMat);
int Reduce(Matrix & augMat);
Matrix Solve(Matrix & augMat);
Matrix GaussElim();
 
int main()
{
  cout << "This program solves a linear system "
          "using Gaussian Elimination. \n";
 
  Matrix solutionVector = GaussElim();
  
  cout << "\nThe solution vector (x1, x2, ...)  for this system is:\n\n"
       << solutionVector << endl;
}

/* GaussElim() performs the Gaussian Elimination algorithm
 *
 * Input:  the coefficients of a linear system and its constant vector
 * Return: the solution to the linear system
 * Joel Adams, Fall 1996 at Calvin College.
 ************************************************************************/
  
Matrix GaussElim()
{
   Matrix augmentedMatrix;
  
   ReadEquations(augmentedMatrix);
  
   bool isSingular = Reduce(augmentedMatrix);
 
   if (isSingular)
   {
      cerr << "\n*** GaussElim: Coefficient Matrix is (nearly) singular!\n";
      exit (0);
   }
 
   Matrix solutionVector = Solve(augmentedMatrix);
 
   return solutionVector;
}
 
 
void ReadEquations(Matrix & augMat)
{
  int numEquations;
 
  for (;;)
    {
      cout << "\nPlease enter the number of equations in the system: ";
      cin >> numEquations;
      if (numEquations > 1) break;
      cerr << "\n*** At least two equations are needed ...\n";
    }
 
  augMat = Matrix(numEquations,       // numEquations rows
                  numEquations+1);    // numEquations+1 columns
 
  cout << "\nPlease enter the coefficient matrix row-by-row...\n";
  for (int r = 0; r < numEquations; r++)
    for (int c = 0; c < numEquations; c++)
      cin >> augMat[r][c];
 
  cout << "\nPlease enter the constant vector...\n";
  for (int r = 0; r < numEquations; r++)
    cin >> augMat[r][numEquations];
}
 
 
inline double Abs(double val)
{
  return (val < 0) ? -(val) : val;
}
 
inline void Swap(double & a, double & b)
{
  double t = a; a = b; b = t;
}
  
int Reduce(Matrix & augMat)
{
  const double EPSILON = 1.0E-6;
  bool isSingular = false;
  int i = 0,
      j,
      k,
      numRows = augMat.Rows(),
      pivotRow;
  double quotient,
         absolutePivot;
  
  while ((!isSingular) && (i < numRows))
  {
    absolutePivot = Abs(augMat[i][i]);
    pivotRow = i;
    for (k = i+1; k < numRows; k++)
      if (Abs(augMat[k][i]) > absolutePivot)
      {
         absolutePivot = Abs(augMat[k][i]);
         pivotRow = k;
      }
    isSingular = absolutePivot < EPSILON;
    if (!isSingular)
    {
       if (i != pivotRow)
         for (j = 0; j <= numRows; j++)
           Swap(augMat[i][j], augMat[pivotRow][j]);
 
       for (j = i+1; j < numRows; j++)
       {
          quotient = -augMat[j][i] / augMat[i][i];
          for (k = i; k <= numRows; k++)
            augMat[j][k] = augMat[j][k] + quotient * augMat[i][k];
        }
    }
    i++;
  }
  return isSingular;
}
         
Matrix Solve(Matrix & augMat)
{
  Matrix solutionVector(1, augMat.Rows());
  int n = augMat.Rows()-1;
 
  solutionVector[0][n] = augMat[n][n+1] / augMat[n][n];
  
  for (int i = n-1; i >= 0; i--)
  {
    solutionVector[0][i] = augMat[i][n+1];
      
    for (int j = i+1; j <= n; j++)
      solutionVector[0][i] -= augMat[i][j] * solutionVector[0][j];
 
    solutionVector[0][i] /= augMat[i][i];
  }
  
  return solutionVector;
}