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

  This class implement the Conjugate Gradient Squared method, preconditionned or not.

  CGS is suitable for non symetric systems
  For more explaination about Conjugate Gradient Squared method,
see http://mathworld.wolfram.com/ConjugateGradientSquaredMethod.html.

Base class : IterativeSolver

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

#ifndef CGSSOLVER_H
#define CGSSOLVER_H

#include "IterativeSolver.h++"

class CgsSolver : public IterativeSolver  {
public:
/// Constructor:
CgsSolver(){IterativeSolver();};
CgsSolver(int_t M, real_t eps=1.e-6) {MaxIter=M; epsilon=eps;}    
CgsSolver(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 CGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("CgsSolver::operator(A,B)");
     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 CGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     preconditioning=false;
     Preconditioner<O> NoPrec;
     me_TRACE->push("CgsSolver::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 CGS in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("CgsSolver::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 CGS ");
     #endif
     me_TRACE->push("CgsSolver::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 CgsRC((real_t) 0,A,B,X0,PC);
                                else  return CgsRC((complex_t) 0,A,B,X0,PC);
}


/// solve A.X=B using CGS method with initial guess X0 and may be use with PC as a preconditionner.
//=*((Prec*)(NULL))
 template < typename T, class O , class Vec,  class Prec >
  Vec   CgsRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
  me_TRACE->push("CgsSolver::CgsRC");

  Vec R0=B;
  R0 -= A*X0;
  Vec X(X0);
  Vec U(X0);
  Vec Utilda(X0);
  Vec Q(X0);
  Vec V(X0);
  Vec P(X0);
  Vec Phat(X0);
  Vec Rtilda=R0;
  Vec UpQ(X0);
  int_t iter=0;
  T beta,rho1,alpha;
  T rho0=R0.norm2();
  
  while ((abs(rho0) > 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)
       {
           U=R0;
           P=U;
       }
       else
       {
           beta=rho1/rho0;
           Q*=beta;
           U=(R0+Q);
           P*=pow(beta,2); P+=U+Q;
       }
       if (abs(rho1) < 1.e-40)  {
             ERR_DATA <<abs(rho1);
             me_error(ERR_SOLVER,"breakdown",ERR_DATA);
       }
       if (preconditioning) Phat=PC.solve(P);
       else Phat=P;
           
       V=(A*Phat);
       equal(alpha , rho1/Rtilda.dotRC(V) );   //alpha=rho1/Rtilda.dot(V);
       Q=(U-alpha*V);   
       if (preconditioning) {UpQ=U+Q;
                             Utilda=PC.solve(UpQ);}      
       else Utilda=U+Q; 
              
       X+=(alpha*Utilda);
       R0+=(-alpha*(A*Utilda));
       rho0=rho1;
       iter++;  
  }
       
  if (abs(rho0) < epsilon) print("CGS",iter,abs(rho0));
  else {
         ERR_DATA << "CGS" << MaxIter << epsilon;
         me_error(ERR_SOLVER,"echec",ERR_DATA);
  }    
  me_TRACE->pop(); 
  return X;
}


};




#endif