fopr_Clover_eo.h

Go to the documentation of this file.

#ifndef FOPR_CLOVER_EO_INCLUDED
#define FOPR_CLOVER_EO_INCLUDED

#include <valarray>
#include <vector>
#include <string>

#include "fopr_Wilson_eo.h"
#include "fopr_CloverTerm_eo.h"

#include "solver_CG.h"

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

//- parameters class
class Parameters_Fopr_Clover_eo : virtual public Parameters
{
 public:
  Parameters_Fopr_Clover_eo();
};
//- end


class Fopr_Clover_eo : public Fopr_eo
{
 private:
  int                m_Nvol, m_Nvol2, m_Ndim, m_Nc, m_Nd, m_NinF;
  double             m_kappa, m_cSW;
  std::valarray<int> m_boundary;
  std::string        m_mode;

  Fopr_Wilson_eo     *m_fopr_w;
  Fopr_CloverTerm_eo *m_fopr_csw;

  Index_eo m_idx;
  Field_G  *m_Ueo;

  ShiftField_eo m_shift_eo;

  Field_F  *m_fee_inv;
  Field_F  *m_foo_inv;
  Vec_SU_N v1, v2;

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

 public:
  Fopr_Clover_eo(std::string repr)
  {
    m_Nvol  = CommonParameters::Nvol();
    m_Nvol2 = m_Nvol / 2;
    m_Nc    = CommonParameters::Nc();
    m_Ndim  = CommonParameters::Ndim();
    m_Nd    = CommonParameters::Nd();
    m_boundary.resize(m_Ndim);

    m_fee_inv = new Field_F(m_Nvol2, m_Nc * m_Nd);
    m_foo_inv = new Field_F(m_Nvol2, m_Nc * m_Nd);

    m_Ueo = new Field_G(m_Nvol, m_Ndim);

    m_fopr_w   = new Fopr_Wilson_eo(repr);
    m_fopr_csw = new Fopr_CloverTerm_eo(repr);
  }

  ~Fopr_Clover_eo()
  {
    delete m_fopr_w;
    delete m_fopr_csw;

    delete m_Ueo;

    delete m_foo_inv;
    delete m_fee_inv;
  }

  void set_parameters(const Parameters& params);
  void set_parameters(const double kappa, const double cSW,
                      const std::valarray<int> bc);
  void set_config(Field *U);

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

    if (m_mode == "D") {
      m_mult     = &Fopr_Clover_eo::D;
      m_mult_dag = &Fopr_Clover_eo::Ddag;
      m_preProp  = &Fopr_Clover_eo::prePropD;
      m_postProp = &Fopr_Clover_eo::postPropD;
    } else if (m_mode == "Ddag") {
      m_mult     = &Fopr_Clover_eo::Ddag;
      m_mult_dag = &Fopr_Clover_eo::D;
      m_preProp  = &Fopr_Clover_eo::prePropDag;
      m_postProp = &Fopr_Clover_eo::postPropDag;
    } else if (m_mode == "DdagD") {
      m_mult     = &Fopr_Clover_eo::DdagD;
      m_mult_dag = &Fopr_Clover_eo::DdagD;
    } else if (m_mode == "H") {
      m_mult     = &Fopr_Clover_eo::H;
      m_mult_dag = &Fopr_Clover_eo::H;
    } else {
      vout.crucial("Fopr_Clover_eo: undefined mode = %s\n", mode.c_str());
      abort();
    }
  }

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

  const Field mult(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    mult(v, f);
    return v;
  }

  const Field mult_dag(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    mult_dag(v, f);
    return v;
  }

  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);
  }

  //- 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&);

  const Field_F mult_csw_inv(const Field_F&, const int ieo);
  const Field_G trSigmaInv(const int mu, const int nu);

  const Field MeoMoe(const Field& f)
  {
    Field_F v((Field_F)f);

    v -= Meo(Meo(f, 1), 0);
    return (Field)v;
  }

  const Field D(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    D(v, f);
    return v;
  }

  void D(Field& v, const Field& f)
  {
    v  = f;
    v -= (Field)Meo(Meo(f, 1), 0);
  }

  const Field Ddag(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    Ddag(v, f);
    return v;
  }

  void Ddag(Field& v, const Field& f)
  {
    v  = f;
    v -= (Field)Mdageo(Mdageo(f, 1), 0);
  }

  const Field DdagD(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    DdagD(v, f);
    return v;
  }

  void DdagD(Field& v, const Field& f)
  {
    v = Ddag(D(f));
  }

  const Field H(const Field& f)
  {
    Field v(f.nin(), f.nvol(), f.nex());

    H(v, f);
    return v;
  }

  void H(Field& v, const Field& f)
  {
    v = mult_gm5(D(f));
  }

  inline const Field_F Meo(const Field_F& f, const int ieo)
  {
    // ieo=0: even <-- odd
    // ieo=1: odd  <-- even
    return mult_csw_inv(m_fopr_w->Meo(f, ieo), ieo);
  }

  inline const Field_F Meo_gm5(const Field_F& f, const int ieo)
  {
    return mult_gm5(Meo(f, ieo));
  }

  inline const Field_F Mdageo(const Field_F& f, const int ieo)
  {
    return mult_csw_inv(m_fopr_w->Mdageo(f, ieo), ieo);
  }

  inline const Field mult_gm5(const Field& f)
  {
    return m_fopr_w->mult_gm5(f);
  }

  inline void mult_isigma(Field_F& w, const Field_F& f,
                          const int mu, const int nu)
  {
    return m_fopr_csw->mult_isigma(w, f, mu, nu);
  }

  inline std::vector<double> csmatrix(const int& site)
  {
    return m_fopr_csw->csmatrix(site);
  }

  int field_nvol() { return CommonParameters::Nvol() / 2; }
  int field_nin() { return 2 * CommonParameters::Nc() * CommonParameters::Nd(); }
  int field_nex() { return 1; }

 private:
  void solve_csw_inv();
};
#endif