// $Id: bicgstab.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 __bicgstab_h
#define __bicgstab_h
#include <math.h>
#include <stdio.h>

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

////////// 
//
// Bi-cgstab Verfahren
//
template<class MATRIX,class VECTOR,class MEM>
inline int bicgstab(MATRIX& A, VECTOR& x,  VECTOR& b,
		    MEM& mem, IterationInfo& info)
{
  int it, reached;
  double res;
  VECTOR& r = mem.v(0);
  VECTOR& rbar = mem.v(1);
  VECTOR& p = mem.v(2);
  VECTOR& y = mem.v(3);
  VECTOR& z = mem.v(4);
  VECTOR& s = mem.v(5);
  VECTOR& t = mem.v(6);
  VECTOR& v = mem.v(7);
  double alpha,beta,omega,rho,rhorho, sprod;

  res = A.residual(r,x,b);
  reached = info.check_residual(0,res);
  rbar.equ(1.,r);
  v.reinit(v);
  p.reinit(v);
  alpha=omega=rho=1;
  for (it=1;!reached;it++)
  {
    rhorho = rbar*r;
    beta = rhorho*alpha/(rho*omega);
    rho = rhorho;
    p.sadd(beta, 1.,r, -beta*omega,v);
    A.precondition(y,p);
    A.vmult(v,y);
    rhorho = rbar*v; 

    if (fabs(rhorho)<1.e-19) return -1;

    alpha = rho/rhorho;
    s.equ(1., r, -alpha, v);
    A.precondition(z,s);
    A.vmult(t,z);
    rhorho = t*s;
    sprod = t*t;
    omega = rhorho/sprod;
    x.add(alpha, y, omega, z);
    r.equ(1.,s,-omega,t);
    res = sqrt(r*r);
    reached = info.check_residual(it,res);
  }
  if (reached<0) return 1;
  return 0;
}

////////// 
//
// Bi-cgstab Verfahren mit exakter Residuumsberechnung
//
template<class MATRIX,class VECTOR,class MEM>
inline int bicgstab_ex(MATRIX& A, VECTOR& x, const VECTOR& b,
		    MEM& mem, IterationInfo& info)
{
  int it, reached;
  double res;
  VECTOR& r = mem.v(0);
  VECTOR& rbar = mem.v(1);
  VECTOR& p = mem.v(2);
  VECTOR& y = mem.v(3);
  VECTOR& z = mem.v(4);
  VECTOR& s = mem.v(5);
  VECTOR& t = mem.v(6);
  VECTOR& v = mem.v(7);
  double alpha,beta,omega,rho,rhorho, sprod;

  res = A.residual(r,x,b);
  reached = info.check_residual(0,res);
  rbar.equ(1.,r);
  v.reinit(v);
  p.reinit(v);
  alpha=omega=rho=1;
  for (it=1;!reached;it++)
  {
    rhorho = rbar*r;
    beta = rhorho*alpha/(rho*omega);
    rho = rhorho;
    p.sadd(beta, 1.,r, -beta*omega,v);
    A.precondition(y,p);
    A.vmult(v,y);
    rhorho = rbar*v; 

    if (fabs(rhorho)<1.e-19) return -1;

    alpha = rho/rhorho;
    s.equ(1., r, -alpha, v);
    A.precondition(z,s);
    A.vmult(t,z);
    rhorho = t*s;
    sprod = t*t;
    omega = rhorho/sprod;
    x.add(alpha, y, omega, z);
    r.equ(1.,s,-omega,t);
    res = A.residual(t,x,b);
    reached = info.check_residual(it,res);
  }
  if (reached<0) return 1;
  return 0;
}
    
#endif