#ifndef __MDS_H
#define __MDS_H

#include <iostream>
#include <vector>
#include "vector.hh"     // depends on which "Vector" you use
#include "smatrix.hh"    // depends on which sparse matrix you use

using namespace std;

void Restrict (const Vector& fineRes, Vector& coarseRes);

void Prolongate (const Vector& coarseVec, Vector& fineVec);


class MDSPreconditioner
{
public:
  // constructor: initialize with (total) number of levels
  MDSPreconditioner (int numLevels);

  // destructor
  ~MDSPreconditioner ();
  
  // set the diagonal for Jacobi preconditioner for a certain level
  void InitDiagonal (int level, const SMatrix& K);
  
  // apply the preconditioner C^{-1}
  inline
  void Solve (const Vector& r, Vector& w) const
  {
    ApplyCL (numLevels_-1, r, w);
  }

private:

  // recursive helper routine
  void ApplyCL (int level, const Vector& r, Vector& w) const;

  //--- members

  int numLevels_;                // number of levels
  std::vector<Vector> jacobi_;   // diagonal for Jacobi preconditioner on each level
  
}; // class MDSPreconditioner

#endif