cg.hh

Go to the documentation of this file.

 
#ifndef cg_hh
#define cg_hh
 
#include <string>
 
#include "basics/typedefs.hh"
#include "basics/parallel.hh"
#include "compositions.hh"
 
#ifdef HAS_MPI
#include <mpi.h>
#endif
 
namespace concepts {
 
  // ******************************************************************* CG **
 
  template<class F>
  class CG : public VecOperator<F> {
  public:
    CG(Operator<F>& A, Real maxeps, int maxit = 0, uint relres = false,
       bool throwing = true ) 
      : VecOperator<F>(A.dimY(), A.dimX()), W_(0), A_(A)
      , maxeps_(maxeps*maxeps), maxit_(maxit), eps_(1.0), it_(0)
      , relres_(relres), throwing_(throwing) {
#ifdef HAS_MPI
      // Parallel compiled environment 
      if (isParallelRunning ())
      {
        // Parallel running environment
        MPI_Comm_size(MPI_COMM_WORLD, &N_);
        MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
        parallel_ = (N_>1);
        if (parallel_) J_ = create_index_vector();
      }
      else
      {
        // Sequential running environment
        N_ = 1; rank_ = 0;
        parallel_ = false;
      }
#else
      // Sequential compiled environment
      N_ = 1; rank_ = 0;
      parallel_=false;
#endif
    }
 
    CG(Operator<F>& A, Operator<F>& Minv, Real maxeps, int maxit = 0, 
       bool relres = 0, bool throwing = true)
      : VecOperator<F>(A.dimY(), A.dimX()), W_(&Minv), A_(A)
      , maxeps_(maxeps*maxeps), maxit_(maxit), eps_(1.0), it_(0)
      , relres_(relres), throwing_(throwing) {
#ifdef HAS_MPI
      // Parallel compiled environment 
      if (isParallelRunning ())
      {
        // Parallel running environment
        MPI_Comm_size(MPI_COMM_WORLD, &N_);
        MPI_Comm_rank(MPI_COMM_WORLD, &rank_);
        parallel_ = (N_>1);
        if (parallel_) J_ = create_index_vector();
      }
      else
      {
        // Sequential running environment
        N_ = 1; rank_ = 0;
        parallel_ = false;
      }
#else
      // Sequential compiled environment
      N_ = 1; rank_ = 0;
      parallel_=false;
#endif
    }
 
    uint iterations() const { return it_; }
 
    Real epsilon() const { return eps_; }
  protected:
    std::ostream& info(std::ostream& os) const;
 
  private:
    Operator<F>*  W_;
    Operator<F>&  A_;
 
    Real maxeps_;
    uint maxit_;
    Real eps_;
    uint it_;
    bool relres_;
    bool throwing_;
    
    int N_;
    
    int rank_;
 
    bool parallel_;
 
    Vector<Real> J_;
    
    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX);
    virtual void apply_();
    
    Vector<Real> create_index_vector(); 
    
#ifdef HAS_MPI
    Real parallel_l2_2(const std::map<int,Vector<F>> &K_ij, const Vector<F>& r_i,std::map<int,Vector<F>> &map_tilde_r_ji);
#endif
  };
 
} // namespace concepts
 
#endif // cg_hh