#include <iostream.h>

class Vector
{
public:
   double 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:
   double A[3][3];

   virtual 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;
   }
};


class SparseMatrix : public DenseMatrix
{
public:
   int nnz;
   double val[6];
   int row[6];
   int col[6];

   virtual 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;
   }
};





double distance(Vector & a, Vector & b)
{
   double 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;
   double 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 P1, Q1;
   Vector sol, b;

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

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

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

   sol = Jacobi(&P1, &Q1, b);

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

   cout << "\n";

   SparseMatrix P2, Q2;

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

   Q2.nnz = 3;
   Q2.val[0] = 0.1; Q2.val[1] = 0.1; Q2.val[2] = 0.1;
   Q2.row[0] = 0; Q2.row[1] = 1; Q2.row[2] = 2;
   Q2.col[0] = 0; Q2.col[1] = 1; Q2.col[2] = 2;

   sol = Jacobi(&P2, &Q2, b);

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

}