typedef valarray<double> Vector;

class MDSPreconditioner
{
 public:

  MDSPreconditioner (int lvls)
    : levels(lvls), jacobi(lvls) // (reserve lvls entries in jacobi)
  { ; }

  void setJacobi (int l, const Vector& diag)
  { /* set jacobi[l] to diag; */ }

  void solve (const Vector& r, Vector& w)
  { /* call the routine MDS (...) on finest level */ }

 private:
  int levels;
  vector<JacobiPreconditioner> jacobi;
  
  void MDS (int l, const Vector& r, Vector& w);
    // (this is essentially the routine from Tutorial 7)
  {
    // jacobi[l].solve (r, w);
    if (l != 0)
    {
       Restrict (r, r_coarse);
       MDS (l-1, r_coarse, w_coarse);
       Prolongate (w_coarse, w_fine);
       w += w_fine;
    }
  }
}; // class MDSPreconditioner

void main ()
{
  ...
  MDSPreconditioner M(levels);
  // create coarsest mesh

  for (l=0; l<levels; l++)
  {
    // assemble stiffness matrix K, load vector f, implement B.C.
    // M.setJacobi (l, ...)
    // if (l != levels-1) refineUniform (...)
  }

  CG (K ... f ... M ...);
  ...
} // main