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

  This class implement the QMR method, preconditionned or not.
  QMR is not suitable for non symetric systems

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

Base class : IterativeSolver

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

#ifndef QMRSOLVER_H
#define QMRSOLVER_H

#include "IterativeSolver.h++"

class QmrSolver : public IterativeSolver  {
public:

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

/// solve A.X=B using QMR method with initial guess X0 and may be use with PC as a preconditionner.
 template < typename T, class O , class Vec,  class Prec >
  Vec   QmrRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
  me_TRACE->push("QmrSolver::QmrRC");
  T no=ty;  no*=0.;      // to erase advertissement
  Vec R0=B - A*X0;
  Vec X(X0);
  Vec V(R0);
  Vec W(R0);
  Vec P(X0);
  Vec Q(X0);
  Vec Ptilda(X0);
  Vec D(X0);
  Vec S(X0);
  real_t ksi;
  T gamma0=1.;
  T etha=-1.;
  real_t rho0,rho1;
  T      gamma1,delta,eps,beta,tetha1,tetha0;
  int_t iter=0;
  real_t rhoo=R0.norm2();

  if(!preconditioning) {
// solve without preconditionning:
              rho0=V.norm2(); ksi=W.norm2();
              while ((rhoo > epsilon) && (iter<MaxIter)) {
                    if (rho0 < 1.e-40)  {
                           ERR_DATA <<abs(rho0);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }
                    V/=rho0;
                    W/=ksi;
                    equal(delta , V.dotRC(W) );    //delta=V.dot(W)
                    if (abs(delta) < 1.e-40)  {
                           ERR_DATA <<abs(delta);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }

                    if (iter == 0) {
                        P=V;
                        Q=W;
                    }
                    else{
                        P*=(-ksi*delta/eps);  P+=V;
                        Q*=(-rho0*delta/eps);  Q+=W;
                    }

                    Ptilda=A*P;
                    equal(eps , Q.dotRC(Ptilda) );   //eps=Q.dot(Ptilda)
                    if (abs(eps) < 1.e-40)  {
                           ERR_DATA <<abs(eps);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }

                    beta=eps/delta;
                    if (abs(beta) < 1.e-40)  {
                           ERR_DATA <<abs(beta);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }

                    V*=(-beta); V+=Ptilda;
                    rho1=V.norm2();
                    W*=(-beta);  W+=Q*A;   //W+=A^t*Q 
                    ksi=W.norm2();

                    tetha1=(rho1)/(gamma0*abs(beta));
                    gamma1=T(1.)/ sqrt(T(1.)+pow(tetha1,2));
                    if (abs(gamma1) < 1.e-40)  {
                           ERR_DATA <<abs(gamma1);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }
                    etha*=(-rho0*pow(gamma1/gamma0,2)/beta);

                    if (iter == 0) {
                        D=etha*P;
                        S=etha*Ptilda;
                    }
                    else{
                        D*=pow(tetha0*gamma1,2);  D+=etha*P;
                        S*=pow(tetha0*gamma1,2);  S+=etha*Ptilda;
                    }

            /// update
                    X+=D;
                    R0-=S;
                    rho0=rho1;
                    gamma0=gamma1;
                    tetha0=tetha1;
                    iter++;
                    rhoo=R0.norm2();    
              }
 
      if (abs(rhoo) < epsilon) print("QMR",iter,abs(rhoo));
      else {              
             ERR_DATA <<"QMR"<<MaxIter<<epsilon;
             me_error(ERR_SOLVER,"echec",ERR_DATA);
      }
   }
   else { 
// solve with preconditionner :
              Vec Y(B);
              Vec Yhat(B);
              Vec Z(B);
              Vec Zhat(B);
              
              Y=PC.solve(V); 
              Z=PC.transposedSolve(W);
              
              rho0=Y.norm2(); ksi=Z.norm2();

              while ((rhoo > epsilon) && (iter<MaxIter)) {
                    if (rho0 < 1.e-40)  {
                           ERR_DATA <<abs(rho0);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }
                    V/=rho0; Y/=rho0;
                    W/=ksi; Z/=ksi;
                    equal(delta , Y.dotRC(Z) );     //delta=Y.dot(Z)
                    if (abs(delta) < 1.e-40)  {
                           ERR_DATA <<abs(delta);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }                    
                    
                   Yhat=PC.solve(Y);
                   Zhat=PC.transposedSolve(Z); /// W ou Z????                                                                
                    if (iter == 0) {
                        P=Yhat;
                        Q=Zhat;
                    }
                    else{
                        P*=(-ksi*delta/eps);  P+=Yhat;
                        Q*=(-rho0*delta/eps);  Q+=Zhat;
                    }

                    Ptilda=A*P;
                    equal(eps , Q.dotRC(Ptilda) );   //eps=Q.dot(Ptilda)
                    if (abs(eps) < 1.e-40)  {
                           ERR_DATA <<abs(eps);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }

                    beta=eps/delta;
                    if (abs(beta) < 1.e-40)  {
                           ERR_DATA <<abs(beta);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }

                    V*=(-beta);  V+=Ptilda;                      
 
                    W*=(-beta);  W+=Q*A;  //W+=A^t*Q
                    Y=PC.solve(V);
                    Z=PC.transposedSolve(W); 
                    
                    rho1=Y.norm2();                    
                    ksi=Z.norm2();
                    tetha1=rho1/(gamma0*abs(beta));
                    gamma1=T(1.)/(sqrt(T(1.)+pow(tetha1,2)));
                    if (abs(gamma1) < 1.e-40)  {
                           ERR_DATA <<abs(gamma1);
                           me_error(ERR_SOLVER,"breakdown",ERR_DATA);
                    }
                    etha*= (-rho0*pow(gamma1/gamma0,2)/beta);

                    if (iter == 0) {
                        D=etha*P;
                        S=etha*Ptilda;
                    }
                    else{
                        D*=pow(tetha0*gamma1,2);  D+=etha*P;
                        S*=pow(tetha0*gamma1,2);  S+=etha*Ptilda;
                    }
            /// update
                    X+=D;
                    R0-=S;
                    rho0=rho1;
                    gamma0=gamma1;
                    tetha0=tetha1;
                    iter++;
                    rhoo=R0.norm2();
              }
     if (abs(rhoo) < epsilon) print("Qmr preconditionned",iter,abs(rhoo));
     else {
             ERR_DATA <<"QMR preconditionned"<<MaxIter<<epsilon;
             me_error(ERR_SOLVER,"echec",ERR_DATA);             
      }  
   }// end if      
   me_TRACE->pop();      
  return X;

}


};

#endif