#ifndef OP_MG_H
#define OP_MG_H

// Multigrid using arrays of operators
#include "my-mg.h"
#include "operator.h"

template<typename Types>
class OpMG : public MY_MG::multigrid<Types> {
 public:
  // typedef's from Types
  typedef typename Types::value_type value_type; // Scalar type
  typedef typename Types::int_type   int_type;   // Integer type
  typedef typename Types::size_type  size_type;  // Matrix size type
  typedef typename Types::Vector     Vector;     // Vector type
  typedef typename Types::Grid       Grid;       // Support vectors
  
  typedef Op<Vector> *PtrOp;
 private:
  // These are the arrays of (pointers to) operators that 
  // define the multigrid preconditioner.
  // We use pointers to Op's to allow polymorphism
  // via virtual functions
  PtrOp *A, *pre_rel, *post_rel, *interp, *restrikt, top_solve;
 public:
  // main constructor
  OpMG(int_type n_levels, int_type _dims[], PtrOp _A[], 
       PtrOp _pre_rel[], PtrOp _post_rel[], PtrOp _interp[], PtrOp _restrikt[],
       PtrOp _top_solve) : MY_MG::multigrid<Types>(n_levels)
  {
    int_type _nu[] = { 1, 1 };
    resize(_dims);
    set_nu(_nu);
    A        = new PtrOp[n_levels];
    pre_rel  = new PtrOp[n_levels];
    post_rel = new PtrOp[n_levels];
    interp   = new PtrOp[n_levels-1];
    restrikt = new PtrOp[n_levels-1];
    for ( int i = 0; i < n_levels-1; ++i )
      {
        A[i]        = _A[i];
        pre_rel[i]  = _pre_rel[i];
        post_rel[i] = _post_rel[i];
        interp[i]   = _interp[i];
        restrikt[i] = _restrikt[i];
      }
    A[n_levels-1] = _A[n_levels-1];
    top_solve = _top_solve;
  }
  
  // Destructor
  ~OpMG()
  {
    delete[] A;
    delete[] pre_rel;
    delete[] post_rel;
    delete[] interp;
    delete[] restrikt;
  }
  
  // now here we put the essential (virtual) functions
  // defining the multigrid preconditioner
  Vector &pre_relax (Vector &v, const Vector &b, int_type level)
  {  return (*pre_rel[level])(b, v);  }
  Vector &post_relax(Vector &v, const Vector &b, int_type level)
  {  return (*post_rel[level])(b, v);  }
  Vector &solve     (Vector &v, const Vector &b, int_type level)
  {  return (*top_solve)(b, v);  }
  // Note that "out" is on level "level", while "v" is on "level+1"
  Vector &interpolate(const Vector &v, int_type level, Vector &out)
  {  return (*interp[level])(v, out);  }
  // Note that "out" is on level "level", while "v" is on "level-1"
  Vector &restriction(const Vector &v, int_type level, Vector &out)
  {  return (*restrikt[level-1])(v, out);  }
  Vector &apply_operator(const Vector &v, int_type level, Vector &out)
  {  return (*A[level])(v, out);  }
};

#endif // OP_MG_H