#include <iostream.h>

#define N 3

float A[3][3] = { 
//	        { 2, 2, 4 },	// This matrix will produce 
//	        { 2, 4, 6 },	// virtually no rounding
//	        { 2, 4, 5 }	// error

	        { 2, 2, 4 },		// This matrix will produce
	        { 2, 2.0000002, 6 },	// a very large rounding
	        { 3, 3.0000003, 5 }	// error
	      };

float b[3] = { 2, 4, 3 };


void Print()
{
   int i, j;

   for (i = 0; i < N; i = i+1 )
   {
      for (j = 0; j < N; j = j+1 )
	 cout << A[i][j] << "\t";
      cout << "=\t" << b[i] << "\n";
   }
   cout << "\n";
}


int main(int argc, char *argv[])
{
   int i, j, k, l;
   float pivot, v;

   float A_Orig[3][3];
   float b_Orig[3];


   for (i = 0; i < N; i = i+1 )
   {
      for (j = 0; j < N; j = j + 1)
         A_Orig[i][j] = A[i][j];
      b_Orig[i] = b[i];
   }


   for ( i = 0; i < N; i = i + 1 )
   {
      Print();

      pivot = A[i][i];
      for (j = i; j < N; j = j + 1 )
         A[i][j] = A[i][j]/pivot;
      b[i] = b[i]/pivot;

      Print();
      for (k = 0; k < N; k = k + 1)
         if ( k != i )
         {
            pivot = A[k][i];
            for (l = i; l < N; l = l + 1)
	        A[k][l] = A[k][l] - pivot*A[i][l];
	    b[k] = b[k] - pivot*b[i];
         }

      Print();

      cout << "=========\n\n";
   }


   for (i = 0; i < N; i = i+1 )
   {
      v = 0.0;
      for (j = 0; j < N; j = j + 1)
         v = v + A_Orig[i][j]*b[j];
      cout << "solution value = " << v << " - should be = " 
	   << b_Orig[i] << "\n";
   }


}