/**
================================================================================
=                       Class GmresSolver header file                         =
================================================================================
Author : C.Chambeyron (24 november 2004)
Version: 1
================================================================================
 
  This class implement the GMRES method, preconditionned or not.
The preconditionner is a default argument.
  GMRES is suitable for non symetric systems
 
 For more explaination about BiConjugate Gradient Stabilized method,
see http://mathworld.wolfram.com/GMRES.html.
 
Base class : IterativeSolver
 
================================================================================
**/


#ifndef GMRESSOLVER_H
#define GMRESSOLVER_H

#include "IterativeSolver.h++"

class GmresSolver : public IterativeSolver
{
public:

  int_t DimKrylov;       //dimension of Krylov space

/// Constructor:
GmresSolver(){IterativeSolver();DimKrylov=100;};  
GmresSolver(int_t M, real_t eps=1.e-6,int_t d=100)
        {MaxIter=M; epsilon=eps;DimKrylov=d;}
GmresSolver( real_t eps,int_t M,int_t d=100)
        {MaxIter=M; epsilon=eps;DimKrylov=d;}  
   

template < class O, class Vec  >
  Vec   operator()(O &A,Vec &B)
  { 
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver GMRES in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     Vec X0(B);
     X0.setEntries(0.);
     preconditioning=false;
     Preconditioner<O> NoPrec;
    return  (*this)(A,B,X0,NoPrec);
  }

template < class O, class Vec  >
  Vec   operator()(O &A,Vec &B, Vec &X0)
  {
 #ifdef ME_DEBUG_ON
    me_TRACE->log(" Solver GMRES in  ",MaxIter," iterations max and tolerance is ",epsilon);
 #endif
 preconditioning=false;
 Preconditioner<O> NoPrec;
 return  (*this)(A,B,X0,NoPrec);
  }

  
template < class O, class Vec, class Prec >
  Vec   operator()(O &A,Vec &B,Prec &PC)
  { 
    #ifdef ME_DEBUG_ON
       me_TRACE->log(" Solver GMRES in  ",MaxIter," iterations max and tolerance is ",epsilon);
    #endif
    Vec X0(B);
    X0.setEntries(0.);
    return  (*this)(A,B,X0,PC);
  }

template < class O,  class Vec,  class Prec >
  Vec   operator()(O &A,Vec &B, Vec &X0,Prec &PC=*((Prec*)(NULL)))
  {
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver GMRES ");
     #endif
     if (preconditioning && !PC.Type_Of_Preconditioner )   PC.precond_Matrix=&A;
     if (A.type()==2 && B.type()==2)  return GmresRC( real_t(0),A,B,X0,PC );
                                 else return GmresRC(complex_t(0),A,B,X0,PC);
  }
  
/// solve A.X=B using GMRES method with initial guess X0 
///and may be use with PC as a preconditionner.
  template < typename T, class O , class Vec,  class Prec >
  Vec   GmresRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC )
  {
    me_TRACE->push("GmresSolver::GmresRC");
    
    /// Default value of Krylov space dimension if not defined
    if (DimKrylov == 0) DimKrylov=100;
    Vec**      V=new Vec* [DimKrylov+1];
    for (int_t i=0;i<DimKrylov+1;i++)
         V[i]=new Vec (X0);               ///DimKrylov+1 MeVectors, size B.dim (~NVECT)
    Vec    X(X0);
    T*    rs=new T[DimKrylov+1]; 
    T*    Hess=new T[(DimKrylov+1)*DimKrylov];  
         ///Hess is an Hessemberg matrix represented by a T-tab
    T*    c=new T[DimKrylov];
    T*    s=new T[DimKrylov];
    T     t,gam;    

    int_t iter=0;
    *V[0]=B;
    *V[0] -= A*X0;
    real_t rho=(*V[0]).norm2();
    real_t rhoTmp,normOfB=B.norm2();
   // real_t eps=epsilon*B.norm2();
   
    if (rho <= epsilon) {print("GMRES",iter,abs(rho)); return X;}
    while((rho > epsilon) && (iter<MaxIter)  )
    { 
      if (preconditioning) *V[0]=PC.solve(*V[0]);  
      rho=(*V[0]).norm2();     
      rhoTmp=normOfB;//1.;
      rs[0]=rho;
      *V[0] *= (1./rho);
      int_t i=0;
      while( (i < DimKrylov) && (rhoTmp > epsilon))
      {  
        (*V[i+1])= A*(*V[i]);           
        /// solve PC*Y=A*V(i) and V(i)=Y   
        if (preconditioning) *V[i+1]=PC.solve(*V[i+1]);         
        
        /// Arnoldi Method
        /// The inner product coefficients are stored in the upper Hessemberg matrix 
        for(int_t j=0;j<=i;j++)
        {
          equal(Hess[j*(DimKrylov+1)+i] , (*V[i+1]).dotC((*V[j])) );
             ///Hess(j,i)=(*V[i+1]).hermDot(*V[j])
          (*V[i+1]) -= ( Hess[j*(DimKrylov+1)+i]*(*V[j]) ); 
        }
        
        Hess[(i+1)*(DimKrylov+1)+i]=(*V[i+1]).norm2(); 
       (*V[i+1])/=Hess[(i+1)*(DimKrylov+1)+i];   
     
        /// Givens Rotation
        if(i>=1)
        {
          for(int_t k=1;k<=i;k++)
          {
             t=Hess[(k-1)*(DimKrylov+1)+i];
             Hess[(k-1)*(DimKrylov+1)+i]= c[k-1]*t + s[k-1]*Hess[k*(DimKrylov+1)+i];
             Hess[k*(DimKrylov+1)+i] =-s[k-1]*t + c[k-1]*Hess[k*(DimKrylov+1)+i];
          }
        }
        gam=sqrt( pow(Hess[i*(DimKrylov+1)+i],2) +
                  pow(Hess[(i+1)*(DimKrylov+1)+i],2) );
        c[i]   =  Hess[i*(DimKrylov+1)+i]/gam;
        s[i]   =  Hess[(i+1)*(DimKrylov+1)+i]/gam;
        rs[i+1] = -s[i]* rs[i];
        rs[i]  *= c[i];
 
        Hess[i*(DimKrylov+1)+i] = c[i]*Hess[i*(DimKrylov+1)+i] + s[i]*Hess[(i+1)*(DimKrylov+1)+i];
        rhoTmp = abs(rs[i+1]);
        i++;
      } 
      ///end while -> i=DimKrylov or rho<epsilon*||b||
      /// update solution
      i--;
      rs[i]/=Hess[i*(DimKrylov+1)+i];
      for(int_t ii=1;ii<=i;ii++)
      {
        int_t k=i-ii;
        t = rs[k];
        for(int_t j=k+1;j<=i;j++)
            t -= Hess[k*(DimKrylov+1)+j]*rs[j];
        rs[k]=t/Hess[k*(DimKrylov+1)+k];
      }
      for(int_t j=0;j<=i;j++)   X += rs[j]*(*V[j]); 
      
      *V[0]=B - A*X;
      rho=(*V[0]).norm2();
      /// if no CV , restart with X0=X
      if (rho > epsilon)  rs[0]=rho; 
      iter++;
    }   /// end while -> iter=MaxIter or rho<epsilon

    if (abs(rho) < epsilon) print("GMRES",iter,abs(rho));
    else {
       X.print(50);
       ERR_DATA <<"GMRES"<<MaxIter<<epsilon;
       me_error(ERR_SOLVER,"echec",ERR_DATA);
    }    
    /// delete V:
    for (int_t i=0;i<DimKrylov+1;i++)
         delete  V[i];
    delete [] V;    
    
    delete rs; 
    delete Hess;
    delete c;
    delete s;
         
    me_TRACE->pop(); 
    return X;

  }

};



#endif