/**
================================================================================
=                       Class BicgsSolver header file                         =
================================================================================
Author : C.Chambeyron (24 november 2004)
Version: 1
================================================================================

  This class implement th BiConjugate Gradient Stabilized method, preconditionned or not.

  Bicgs is suitable for non symetric systems
 For more explaination about BiConjugate Gradient Stabilized method,
see http://mathworld.wolfram.com/BiconjugateGradientStabilizedMethod.html.

Base class : IterativeSolver

================================================================================
**/

#ifndef BICGSSOLVER_H
#define BICGSSOLVER_H


#include "IterativeSolver.h++"

class BicgsSolver : public IterativeSolver  {
public:
/// Constructor:
BicgsSolver(){IterativeSolver();};
BicgsSolver(int_t M, real_t eps=1.e-6) {MaxIter=M; epsilon=eps;}    
BicgsSolver(real_t eps, int_t M=100)   {MaxIter=M; epsilon=eps;}


template < class O, class Vec  >
  Vec   operator()(O &A,Vec &B)
  { 
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver BICGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("BICgsSolver::operator(A,B,X0)");
     me_TRACE->pop();
     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 BiCGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("BicgsSolver::operator(A,B,X0)");
     me_TRACE->pop(); 
     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 BICGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("BICgsSolver::operator(A,B,X0)");
     me_TRACE->pop();
     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 BiCGS ");
     #endif
     me_TRACE->push("BicgsSolver::operator(A,B,X0,(PC) )");
     me_TRACE->pop();
     if (preconditioning && !PC.precond_Matrix )   PC.precond_Matrix=&A;
     if (A.type()==2 && B.type()==2)  return BicgsRC(real_t(0.),A,B,X0,PC);
                                 else return BicgsRC(complex_t(0.),A,B,X0,PC);
}
/// solve A.X=B using BiCGS method with initial guess X0 and may be use with PC as a preconditionner.
 template < typename T, class O , class Vec,  class Prec >
  Vec   BicgsRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
  me_TRACE->push("BicgsSolver::BicgsRC");
  Vec R0=B;
  R0 -= (A*X0);
  Vec Rtilda=R0;
  Vec X(X0);
  Vec V(X0);
  Vec P(X0);
  Vec Ptilda(X0);
  Vec S(X0);
  Vec Stilda(X0);
  Vec TT(X0);

  int_t iter=0;
  real_t rho=R0.norm2();
  T beta,rho1,alpha,omega;
  T rho0=0.;
  T RtildaV;  
  
  while ((abs(rho) > epsilon) && (iter<MaxIter)) {
       equal (rho1 , Rtilda.dotRC(R0));       //rho1=Rtilda.dot(R0);
       if (abs(rho1) <= 1.e-40)  {
             ERR_DATA <<abs(rho1);
             me_error(ERR_SOLVER,"breakdown",ERR_DATA);
       }

       if (iter == 0) P=R0;
       else
       {
           beta=(rho1*alpha)/(rho0*omega);
           P*=beta;
           P+=R0-(beta*omega)*V;
       }
       /// solve PC*Ptilda=P    
       if (preconditioning) Ptilda=PC.solve(P);  
       else Ptilda=P;  
       
       V=A*Ptilda;
       equal(RtildaV , Rtilda.dotRC(V));   //RtildaV=Rtilda.dot(V);
       if (abs(RtildaV) <= 1.e-40) {
             ERR_DATA <<abs(RtildaV);
             me_error(ERR_SOLVER,"breakdown",ERR_DATA);
       }
       alpha=rho1/RtildaV;
       S=R0-alpha*V;
       if (S.norm2()<epsilon)
       {
          X+=alpha*Ptilda;
          print("BiCGS",iter,S.norm2());
          return X;  
       }
       if (preconditioning) Stilda=PC.solve(S); 
       else Stilda=S;
       /// solve PC*Stilda=S
       TT=A*Stilda;
       equal(omega , TT.dotRC(S)/TT.dotRC(TT));     //omega=TT.dot(S)/TT.dot(TT);
       if (abs(omega) <= 1.e-40) {
             ERR_DATA <<abs(omega);
             me_error(ERR_SOLVER,"breakdown",ERR_DATA);
       }
       X+=alpha*Ptilda+omega*Stilda;
       R0=S-omega*TT;
       rho=R0.norm2();
       rho0=rho1;
       iter++;
  } 
  if (abs(rho) < epsilon) print("BiCGS",iter,abs(rho));
  else {
      ERR_DATA<<"BICGS"<<MaxIter<<epsilon;
      me_error(ERR_SOLVER,"echec",ERR_DATA);
  }
    
  me_TRACE->pop(); 
  return X;
}


};




#endif