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

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

#ifndef __qmr_h
#define __qmr_h
#include <math.h>
#include <stdio.h>

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

////////// 
//
// QMR-Verfahren
//
template<class MATRIX,class VECTOR,class MEM>
inline int qmr(MATRIX& A, VECTOR& x, const VECTOR& b,
		    MEM& mem, IterationInfo& info)
{
  int it, reached;
  double res;
  VECTOR&r = mem.v(0);
  VECTOR&y = mem.v(2);
  VECTOR&wtil = mem.v(3);
  VECTOR&z = mem.v(4);
  VECTOR&v = mem.v(5);
  VECTOR&w = mem.v(6);
  VECTOR&ytil = mem.v(7);
  VECTOR&p = mem.v(8);
  VECTOR&ptil = mem.v(9);
  VECTOR&q = mem.v(10);
  VECTOR&qtil = mem.v(11);
  VECTOR&d = mem.v(12);

  double rho, rho2, xi, gamma, gamma2, beta, eta, eps, theta, theta2, delta;

  res = A.residual(r,x,b);
  reached = info.check_residual(0,res);
  y = r;
  rho = sqrt(y*y);
  wtil = r;

  A.precondition(z,wtil);
  xi = sqrt(z*z);
  gamma = 1.;
  eta = -1.;
  theta = 0.;

  for (it=1;!reached;it++)
    {
      if ((fabs(rho)<1.e-19) || (fabs(xi)<1.e-19)) return -1;
      y.equ(1./rho,y);
      v = y;
      w.equ(1./xi,wtil);
      z.equ(1./xi,z);

      delta = z*y;
      if (fabs(delta)<1.e-19) return -1;
      A.precondition(ytil,y);
      if (it==1)
	{
	  p = ytil;
	  q = z;
	}
      else
	{
	  p.sadd(-xi*delta/eps,1.,ytil);
	  q.sadd(-rho*delta/eps,1.,z);
	}
      A.vmult(ptil,p);
      eps = q*ptil;
      if (fabs(eps)<1.e-19) return -1;
      beta = eps/delta;
      if (fabs(beta)<1.e-19) return -1;

      y.equ(1.,ptil,-beta,v);
      rho2 = sqrt(y*y);
      A.Tvmult(wtil,q);
      wtil.add(-beta,w);
      A.Tprecondition(z,wtil);
      xi = sqrt(z*z);
      theta2 = rho2/(gamma*fabs(beta));
      gamma2 = 1./sqrt(1+theta2*theta2);

      eta *= -rho*gamma2*gamma2/(beta*gamma*gamma);

      if (it==1)
	{
	  d.equ(eta,p);
	}
      else
	{
	  double weissnich = theta*theta * gamma2*gamma2;
	  d.sadd(weissnich,eta,p);
	}
      x.add(d);

      res = A.residual(r,x,b);
      reached = info.check_residual(it,res);

      gamma = gamma2;
      rho = rho2;
      theta = theta2;
    }
  if (reached<0) return 1;
  return 0;
}
#endif