/* GaussianMatrix.cpp defines non-trivial GaussianMatrix operations.
 *
 */

#include "GaussianMatrix.h"

void GaussianMatrix::readCoefficientMatrix(istream& in)
{
   int numEquations = rows();
   
   for (int r = 0; r < numEquations; r++)
     for (int c = 0; c < numEquations; c++)
      in >> (*this)[r][c];
}

void GaussianMatrix::readConstantVector(istream& in)
{  
   int numEquations = rows();
   for (int r = 0; r < numEquations; r++)
     in >> (*this)[r][numEquations];
}

inline void exchange(double& a, double& b)
{
  double t = a; a = b; b = t;
}
  
bool GaussianMatrix::reduce()
{
  const double EPSILON = 1.0E-6;
  bool isSingular = false;
  int i = 0,
      j,
      k,
      numRows = rows(),
      pivotRow;
  double quotient,
         absolutePivot;
  
  while ((!isSingular) && (i < numRows))
  {
    absolutePivot = abs( (*this)[i][i] );
    pivotRow = i;
    for (k = i+1; k < numRows; k++)
      if (abs( (*this)[k][i] ) > absolutePivot)
      {
         absolutePivot = abs( (*this)[k][i] );
         pivotRow = k;
      }
      
    isSingular = absolutePivot < EPSILON;
    
    if (!isSingular)
    {
       if (i != pivotRow)
         for (j = 0; j <= numRows; j++)
            exchange( (*this)[i][j], (*this)[pivotRow][j] );
 
       for (j = i+1; j < numRows; j++)
       {
          quotient = -(*this)[j][i] / (*this)[i][i];
          for (k = i; k <= numRows; k++)
            (*this)[j][k] = (*this)[j][k] + quotient * (*this)[i][k];
        }
    }
    i++;
  }
  return isSingular;
}

Matrix GaussianMatrix::solve() const
{
  int n = rows()-1;
  Matrix solutionVector(1, n+1);
 
  solutionVector[0][n] = (*this)[n][n+1] / (*this)[n][n];
  
  for (int i = n-1; i >= 0; i--)
  {
    solutionVector[0][i] = (*this)[i][n+1];
      
    for (int j = i+1; j <= n; j++)
      solutionVector[0][i] -= (*this)[i][j] * solutionVector[0][j];
 
    solutionVector[0][i] /= (*this)[i][i];
  }
  
  return solutionVector;
}