#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 SparseMatrix
{
public:
   int nnz;
   float val[6];
   int row[6];
   int col[6];

   Vector operator*(Vector &v)
   {
      Vector z;
      int k;

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

      for (k = 0; k < nnz; k = k + 1)
         z.x[row[k]] = z.x[row[k]] + val[k]*v.x[col[k]];     

      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(SparseMatrix & P, SparseMatrix & 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[])
{
   SparseMatrix P, Q;
   Vector sol, b;

   P.nnz = 6;
   P.val[0] = -0.1; P.val[1] = -0.2; P.val[2] = -0.1;
   P.val[3] = -0.3; P.val[4] = -0.2; P.val[5] = -0.3;
   P.row[0] = 0; P.row[1] = 0; P.row[2] = 1;
   P.row[3] = 1; P.row[4] = 2; P.row[5] = 2;
   P.col[0] = 1; P.col[1] = 2; P.col[2] = 0;
   P.col[3] = 2; P.col[4] = 0; P.col[5] = 1;

   // 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.nnz = 3;
   Q.val[0] = 0.1; Q.val[1] = 0.1; Q.val[2] = 0.1;
   Q.row[0] = 0; Q.row[1] = 1; Q.row[2] = 2;
   Q.col[0] = 0; Q.col[1] = 1; Q.col[2] = 2;

   // 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" ;
}