// Sets up MTL-style matrices as operators
#ifndef MTL_OP_H
#define MTL_OP_H

#include "operator.h"

template<class Vector>
void zero(Vector &z)
{
  for ( int i = 0; i < z.size(); ++i )
    z[i] = 0;
}

template<class Vector, class Matrix>
class MatrixOp : public Op<Vector> {
 private:
  const Matrix *A;
 public:
  MatrixOp(const Matrix &_A) : A(&_A)  {  }

  Vector &operator()(const Vector &x, Vector &y)
  { zero(y);  mult(*A,x,y);  return y;  }
};


template <class Vector, class Matrix, class PermType>
class LUSolveOp : public Op<Vector> {
 private:
  Matrix &LU;
  PermType   &pivot;
 public:
  LUSolveOp(Matrix &_LU, PermType &_pivot) : LU(_LU), pivot(_pivot)
  {  }
  Vector &operator()(const Vector &b, Vector &x)
  {
    zero(x);
    lu_solve(LU, pivot, b, x);
    return x;
  }
};

template <class Vector, class Matrix>
class GaussSeidelOp : public Op<Vector> {
 private:
  const typename mtl::rows_type<Matrix>::type A;
 public:
  GaussSeidelOp(const Matrix &_A) : A(rows(_A))  {  }

  typedef typename Vector::value_type value_type;
  Vector &operator()(const Vector &b, Vector &x)
  {
    for ( int i = 0; i < A.nrows(); ++i )
      {
        typename Matrix::Row r = A[i];
        typename Matrix::Row::iterator j;
        value_type sum = b[i], diag = 0;
        for ( j = r.begin(); j != r.end(); ++j )
          {
            if ( j.column() == i ) diag = *j;
            else  sum -= (*j)*x[j.column()];
          }
        x[i] = sum / diag;
      }
    return x;
  }
}; // GaussSeidelOp


template <class Vector, class Matrix>
class GaussSeidelRevOp : public Op<Vector> {
 public:
  typedef typename Vector::value_type value_type;
 private:
  const typename mtl::rows_type<Matrix>::type A;
 public:
  GaussSeidelRevOp(const Matrix &_A) : A(rows(_A))
  {  }

  Vector &operator()(const Vector &b, Vector &x)
  {
    int i;
    // Gauss Seidel in reverse order
    for ( i = A.nrows()-1; i >= 0; --i )
      {
        typename Matrix::Row r = A[i];
        typename Matrix::Row::iterator j;
        value_type sum = b[i], diag = 0;
        for ( j = r.begin(); j != r.end(); ++j )
          {
            if ( j.column() == i )
              diag = *j;
            else
              sum -= (*j)*x[j.column()];
          }
        x[i] = sum / diag;
      }
    return x;
  }
}; // GaussSeidelRevOp



#endif // MTL_OP_H