#include <iostream.h>

class Vector
{
public:
   float x[3];

   Vector operator+(Vector &v)
   {
      Vector z;
      int i;
      for (i = 0; i < 3; i = i + 1)
         z.x[i] = x[i] + v.x[i];
      return(z);
   }
};

class DenseMatrix
{
public:
   float A[3][3];

   Vector operator*(Vector &v)
   {
      Vector z;
      int i, j;

      for (i = 0; i < 3; i=i+1)
         z.x[i] = 0;

      for (i = 0; i < 3; i=i+1)
         for (j = 0; j < 3; j=j+1)
            z.x[i] = z.x[i] + A[i][j]*v.x[j];     

      return z;
   }
};

float distance(Vector & a, Vector & b)
{
   float sum = 0;
   int i;

   for (i = 0; i < 3; i = i + 1)
      sum = sum + fabs(a.x[i] - b.x[i]);

   return(sum);
}

Vector Jacobi(DenseMatrix & P, DenseMatrix & Q, Vector b)
{
   Vector v1, v2;
   int i;
   float d;

   for (i = 0; i < 3; i = i + 1)
   {  v1.x[i] = 0;
      v2.x[i] = 1;
   }
                  

   d = 1.0;
   while ( d > 0.0001 )
   {
      cout << v1.x[0] << "\t" << v1.x[1] << "\t" << v1.x[1] << "\n" ;
      v2 = P*v1 + Q*b;
      d = distance(v1, v2);
      v1 = v2;
   }
   return(v2);
}

//
// Solving: [ 10   1   2 ] [ x1 ]   [ 5 ]
//          [  1  10   3 ] [ x2 ] = [ 6 ]
//          [  2   3  10 ] [ x3 ]   [ 7 ]
//
// 
//  -1   [ 0.1  0    0  ]      -1    [
// M   = [  0  0.1   0  ]     M  N = [
//       [  0   0   0.1 ]            [
//
int main(int argc, char *argv[])
{
   DenseMatrix P, Q;
   Vector sol, b;

   P.A[0][0] = 0.0; P.A[0][1] = -0.1; P.A[0][2] = -0.2;
   P.A[1][0] = -0.1; P.A[1][1] = 0.0; P.A[1][2] = -0.3;
   P.A[2][0] = -0.2; P.A[2][1] = -0.3; P.A[2][2] = 0.0;

   Q.A[0][0] = 0.1; Q.A[0][1] = 0.0; Q.A[0][2] = 0.0;
   Q.A[1][0] = 0.0; Q.A[1][1] = 0.1; Q.A[1][2] = 0.0;
   Q.A[2][0] = 0.0; Q.A[2][1] = 0.0; Q.A[2][2] = 0.1;

   b.x[0] = 5; b.x[1] = 6; b.x[2] = 7;

   sol = Jacobi(P, Q, b);

   cout << sol.x[0] << "\t" << sol.x[1] << "\t" << sol.x[1] << "\n" ;
}