// $Id: mg.h,v 1.1.1.1 2000/11/02 08:47:16 fritzi Exp $

// This file is part of the DEAL Library
// DEAL is Copyright(1995) by
// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier


#ifndef __mg_h
#define __mg_h

#ifndef __iteration_h
#include <deal/iteration.h>
#endif

#ifndef __math_h
#include <math.h>
#endif

/*
CLASS
   MGInfo
   */
class MGInfo : public RelativeIterationInfo
{
  protected:
  int l0,l1;
  double res0;

  public:
  //////////
  MGInfo(double residual=1.e-6, int pstep=1, int maxit=10,
	 const char* txt=std_text):
    RelativeIterationInfo(0.,residual,pstep,maxit,txt) { l0=-1; l1=-1; }
  //////////
  int check_residual(int i, double res)
  {
    usediterations = i;
    if (i==1) 
	res0 = res; 
    return RelativeIterationInfo::check_residual(i,res);
  }
  //////////
  double& res() {return aimedresidual;}
  //////////
  int& maxiter() {return maxiterations;}
  //////////
  void minlevel(int i) {l0=i;}
  //////////
  int minlevel() {return l0;}
  //////////
  void maxlevel(int i) {l1=i;}
  //////////
  int maxlevel() {return l1;}
  //////////
  double rho() { return pow(residual()/res0,1./(usediter()-1)); }
};

/**********************************************************/

/*
CLASS
   RelativeMGInfo
   */
class RelativeMGInfo : public MGInfo
{
  protected:
  int l0,l1;
  double factor,global_aimed;
  double res0;

public:
  //////////
  RelativeMGInfo(double f=1e-3, double residual=1.E-6, 
		 int pstep=1, int maxit=10,
	 const char* txt=std_text):
    MGInfo(residual,pstep,maxit,txt) 
  { factor=f; l0=-1; l1=-1; global_aimed = residual;}
  //////////
  void minlevel(int i) {l0=i;}
  //////////
  int minlevel() {return l0;}
  //////////
  void maxlevel(int i) {l1=i;}
  //////////
  int maxlevel() {return l1;}
  //////////
  void relres(double d) {factor=d;}
  //////////
  double& reduction() {return factor;}
  double& globres  () {return global_aimed;}
  void    globres  (double d) {global_aimed = d;}
  //////////
  int check_residual(int i, double res)
  {
    usediterations = i;
    if (i==1)
    {
      res0 = res; 
      double newres = res*factor;
      aimedresidual = (newres > global_aimed) ? newres : global_aimed;
    }
    return IterationInfo::check_residual(i,res);
  }
  //////////
  double rho() { return pow(residual()/res0,1./(usediter()-1)); }
};

////////////////////////////////////////////////////////////
template<class MATRIX, class VECTOR, class MEM>
inline double mgstep( int l, MATRIX& A, 
	     VECTOR& u, VECTOR& f, MEM& mem, MGInfo& info, double& res)
{
  VECTOR& w = mem.v(0);
  VECTOR& d = mem.v(1);

  if(l==info.minlevel())
  {
    res += A.smooth(l,u,f,d,'e');
  }
  else
  {
    res += A.smooth(l,u,f,d,'v');
    d.v(l-1) = 0.;
    d.v(l) = 0.;
    A.mgresidual(l, d, u, f);
    A.restrict(l,d.v(l-1),d.v(l));
    A.mgmult (l-1, d, w);
    f.v(l-1) = d.v(l-1);
    mgstep(l-1, A, u, f, mem, info, res);
    d.v(l-1) = u.v(l-1);
    d.v(l-1).add(-1., w.v(l-1));
    d.v(l) = 0.;
    A.prolongate(l,d.v(l),d.v(l-1));
    u.v(l).add(d.v(l));
    res += A.smooth(l,u,f,d,'n');
  }
  return res;
}

////////////////////////////////////////////////////////////
template<class MATRIX, class VECTOR, class MEM>
inline void mg( MATRIX& A, VECTOR& u, VECTOR& f, MEM& mem, MGInfo& info)
{
  int i,reached=0;
  double res=1.;
  VECTOR& w = mem.v(0);

  if(info.minlevel()==-1)
  {
    info.minlevel(0);
  }
  if(info.maxlevel()==-1)
  {
    info.maxlevel(A.maxlevel());
  }

  for(i=1;!reached;i++)
  {
    res = 0.;
    w = u;
    res = mgstep(info.maxlevel(),A,u,f,mem,info,res);
    reached = info.check_residual(i,res);
  }
  info.minlevel(-1);
  info.maxlevel(-1);
}

////////////////////////////////////////////////////////////
template<class MATRIX, class VECTOR, class MEM>
inline void mg_nores( MATRIX& A, VECTOR& u, VECTOR& f, MEM& mem, MGInfo& info)
{
  int i,reached=0;
  double res=1.;
  VECTOR& w = mem.v(0);

  if(info.minlevel()==-1)
  {
    info.minlevel(0);
  }
  if(info.maxlevel()==-1)
  {
    info.maxlevel(A.maxlevel());
  }

  for(i=0;!reached;i++)
  {
    reached = info.check_residual(i,1.e300);
    w = u;
    res = mgstep(info.maxlevel(),A,u,f,mem,info);
  }
  info.minlevel(-1);
  info.maxlevel(-1);
}
#endif