// $Id: gmres.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 __gmres_h
#define __gmres_h
#include <stdio.h>

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

#ifndef __dfmatrix_h
#include <deal/dfmatrix.h>
#endif

//////////
template<class MATRIX, class VECTOR, class MEM>
inline int gmres(MATRIX& A, VECTOR& x, VECTOR& b,
		      MEM& mem, IterationInfo& info)
{
  int kmax = mem.n();
  dFMatrix H(kmax+1, kmax), H1(kmax+1, kmax);
  dVector y(kmax), b0(kmax+1);
  int i,k,reached=0,iter;
  double rho,beta;
  
  double res = A.residual(mem.v(0),x,b);
  reached = info.check_residual(0,res);

  rho = sqrt(mem.v(0) * mem.v(0));
  beta = rho;
  
  mem.v(0).equ(1./rho, mem.v(0));
  
  for (k=0 ; k<kmax-1 && (!reached) ; k++)
  {
    A.vmult(mem.v(k+1), mem.v(k));
    
    H.reinit(k+2,k+1);
    if (k>0)  H.fill(H1);
    
    for (i=0 ; i<=k ; i++)
    {
      H(i,k) = mem.v(k+1) * mem.v(i);
      mem.v(k+1).add(-H(i,k),mem.v(i));
    }
    
    double s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    // Re-orthogonalization
    
    //printf("\n");
    for (i=0 ; i<=k ; i++)
    {
      double htmp = mem.v(k+1) * mem.v(i);
      //printf(" %e ",htmp);
      H(i,k) += htmp;
      mem.v(k+1).add(-htmp,mem.v(i));
    }
    //printf("\n");
    
    s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    mem.v(k+1).equ(1./s, mem.v(k+1));
    
    // Least - Squares
    
    y.reinit(k+1);
    b0.reinit(k+2);
    b0(0) = beta;
    H1 = H;
    rho = H.least_squares(y,b0);
    reached = info.check_residual(k,rho);
  }
  
  for (i=0 ; i<k ; i++)
      x.add(y(i), mem.v(i));
  
  if (reached<0) return 1;
  return 0;
}

//////////
template<class MATRIX, class VECTOR, class MEM>
inline int pgmres(MATRIX& A, VECTOR& x, VECTOR& b,
		      MEM& mem, IterationInfo& info)
{
  int kmax = mem.n()-1;
  dFMatrix H(kmax+1, kmax), H1(kmax+1, kmax);
  dVector y(kmax), b0(kmax+1);
  int i,k,reached=0,iter;
  double rho,beta;
  
  VECTOR& p = mem.v(kmax);
  
  double res = A.residual(p,x,b);
  reached = info.check_residual(0,res);
  A.precondition(mem.v(0),p);
  
  rho = sqrt(mem.v(0) * mem.v(0));
  beta = rho;
  
  mem.v(0).equ(1./rho, mem.v(0));
  
  for (k=0 ; k<kmax-1 && (!reached) ; k++)
  {
    A.vmult(p, mem.v(k));
    A.precondition(mem.v(k+1),p);
    
    H.reinit(k+2,k+1);
    if (k>0)  H.fill(H1);
    
    for (i=0 ; i<=k ; i++)
    {
      H(i,k) = mem.v(k+1) * mem.v(i);
      mem.v(k+1).add(-H(i,k),mem.v(i));
    }
    
    double s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    // Re-orthogonalization
    
    //printf("\n");
    for (i=0 ; i<=k ; i++)
    {
      double htmp = mem.v(k+1) * mem.v(i);
      //printf(" %e ",htmp);
      H(i,k) += htmp;
      mem.v(k+1).add(-htmp,mem.v(i));
    }
    //printf("\n");
    
    s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    mem.v(k+1).equ(1./s, mem.v(k+1));
    
    // Least - Squares
    
    y.reinit(k+1);
    b0.reinit(k+2);
    b0(0) = beta;
    H1 = H;
    rho = H.least_squares(y,b0);
    reached = info.check_residual(k,rho);
  }
  
  for (i=0 ; i<k ; i++)
      x.add(y(i), mem.v(i));
  
  if (reached<0) return 1;
  return 0;
}

//////////
template<class MATRIX, class VECTOR, class MEM>
inline int pgmres0(MATRIX& A, VECTOR& x, VECTOR& b,
		      MEM& mem, IterationInfo& info)
{
  int kmax = mem.n()-1;
  dFMatrix H(kmax+1, kmax), H1(kmax+1, kmax);
  dVector y(kmax), b0(kmax+1);
  int i,k,reached=0,iter;
  double rho,beta;
  
  VECTOR& p = mem.v(kmax);
  
  double res = A.residual0(p,x,b);
  reached = info.check_residual(0,res);
  A.precondition0(mem.v(0),p);
  
  rho = sqrt(mem.v(0) * mem.v(0));
  beta = rho;
  
  mem.v(0).equ(1./rho, mem.v(0));
  
  for (k=0 ; k<kmax-1 && (!reached) ; k++)
  {
    A.vmult0(p, mem.v(k));
    A.precondition0(mem.v(k+1),p);
    
    H.reinit(k+2,k+1);
    if (k>0)  H.fill(H1);
    
    for (i=0 ; i<=k ; i++)
    {
      H(i,k) = mem.v(k+1) * mem.v(i);
      mem.v(k+1).add(-H(i,k),mem.v(i));
    }
    
    double s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    // Re-orthogonalization
    
    //printf("\n");
    for (i=0 ; i<=k ; i++)
    {
      double htmp = mem.v(k+1) * mem.v(i);
      //printf(" %e ",htmp);
      H(i,k) += htmp;
      mem.v(k+1).add(-htmp,mem.v(i));
    }
    //printf("\n");
    
    s = sqrt(mem.v(k+1) * mem.v(k+1));
    H(k+1,k) = s;
    
    mem.v(k+1).equ(1./s, mem.v(k+1));
    
    // Least - Squares
    
    y.reinit(k+1);
    b0.reinit(k+2);
    b0(0) = beta;
    H1 = H;
    rho = H.least_squares(y,b0);
    reached = info.check_residual(k,rho);
  }
  
  for (i=0 ; i<k ; i++)
      x.add(y(i), mem.v(i));
  
  if (reached<0) return 1;
  return 0;
}

#endif