#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"

#include <itl/interface/mtl.h>
#include "itl/preconditioner/ssor.h"
#include "itl/krylov/qmr.h"


/*
  In thsi example, we show how to use SSOR, the output should be:

iteration 0: resid 2.23607
iteration 1: resid 5.65543
iteration 2: resid 1.04817
iteration 3: resid 0.136963
iteration 4: resid 1.32456e-13
finished! error code = 0
4 iterations
1.32456e-13 is actual final residual.
5.92363e-14 is actual relative tolerance achieved.
Relative tol: 1e-06  Absolute tol: 0
     1      2      0      0      3   x          0.205847  =       1
     4      5      0      6      0   x         0.0419363  =       1
     0      7      8      0      9   x         -0.178049  =       1
     0      0     10     11     12   x       -0.00551162  =       1
     0      0     13      0     14   x           0.23676  =       1
*/

using namespace mtl;
using namespace itl;

typedef  double Type;

//begin
typedef matrix< Type, rectangle<>, 
	        array< compressed<> >, 
                row_major >::type Matrix;

//end

int main (int , char* []) 
{
  using std::cout;
  using std::endl;

  int max_iter = 50;
  //begin
  Matrix A(5, 5);
  //end
  A(0, 0) = 1.0;
  A(0, 1) = 2.0;
  A(0, 4) = 3.0;
  A(1, 0) = 4.0;
  A(1, 1) = 5.0;
  A(1, 3) = 6.0;
  A(2, 1) = 7.0;
  A(2, 2) = 8.0;
  A(2, 4) = 9.0;
  A(3, 2) = 10.;
  A(3, 3) = 11.;
  A(3, 4) = 12.;
  A(4, 2) = 13.; 
  A(4, 4) = 14.;

  //begin
  dense1D<Type> x(A.nrows(), Type(0));
  dense1D<Type> b(A.ncols());
  for (dense1D<Type>::iterator i=b.begin(); i!=b.end(); i++)
    *i = 1.;

  //SSOR preconditioner
  SSOR<Matrix> precond(A);
  
  noisy_iteration<double> iter(b, max_iter, 1e-6);
  qmr(A, x, b, precond.left(), precond.right(), iter);
  //end

  //verify the result
  dense1D<Type> b1(A.ncols());
  itl::mult(A, x, b1);
  itl::add(b1, itl::scaled(b, -1.), b1);

  if (itl::two_norm(b1) < 1.e-6) { //output
    for (int i=0; i<5; i++) {
      for (int j=0; j<5; j++) {
	cout.width(6);
	cout << A(i, j) << " ";
      }
      cout << "  x  ";
      cout.width(16);
      cout << x[i] << "  =  ";
      cout.width(6);
      cout << b[i] << endl;
    }
  } else {
    cout << "Residual " << iter.resid() << endl;
  }

  return 0;
  }