/**
============================================
                          SorSolver.h  -  description
                             -------------------
    begin                : lun jan 10 2005
===========================================
**/

#ifndef SORSOLVER_H
#define SORSOLVER_H

#include "IterativeSolver.h++"
#include "me_Trace.h"
#include "UtilDef.h++"
#include "Preconditioner.h++"
#include "me_Config.h++"


class SorSolver : public IterativeSolver 
 {
public:
  real_t Omega;       ///relaxation factor of SOR method
//constructors:
SorSolver(){IterativeSolver();Omega=1.;};
     
SorSolver(int_t M, real_t eps=1.e-6,real_t w=1.)
        {MaxIter=M; epsilon=eps;Omega=w;}
     
SorSolver(real_t eps,real_t w, int_t M=100)
        {MaxIter=M; epsilon=eps;Omega=w;}
           
SorSolver( real_t eps,int_t M,real_t w=1.)
        {MaxIter=M; epsilon=eps;Omega=w;}  
      
template < class O, class Vec  >
Vec   operator()(O &A,Vec &B)
{ 
     #ifdef ME_DEBUG_ON
        me_TRACE->log(" Solver Sor in  ",MaxIter," iterations max and tolerance is ",epsilon);
     #endif
     me_TRACE->push("SorSolver::operator(A,B,X0)");
     me_TRACE->pop();
     Vec X0(B);
     X0.setEntries(0.); 
     preconditioning=false;
     Preconditioner<O> NoPrec;
     return (*this)(A,B,X0,NoPrec /* *((O*)(NULL)) */);
}         
   
template < class O, class Vec  >
Vec   operator()(O &A,Vec &B, Vec &X0)
{      
     #ifdef ME_DEBUG_ON
       me_TRACE->log(" Solver Sor 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 Sor 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 Sor ");
     #endif
     if (preconditioning && !PC.Type_Of_Preconditioner )   PC.precond_Matrix=&A;
     if (A.type()==2 && B.type()==2)  return SorRC(real_t(0.),A,B,X0,PC);
                                 else return SorRC(complex_t(0.),A,B,X0,PC);
}

/// solve A.X=B using SOR method with initial guess X0 
///and may be use with PC as a preconditionner.
template < typename T, class O , class Vec,  class Prec >
Vec   SorRC(const T ty,O &A,Vec &B, Vec &X0,Prec &PC)
{
  me_TRACE->push("SorSolver::SorRC");
  Vec R0=B;
  R0-=A*X0;
  Vec X(X0);
  Vec  Y(X0);
  Vec TT(X0);
  int_t iter=0;
  real_t rho=R0.norm2();
  
  if ((Omega<=0) || (Omega>=2)) {cout<<"Omega doit etre dans ]0;2["<<endl; abort();}
  while ((rho > epsilon) && (iter<MaxIter)) {    
/// find X(iter+1)= [D/w+L]^-1.[(1/w-1)D-U] *X(iter) + B  
        Y=A.upperPartMultiplication(X,1.-1./Omega);
        Y-=B;
        Y*=(-1.);
        X=A.lowerPartSolver(Y,1./Omega);
/// update:        
        R0=B - A*X;
        rho=R0.norm2();
        iter++;
  }
  if (abs(rho) < epsilon) print("Sor",iter,abs(rho));
  else {
       ERR_DATA <<"SOR"<<MaxIter<<epsilon;
       me_error(ERR_SOLVER,"echec",ERR_DATA);
  }     
  me_TRACE->pop(); 
  return X;
}

};    ///end of class SorSolver
#endif