/**
================================================================================
=                       Class BicgSolver header file                         =
================================================================================
Author : C.Chambeyron (24 november 2004)
Version: 1
================================================================================
  This class implement the BiConjugate Gradient method, preconditionned or not.  
  BiCG is not suitable for non symetric systems

 For more explaination about Conjugate Gradient Squared method,
see http://mathworld.wolfram.com/BiconjugateGradientMethod.html.

Base class : IterativeSolver

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

#ifndef BICGSOLVER_H
#define BICGSOLVER_H

#include "IterativeSolver.h++"
#include "UtilDef.h++"


class BicgSolver : public IterativeSolver  {
public:
/// Constructors:
BicgSolver(){IterativeSolver();};
BicgSolver(int_t M, real_t eps=1.e-6) {MaxIter=M; epsilon=eps;}    
BicgSolver(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 BICG in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("BiCgSolver::operator(A,B,X0)");
     me_TRACE->pop();
     preconditioning=false;
     Preconditioner<O> NoPrec;
     Vec X0(B);
     X0.setEntries(0.);
     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 BiCG in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     preconditioning=false;
     Preconditioner<O> NoPrec;
     me_TRACE->push("BiCGSolver::operator(A,B,X0)");
     me_TRACE->pop();
     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 BICG in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("BicgSolver::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)
{
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver BiCG ");
     #endif
     me_TRACE->push("BiCGSolver::operator(A,B,X0,(PC) )");
     me_TRACE->pop();
     if (preconditioning && !PC.Type_Of_Preconditioner )   PC.precond_Matrix=&A; 
     if (A.type()==2 && B.type()==2)   return BicgRC(real_t(0.),A,B,X0,PC);
                                  else return BicgRC(complex_t(0.),A,B,X0,PC);
}

/// solve A.X=B using BiCG method with initial guess X0 and may be use with PC as a preconditionner.
 template < typename T, class O , class Vec,  class Prec >
  Vec   BicgRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
   me_TRACE->push("BiCGSolver::BicgRC");   
   Vec X(X0);
   Vec Q(X0);
   Vec Qtilda(X0);
   Vec P(X0);
   Vec Ptilda(X0);
   Vec R0 = B - A*X0;
   Vec R0tilda(R0);
   Vec Z(X0);
   Vec Ztilda(X0); 
   T beta;
   T alpha;
   int_t iter=0;
   T rhoold;
   T rho; 
   T gamma;       equal(gamma , R0.dotRC(R0) );                   //gamma=R0.dot(R0);   
   T gammatilda;  equal(gammatilda , R0tilda.dotRC(R0tilda) );    // gammatilda=R0tilda.dot(R0tilda);
   
   while( abs(gammatilda)>abs(epsilon) && abs(gamma)>abs(epsilon) && (iter<MaxIter))
   {      
        if (preconditioning) {Z=PC.solve(R0);
                              Ztilda=PC.transposedSolve(R0tilda);}  
        else {Z=R0; Ztilda=R0tilda;}   
           
        equal(rho , Z.dotRC(R0tilda) );   //rho=Z.dot(R0tilda)
        if (abs(rho) < 1.e-40) {
             ERR_DATA <<abs(rho);
             me_error(ERR_SOLVER,"breakdown",ERR_DATA);
        }
        if (iter==0) {
           P=Z;
           Ptilda=Ztilda;
        }
        else{
           beta=rho/rhoold;
           P*=beta;
           P+=Z;
           Ptilda*=beta;
           Ptilda+=Ztilda;
        }
        Q=A*P;
        Qtilda=Ptilda*A; //A^t*Ptilda;

        equal(alpha , rho/Ptilda.dotRC(Q) );   //alpha=rho/Ptilda.dot(Q);
        X+=alpha*P;
        R0-=alpha*Q;
        R0tilda-=alpha*Qtilda;
        equal(gamma , R0.dotRC(R0) );    //gamma=R0.dot(R0);

        rhoold=rho;
        iter++;
    }

    if (abs(gamma) < epsilon) print("BiCG",iter,abs(gamma));
    else  {
        ERR_DATA <<"BiCG"<<MaxIter<<epsilon;
        me_error(ERR_SOLVER,"echec",ERR_DATA);
    }   
    
  me_TRACE->pop();  
  return X;
}


};




#endif