// $Id: cr.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 __cr_h
#define __cr_h

#include <math.h>

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

//////////
// GROUP: Iterative Linear Solvers
template <class MT,class VT,class MEM>
void cr(MT& M,VT& x,VT& b,MEM& mem,IterationInfo& info)
{
  int l,reached=0;
  double f,alfa0=0.,alfa1,lambda,om1,om0;
  VT& h  =mem.v(0);
  VT& p0 =mem.v(1);
  VT& p1 =mem.v(2);
  VT& wa0=mem.v(3);
  VT& wa1=mem.v(4);
  VT& a0 =mem.v(5);
  VT& a1 =mem.v(6);
  VT& aa =mem.v(7);
  VT& r  =mem.v(8);

  M.residual(r,x,b);
  M.precondition(p1,r);

  M.vmult(a1,p1);

  for (l=0;!reached;l++)
  {
    M.precondition(wa1,a1);  
    M.vmult(aa,wa1);

    om1 = a1*wa1;
    if (om1)
    { 
      lambda = (r*wa1)/om1;
      alfa1  = (aa*wa1)/om1;
    }
    else 
    {
      lambda  = 0.;
    }

    if (l) alfa0 = (aa*wa0)/om0; 

    x.sadd( 1., lambda,p1);
    r.sadd( 1.,-lambda,a1);
    f = sqrt(r*r);
    reached = info.check_residual(l,f);    
    if (!lambda) reached = 1;

    h.equ(1.,p1);
    p1.sadd(-alfa1,1.,wa1);
    p1.sadd(1.,-alfa0,p0);
    p0.equ(1.,h);

    h.equ(1.,a1);
    a1.sadd(-alfa1,1.,aa);
    a1.sadd(1.,-alfa0,a0);
    a0.equ(1.,h);

    wa0 = wa1;

    om0 = om1;
  }
}
#endif