/*

                  Library program for solving simultaneous linear equations
         Programmer:          Larry Caretto
         Version:             2.02
         Date:                March 18, 2003

         Instructions for use

The boolean function GaussianElim solves the following set of simultaneous linear equations:

               a[i][1]*x[1] + a[i][2]*x[2] + ..........+ a[i][N] = b[i]   i = 1, N

N-squared values of the coefficients a[i][j] and N values of the b[i] coefficients are
passed as parameters to the routine by use of the arrays a and b.

      a is a two-dimensional, type double, array containg the coefficients a[i][j] of
        the linear equations
      x is a one-dimensional, type double, array containing the unknowns to be found and
      b is a one-dimensional, type double, array containing the coefficients of the right hand
        sides of the simultaneous linear equations

The function is invoked by a statement like the following:

                     errorFlag = GaussianElim( a, b, x, N )

where

   a is the double array variable containing the equation coefficients.  It must
     have a dimension of [MAX_VAR][MAX_VAR].  The value of MAX_VAR must be set in
     a const statement in the code that calls GaussianElim and in this code.  The
     values of MAX_VAR in the two codes must match.  The number of variables to be
     found must be < MAX_VAR

   b is the double one-dimensional array which contains the right hand sides of
     the equation set

   x is the double one-dimensional array which contains the solutions

   N is an integer variable giving the number of equations; this must be < MAX_VAR
     Note that the requirement that the number of variables must be less 
than MAX_VAR
     is necessary becasue an extra column is added to the a array during 
the solution

   The boolean value returned by GaussianElim is an error flag.  It has the following meaning:

        TRUE - the input a array does not have a valid unique solution; this is
               known as a singular matrix.  Data in the x array are meaningless.

        FALSE - the x array contains a valid solution.

The following C++ program illustrates the use of the GaussianElim routine.  Note
that the cmath library must be included for GaussianElim which requires the
use of the function fabs.

   #include <iostream>
   #include <cmath>
   #include <iomanip>
   using namespace std;

   int const MAX_VAR = 20;

   int main()
   {

      bool GaussianElim( double a[][MAX_VAR], double b[],
                         double x[], int N_eqn );          // function prototype

      double a[MAX_VAR][MAX_VAR],
             b[MAX_VAR],
             x[MAX_VAR];

      int    i, N_eqn = 3;

      a[0][0] = 1;
      a[0][1] = 2;
      a[0][2] = 3;
      a[1][0] = 2;
      a[1][1] = 4;
      a[1][2] = 7;
      a[2][0] = 3;
      a[2][1] = 0;
      a[2][2] = 1;
      b[0] = 8;
      b[1] = 18;
      b[2] = 2;
      if ( GaussianElim( a, b, x, N_eqn ) )
         cout << "\n\n* * * SINGULAR MATRIX,  No unique solution possible * * *";
      else
      {
         cout << "\n\nSolution of Linear Equations: \n";
         for ( i = 0; i < N_eqn; i++)
         cout << "\n   x[" << setw(2) << i << "] = " << x[i];
      }
      a[0][0] = 1;
      a[0][1] = 2;
      a[0][2] = 3;
      a[1][0] = 2;
      a[1][1] = 4;
      a[1][2] = 6;
      a[2][0] = 3;
      a[2][1] = 0;
      a[2][2] = 1;
      if ( GaussianElim( a, b, x, N_eqn )  )
          cout << "\n\n* * * SINGULAR MATRIX,  No unique solution possible * * *";
      else
      {
         cout << "\n\nSolution of Linear Equations: \n";
         for ( i = 0; i < N_eqn; i++)
            cout << "\n   x[" << setw(2) << i << "] = " << x[i];
      }
      cout << endl << endl;
      return EXIT_SUCCESS;

}

This sample program produces the output shown below.  The value of x[0] should be
zero for the first case.  The program output, a small number, is due to roundoff
error.  The value may be different on other computers.  The second case is obtained
by changing a[1][2] to 6; this gives a singular matrix.

Solution of Linear Equations:

   x[ 0] = -5.92119e-16
   x[ 1] = 1
   x[ 2] = 2

* * * SINGULAR MATRIX,  No unique solution possible * * *

Press any key to continue:
*/
#include <cmath>
const int MAX_VAR = 20;   // This must be set to be the same for both this program and
                          // the calling program.  If both programs are placed in one file,
                          // a common global variable can be used.

bool GaussianElim( double a[][MAX_VAR],     // Left-hand-side matrix
                   double b[],                // Right-hand-side column matrix
                   double x[],                // Column matrix for solution
                   int N_eqn )
{

   int       N_rhs = 1;     // Number of right hand side vectors

   void      AugmentMatrix( double a[][MAX_VAR], double b[],
                           int N_eqn, int N_rhs );

   double    GetMax( double a[][MAX_VAR], int N_eqn );   // function to get maximum array element

   bool      swaprows( double a[][MAX_VAR], int pivot, int N_eqn, int N_rhs, 
                       double max_in_matrix ); // Subroutine to swap pivot rows


            //   Start of executable statements

   double max_in_matrix = GetMax( a, N_eqn );   // Determine maximum absolute value in matrix

   AugmentMatrix( a, b, N_eqn, N_rhs );

               //   Reduction to upper triangular matrix starts here
   
   bool singular;   // flag to detect singular matrix (has no solution)

   for ( int pivot = 0; pivot <= N_eqn-1; pivot++ )
   {
         if ( (singular = swaprows( a, pivot , N_eqn, N_rhs, max_in_matrix ) ) )
			     break;
         for ( int row = pivot + 1; row < N_eqn; row++ )
               for ( int column = pivot+1; column < N_eqn + N_rhs; column++ )
                  a[row][column] -= a[row][pivot] * a[pivot][column]
                                  / a[pivot][pivot];
   }

      //   Now do back substitution (unless a singluar matrix is found)

   if ( !singular )
   {
      for ( int rhs = N_eqn; rhs < N_eqn + N_rhs; rhs++ )
         for ( int row = N_eqn-1; row >=0; row-- )
         {
            for ( int column = row+1; column < N_eqn; column++ )
               a[row][rhs] -= a[row][column] * a[column][rhs];

            a[row][rhs] /= a[row][row];
         }

         for ( int row = 0; row < N_eqn; row++ )
            x[row] = a[row][N_eqn];

   }

   return singular;

}

void AugmentMatrix( double a[][MAX_VAR], double b[],
                    int N_eqn, int N_rhs )

{
      //   Add right-hand side vectors as extra columns 
      //   in origina a matrix

   for (int row=0; row<N_eqn; row++)
      for (int column=N_eqn; column<N_eqn+N_rhs; column++)
         a[row][column] = b[row];
}


bool swaprows( double a[][MAX_VAR], int pivot, int N_eqn, int N_rhs,
               double max_in_matrix )
{
/*
      Find row with maximum element in the pivot column and swap this row
      with the current pivot row.  If the resulting pivot element is
      smaller than some desired minimum, the matrix is assumed to be singular.
      In this case the routine returns the value true as an error flag.  If there is
      no error (the calculation can proceed), the function returns the value false.
*/

   const double small = 1e-9;  // The matrix is considered to be a singular matrix if the
                               // maximum element in the pivot column is less than "small"
                               // times the maximum element in the array 

         // Search for maximum element in pivot column starting with pivot row and
         // checking all rows from pivot row to last row.

   double max_element = fabs( a[pivot][pivot] );
   int max_row = pivot;
   for (int row = pivot+1; row < N_eqn; row++ )
      if (  fabs( a[row][pivot] ) > max_element)
      {
         max_row = row;
         max_element = fabs( a[row][pivot] );
      }

      // Maximum element found.  Check to see if this is essentially zero to within
      // specified tolerance for roundoff error.  If it is, return a true error flag.
      // If it is not, swap rows if required and return a error flag value of false.

   if ( max_element < small * max_in_matrix )
         return true;
   else
   {
      if ( max_row != pivot)
         for ( int column = pivot; column < N_eqn + N_rhs; column++)
        {
             double temp = a[max_row][column];
             a[max_row][column] = a[pivot][column];
             a[pivot][column] = temp;
        }
      return false;

   }

}

double   GetMax( double a[][MAX_VAR], int N_eqn )
{
                             // Find matrix element with maximum absolute value

   double max_in_matrix = 0.0;
   for ( int row = 0; row < N_eqn; row++ )
      for ( int column = 0; column < N_eqn; column++ )
         if ( max_in_matrix < fabs( a[row][column] ) )
            max_in_matrix = fabs( a[row][column] );


   return max_in_matrix;
}