fopr_Wilson_eo_impl.h

Go to the documentation of this file.

#ifndef FOPR_WILSON_EO_IMP_BGQ_INCLUDED
#define FOPR_WILSON_EO_IMP_BGQ_INCLUDED

#include "Fopr/fopr_eo.h"
#include "Field/index_eo.h"
#include "Communicator/BGNET/channel.h"

#include "IO/bridgeIO.h"
using Bridge::vout;


namespace Imp_BGQ {
  class Fopr_Wilson_eo : public Fopr_eo
  {
   public:
    static const std::string class_name;

   private:
    int m_Nc, m_Nd, m_Nvc, m_Ndf;
    int m_Nvol, m_Nvol2, m_Ndim;
    int m_Nx, m_Ny, m_Nz, m_Nt, m_Nx2;

    double              m_kappa;     
    std::vector<int>    m_boundary;  
    std::vector<double> m_boundary2; 

    std::string m_mode;              
    std::string m_repr;              

    Index_eo m_index;
    Field_G  *m_Ueo;
    Field_G  *m_U;  

    std::vector<int> m_Leo;
    Field            m_v1, m_v2, m_w1, m_w2; 

    double *vcp1_xp, *vcp2_xp, *vcp1_xm, *vcp2_xm;
    double *vcp1_yp, *vcp2_yp, *vcp1_ym, *vcp2_ym;
    double *vcp1_zp, *vcp2_zp, *vcp1_zm, *vcp2_zm;
    double *vcp1_tp, *vcp2_tp, *vcp1_tm, *vcp2_tm;

    // use async data transfer
    std::vector<int>       m_npe;
    std::vector<Channel *> m_fw_send;
    std::vector<Channel *> m_fw_recv;
    std::vector<Channel *> m_bw_send;
    std::vector<Channel *> m_bw_recv;

    void (Fopr_Wilson_eo::*m_gm5)(Field&, const Field&);
    void (Fopr_Wilson_eo::*m_gm5_self)(Field&);
    void (Fopr_Wilson_eo::*m_mult_tp)(Field&, const Field&, const int ieo);
    void (Fopr_Wilson_eo::*m_mult_tm)(Field&, const Field&, const int ieo);

    void (Fopr_Wilson_eo::*m_mult)(Field&, const Field&);
    void (Fopr_Wilson_eo::*m_mult_dag)(Field&, const Field&);
    void (Fopr_Wilson_eo::*m_D)(Field&, const Field&);
    void (Fopr_Wilson_eo::*m_preProp)(Field&, Field&, const Field&);
    void (Fopr_Wilson_eo::*m_postProp)(Field&, const Field&, const Field&);

   public:

    Fopr_Wilson_eo()
    { init("Dirac"); }

    Fopr_Wilson_eo(std::string repr)
    { init(repr); }

    ~Fopr_Wilson_eo()
    { tidyup(); }

    void set_parameters(const Parameters& params);
    void set_parameters(const double kappa, const std::vector<int> bc);

    void set_config(Field *U);

    void set_config(unique_ptr<Field_G>& U)
    {
      set_config(U.get());
    }

    // method for even odd fermion operator
    void preProp(Field& Be, Field& bo, const Field& b)
    { (this->*m_preProp)(Be, bo, b); }

    void postProp(Field& x, const Field& xe, const Field& bo)
    { (this->*m_postProp)(x, xe, bo); }

    void prePropD(Field&, Field&, const Field&);
    void postPropD(Field&, const Field&, const Field&);
    void prePropDag(Field&, Field&, const Field&);
    void postPropDag(Field&, const Field&, const Field&);

    void mult(Field& v, const Field& f)
    { (this->*m_mult)(v, f); }

    void mult_dag(Field& v, const Field& f)
    { (this->*m_mult_dag)(v, f); }

    void set_mode(std::string mode)
    {
      m_mode = mode;

      if (m_mode == "D") {
        m_mult     = &Fopr_Wilson_eo::D;
        m_mult_dag = &Fopr_Wilson_eo::Ddag;
        m_preProp  = &Fopr_Wilson_eo::prePropD;
        m_postProp = &Fopr_Wilson_eo::postPropD;
      } else if (m_mode == "Ddag") {
        m_mult     = &Fopr_Wilson_eo::Ddag;
        m_mult_dag = &Fopr_Wilson_eo::D;
        m_preProp  = &Fopr_Wilson_eo::prePropDag;
        m_postProp = &Fopr_Wilson_eo::postPropDag;
      } else if (m_mode == "DdagD") {
        m_mult     = &Fopr_Wilson_eo::DdagD;
        m_mult_dag = &Fopr_Wilson_eo::DdagD;
      } else if (m_mode == "DDdag") {
        m_mult     = &Fopr_Wilson_eo::DDdag;
        m_mult_dag = &Fopr_Wilson_eo::DDdag;
      } else if (m_mode == "H") {
        m_mult     = &Fopr_Wilson_eo::H;
        m_mult_dag = &Fopr_Wilson_eo::H;
      } else {
        vout.crucial("Error at %s: input mode is undefined.\n", class_name.c_str());
        exit(EXIT_FAILURE);
      }
    }

    std::string get_mode() const
    { return m_mode; }

    void D(Field& v, const Field& f);
    void Ddag(Field& v, const Field& f);
    void DdagD(Field& v, const Field& f);
    void DDdag(Field& v, const Field& f);
    void H(Field& v, const Field& f);

    // ieo=0: even <-- odd
    // ieo=1: odd  <-- even

    void Meo(Field&, const Field&, const int ieo);
    void Mdageo(Field&, const Field&, const int ieo);
    void MeoMoe(Field&, const Field&);
    void Meo_gm5(Field&, const Field&, const int ieo);

    void mult_gm5(Field&, const Field&);
    void mult_gm5(Field&);
    void gm5_dirac(Field&, const Field&);
    void gm5_chiral(Field&, const Field&);
    void gm5_self_dirac(Field&);
    void gm5_self_chiral(Field&);

    void gm5p(const int mu, Field&, const Field& v);

    int field_nvol() { return m_Nvol2; }
    int field_nin()  { return m_Nvc * m_Nd; }
    int field_nex()  { return 1; }

    double flop_count();

   private:
    void init(const std::string);
    void tidyup();

    void mult_xp(Field&, const Field&, const int ieo);
    void mult_xm(Field&, const Field&, const int ieo);
    void mult_yp(Field&, const Field&, const int ieo);
    void mult_ym(Field&, const Field&, const int ieo);
    void mult_zp(Field&, const Field&, const int ieo);
    void mult_zm(Field&, const Field&, const int ieo);
    void mult_tp_dirac(Field&, const Field&, const int ieo);
    void mult_tm_dirac(Field&, const Field&, const int ieo);
    void mult_tp_chiral(Field&, const Field&, const int ieo);
    void mult_tm_chiral(Field&, const Field&, const int ieo);

    void clear_impl(Field&);
    void scal_impl(Field&, double);

    // member data for threading
    int m_Mz, m_Mt;
    int m_Nthread, m_Ntask;
    int m_Ntask_z, m_Ntask_t;

    struct mult_arg
    {
      int isite;
      int isite_cpx, isite_cpy, isite_cpz, isite_cpt;
      int kz0, kz1, kt0, kt1;
    };
    std::vector<mult_arg> m_arg;

    // member functions for threading
    void setup_thread();
    void mult_xp1_thread(int, double *, const double *, int);
    void mult_xp2_thread(int, double *, const double *, int);
    void mult_xpb_thread(int, double *, const double *, int);
    void mult_xm1_thread(int, double *, const double *, int);
    void mult_xm2_thread(int, double *, const double *, int);
    void mult_xmb_thread(int, double *, const double *, int);
    void mult_yp1_thread(int, double *, const double *, int);
    void mult_yp2_thread(int, double *, const double *, int);
    void mult_ypb_thread(int, double *, const double *, int);
    void mult_ym1_thread(int, double *, const double *, int);
    void mult_ym2_thread(int, double *, const double *, int);
    void mult_ymb_thread(int, double *, const double *, int);
    void mult_zp1_thread(int, double *, const double *, int);
    void mult_zp2_thread(int, double *, const double *, int);
    void mult_zpb_thread(int, double *, const double *, int);
    void mult_zm1_thread(int, double *, const double *, int);
    void mult_zm2_thread(int, double *, const double *, int);
    void mult_zmb_thread(int, double *, const double *, int);
    void mult_tp1_dirac_thread(int, double *, const double *, int);
    void mult_tp2_dirac_thread(int, double *, const double *, int);
    void mult_tpb_dirac_thread(int, double *, const double *, int);
    void mult_tm1_dirac_thread(int, double *, const double *, int);
    void mult_tm2_dirac_thread(int, double *, const double *, int);
    void mult_tmb_dirac_thread(int, double *, const double *, int);
    void mult_tp1_chiral_thread(int, double *, const double *, int);
    void mult_tp2_chiral_thread(int, double *, const double *, int);
    void mult_tpb_chiral_thread(int, double *, const double *, int);
    void mult_tm1_chiral_thread(int, double *, const double *, int);
    void mult_tm2_chiral_thread(int, double *, const double *, int);
    void mult_tmb_chiral_thread(int, double *, const double *, int);
    void scal_thread(int, double *, double);
    void clear_thread(int, double *);
    void gm5_dirac_thread(int, double *, const double *);
    void gm5_chiral_thread(int, double *, const double *);
    void gm5_dirac_thread(int, double *);
    void gm5_chiral_thread(int, double *);
  };
}
#endif