/*
*   JACOBI ITERATIVE ALGORITHM 7.1
*
*   To solve Ax = b given an initial approximation x(0).
*   Example 7.3.1
*
*   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 NN.
*
*    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 A[N][N+1] = {{4.0, 1.11, -1.0, 1.0, 1.12},
                           {1.2, 4.122, -1.0, -1.0, 3.44},
                           {0.0, 1.0, 99.9, 1.0, 2.15},
                           {1.3, -0.11, -1.1, 10.0, 4.16}};

      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<=N; J++) 
	    {
	        if(J == I) continue;
		S = S - A[I-1][J-1] * X1[J-1];
	    }
	    X2[I-1] = (S + A[I-1][N]) / A[I-1][I-1];
         }

	 // find relative error in L_{infinity} norm
	 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;
}