#include <iostream>
#include <cstdlib>
#include <mtl/mtl.h>
#include "mtl-op.h"

using std::ostream;

typedef double value_type;
using namespace mtl;

typedef dense1D<value_type> Vector;
typedef matrix<value_type, rectangle<>, compressed<>, row_major>::type Matrix;

typedef MatrixOp<Vector, Matrix> MatOp;
typedef GaussSeidelOp<Vector, Matrix> GSOp;
typedef GaussSeidelRevOp<Vector, Matrix> GSROp;

ostream &operator<<(ostream &s, Vector &x)
{
  s << "Vector: dim: " << x.size() << std::endl;
  for ( int i = 0; i < x.size(); ++i )
    {
      s << x[i] << "  ";
      if ( i % 6 == 5 )
        s << std::endl;
    }
  s << std::endl;
  return s;
}


int  main(int argc, char *argv[])
{
  int m = 10, n = 10;
  Matrix A(m,n);

  for ( int i = 0; i < n-1; ++i )
    {
      A(i,i) = 2;
      A(i,i+1) = -1;
      A(i+1,i) = -1;
    }
  A(n-1,n-1) = 2;

  print_all_matrix(A);

  Vector x(n), b(n), y(n);

  for ( int i = 0; i < b.size(); ++i )
    b[i] = rand() / value_type(RAND_MAX);
  x = 0;
  x.resize(n);

  MatOp Aop(A);
  GSOp  gs(A);
  GSROp gsr(A);

  std::cout << "b = " << b;
  std::cout << "x = " << x;
  Aop(b,y);
  std::cout << "y = A*b = " << y;
  gs(b,x);
  std::cout << "Result of Guass-Seidel = " << x;
  gsr(b,x);
  std::cout << "... followed by reverse Gauss-Seidel = " << x;
}