/*
*   GAUSS-SEIDEL ITERATAIVE TECHNIQUE ALGORITHM 7.2
*
*   To solve Ax = b given an initial approximation x(0).
*   Example 7.3.3
*
*   INPUT:   the number of equations and unknowns n; the entries
*            A(I,J), 1<=I, J<=n, of the matrix A; the entries
*            B(I), 1<=I<=n, of the inhomogeneous term b; the
*            entries XO(I), 1<=I<=n, of x(0); tolerance TOL;
*            maximum number of iterations N.
*
*    OUTPUT: the approximate solution X(1),...,X(n) or a message
*            that the number of iterations was exceeded.
*/

#include<stdio.h>
#include<math.h>
#include<iostream>
#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] = { {10.0, -1.0, 2.0, 0.0, 6.0}, 
	                     {-1.0, 11.0, -1.0, 3.0, 25.0},
	                     {2.0, -1.0, 10.0, -1.0, -11.0},
	                     {0.0, 3.0, -1.0, 8.0, 15.0}
	                    };
      double X1[N],X2[N];
      double S,ERR,TOL;
      int I,J,NN,K,OK; 

      NN = 200; // Maximum number of iterations
      TOL = 1.0e-10; // Tol for checking convergence

      // Initial guess of the soln.
      for(I = 1; I <=N; I++)
          X1[I-1] = 0.0; 
	   	  
	  /* STEP 1 */
      K = 1;
      OK = false;
      /* STEP 2 */
      while ((!OK) && (K <= NN )) {
         /* err is used to test accuracy - it measures the  
            infinity-norm */
         ERR = 0.0; 
         /* STEP 3 */
         for (I=1; I<=N; I++) 
	 {
            S = 0.0;
            for (J=1; J<=I-1; J++) 
	    {
		S = S - A[I-1][J-1] * X2[J-1];
	    }
	    for(J = I+1; J <= N; J++)
	    {
		S = S - A[I-1][J-1] * X1[J-1];
	    }
            X2[I-1] = (S + A[I-1][N]) / A[I-1][I-1];

            //if (absval(S) > ERR) ERR  = absval(S);
            //X1[I-1] = X1[I-1] + S;
         }

	 // find L_{infinity} error
	 S = absval(X2[0] - X1[0]);
	 ERR = absval(X2[0]);
	 for (I=1; I<=N; I++) 
	 {
	     if(absval(X2[I-1] - X1[I-1]) > S)
		 S = absval(X2[I-1] - X1[I-1]);
	     if(absval(X2[I-1]) > ERR)
		 ERR = absval(X2[I-1]);
	 }

         /* STEP 4 */
	 // relative error
	 ERR = S/ERR;

         if (ERR <= TOL) OK = true;
         /* process is complete */
         else
	 {
            /* STEP 5 */
            K++;
            /* STEP 6 */
	    for(I=1; I<=N; I++) X1[I-1] = X2[I-1];
	 }
      }

      if (!OK) cout<<"Maximum Number of Iterations Exceeded."<<endl;
      /* STEP 7 */
      /* procedure completed unsuccessfully */
      else
      {
	  cout<<"Number of iteration = "<< K <<endl;
	  cout<<"The solution vector is:"<<endl;
	  for (I=1; I<=N; I++)
	      cout<<"X["<<I <<"]= " << X2[I-1] <<endl;
      }

      return 0;
}

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