// $Id: cg.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 __cg_h
#define __cg_h
#ifndef __iteration_h
#include <deal/iteration.h>
#endif
#include <math.h>

//////////
// GROUP: Iterative Linear Solvers
template<class MT, class VT, class MEM>
inline int cg_simple(MT& A, VT& x, VT& b, MEM& mem, IterationInfo& info)
{
  int it, reached=0;
  double res, alpha,beta;
  VT& d =mem.v(0);
  VT& g =mem.v(1);
  VT& Ad=mem.v(2);

  res = A.residual(d,x,b);
  reached = info.check_residual(0,res);
  g.equ(-1.,d);
  beta = g*g;

  for (it=0;!reached;it++)
  {
    A.vmult(Ad,d);
    alpha = Ad*d;
    alpha = beta/alpha;

    x.add(alpha,d);
    g.add(alpha,Ad);

    res = g*g;
    reached = info.check_residual(it,sqrt(res));

    if (reached) continue;
    d.sadd(res/beta,-1.,g);
    beta = res;
  }
  if (reached<0) return 1;
  return 0;
}

//////////
template <class MT,class VT,class MEM>
inline int cg(MT& M,VT& x,VT& b,MEM& mem, IterationInfo& info)
{
  VT& g = mem.v(0);
  VT& h = mem.v(1);
  VT& d = mem.v(2);
  VT& Ad = mem.v(3);

  int  it,reached = 0;
  double res,gh,alpha,beta;

  res = M.residual(g,x,b);

  g.equ(-1.,g);
  reached = info.check_residual(0,sqrt(g*g));

  M.precondition(h,g);

  d.equ(-1.,h);

  gh = g*h;

  for(it=0;!reached;it++)
  {
    M.vmult(Ad,d);

    alpha = d*Ad;
    alpha = gh/alpha;

    g.add(alpha,Ad);
    x.add(alpha,d );
    res = sqrt(g*g);

    reached = info.check_residual(it,res);
    if (reached) continue;
    
    M.precondition(h,g);
    
    beta = gh;
    gh   = g*h;
    beta = gh/beta;
    
    d.sadd(beta,-1.,h);
  }
  if (reached<0) return 1;
  return 0;
}
#endif