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

  This class implement the Conjugate Gradient method, preconditionned or not.
The preconditionner is a default argument.
  CG is not suitable for non symetric define positive  systems

 For more explaination about BiConjugate Gradient Stabilized method,
see http://mathworld.wolfram.com/ConjugateGradientMethod.html.

Base class : IterativeSolver

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

#ifndef CGSOLVER_H
#define CGSOLVER_H

#include "IterativeSolver.h++"

class CgSolver : public IterativeSolver {
public:
/// Constructor:
CgSolver(){IterativeSolver();};
CgSolver(int_t M, real_t eps=1.e-6) {MaxIter=M; epsilon=eps;}    
CgSolver(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 CG in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("CgSolver::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 CG 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 CG 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)
{
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver CG ");
     #endif  
     if (preconditioning && !PC.Type_Of_Preconditioner )   PC.precond_Matrix=&A; 
     if (A.type()==2 && B.type()==2)  return CgRC(real_t(0.),A,B,X0,PC);
                                else  return CgRC(complex_t(0.),A,B,X0,PC);
}

/// solve A.X=B using CG method with initial guess X0 and may be use with PC as a preconditionner.
template < typename T, class O , class Vec,  class Prec >
Vec   CgRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
  me_TRACE->push("CgSolver::CgRC"); 
   Vec X(X0);
   Vec R(B);
   Vec Z(X0);
   Vec P(X0);
   Vec Q(X0);
   T alpha,beta,rho1,rho0;
   int_t iter=0; 

   R-=A*X0;   
   while( (R.norm2()>abs(epsilon)) && (iter<MaxIter))
   { 
     if (preconditioning) Z=PC.solve(R);     
     else Z=R;

     equal(rho1 , R.dotRC(Z));                        //rho1=R.dot(Z);
     if (iter == 0) {P=Z;}
     else { beta=rho1/rho0;
            P*=beta; P+=Z;}                           //P=Z+beta*P
     Q=A*P;
     equal(alpha , rho1/P.dotRC(Q) );                 //alpha=rho1/P.dot(Q)
     X+=alpha*P;
     R-=alpha*Q;
     rho0=rho1;    
     iter++;
    }
    if (iter >= MaxIter) {
       ERR_DATA << "CG" << MaxIter << epsilon;
       me_error(ERR_SOLVER,"echec",ERR_DATA);
    }
    else print("CG",iter,R.norm2());
   
   me_TRACE->pop(); 
   return X;
}

};
#endif