/*
*   GAUSSIAN ELIMINATION WITH BACKWARD SUBSTITUTION ALGOTITHM 6.1
*
*   To solve the n by n linear system
*
*   E1:  A[1,1] X[1] + A[1,2] X[2] +...+ A[1,n] X[n] = A[1,n+1]
*   E2:  A[2,1] X[1] + A[2,2] X[2] +...+ A[2,n] X[n] = A[2,n+1]
*   :
*   .
*   EN:  A[n,1] X[1] + A[n,2] X[2] +...+ A[n,n] X[n] = A[n,n+1]
*
*   INPUT:   number of unknowns and equations n; augmented
*            matrix A = (A(I,J)) where 1<=I<=n and 1<=J<=n+1.
*
*   OUTPUT:  solution x(1), x(2),...,x(n) or a message that the
*            linear system has no unique solution.
*
*    Example 1 of Section 6.1 is solved by this algorithm. 
*/



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

#define ZERO 1.0E-20
#define true 1
#define false 0

#define N  4   /* Number of equations */

static double absval(double);

using namespace std;

int main()
{
	double A[N][N+1] = { {1.0, 1.0, 0.0, 3.0, 4.0},
	                   {2.0, 1.0, -1.0, 1.0, 1.0},
	                   {3.0, -1.0, -1.0, 2.0, -3.0},
	                   {-1.0, 2.0, 3.0, -1.0, 4.0}
	                 };
	double X[N];

   double C,XM,SUM;
   int M, I,J,ICHG,NN,IP,JJ,K,L,KK,OK = true;


      /* STEP 1 */
      NN = N - 1;
      M = N + 1;
      ICHG = 0;
      I = 1;

      while ((OK) && (I <= NN)) 
      {
         /* STEP 2 */
         /* use IP instead of p */
         IP = I;
         while ((absval(A[IP-1][I-1]) <= ZERO) && (IP <= N))
            IP++;
         if (IP == M) OK = false;
         else {
            /* STEP 3 */
            if (IP != I) {
               for (JJ=1; JJ<=M; JJ++) {
                  C = A[I-1][JJ-1];
                  A[I-1][JJ-1] = A[IP-1][JJ-1];
                  A[IP-1][JJ-1] = C;
               }
               ICHG = ICHG + 1;
            }  
            /* STEP 4 */
            JJ = I + 1; 
            for (J=JJ; J<=N; J++) {
               /* STEP 5 */
               /* use XM in place of m(J,I) */
               XM = A[J-1][I-1] / A[I-1][I-1]; 
               /* STEP 6 */
               for (K=JJ; K<=M; K++)
                  A[J-1][K-1] = A[J-1][K-1] - XM * A[I-1][K-1];
               /*  Multiplier XM could be saved in A[J,I].  */
               A[J-1][I-1] = 0.0;
            }  
         }
         I = I + 1;
      }

      if (OK) 
      {
         /* STEP 7 */
         if (absval(A[N-1][N-1]) <= ZERO) OK = false;
         else 
	 {
            /* STEP 8 */
            /* start backward substitution */
            X[N-1] = A[N-1][M-1] / A[N-1][N-1];
            /* STEP 9 */
            for (K=1; K<=NN; K++) {
               I = NN - K + 1;
               JJ = I + 1;
               SUM = 0.0;
               for (KK=JJ ; KK<=N; KK++)  
                  SUM = SUM - A[I-1][KK-1] * X[KK-1];
               X[I-1] = (A[I-1][M-1]+SUM) / A[I-1][I-1];
            }

            /* STEP 10 */
            /* procedure completed successfully */
            //OUTPUT(N, M, ICHG, X, A);
            cout<<"GAUSSIAN ELIMINATION"<<endl;
	    cout<<"The reduced system - output by rows:"<<endl;
	    for (I=0; I<N; I++) 
	    {
                for (J=0; J<M; J++)
		{
	   	    cout<<A[I][J] << " ";
		}
		cout<<endl;
	    }

	    cout<<endl<<"Solution vector:"<<endl;
	    for (I=0; I<N; I++) 
		cout<< X[I] <<" ";
	    cout<<endl <<endl;

	    cout<<ICHG << " row interchanges have been performed"<<endl;
         }  
      }

      if (!OK) 
      {
	  cout<<"System has no unique solution"<<endl;
      }

   return 0;

}


/* Absolute Value Function */
static double absval(double val)
{
   if (val >= 0) return val;
   else return -val;
}