fopr_Wilson_impl.cpp

Go to the documentation of this file.

#include "fopr_Wilson_impl.h"

using std::valarray;
using std::string;

//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::init(string repr)
{
  m_vl = CommonParameters::Vlevel();

  m_Nvol = CommonParameters::Nvol();
  m_Ndim = CommonParameters::Ndim();
  m_boundary.resize(m_Ndim);

  m_U = 0;

  m_repr = repr;

  m_GM.resize(m_Ndim + 1);

  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);

  m_GM[0] = gmset->get_GM(GammaMatrixSet::GAMMA1);
  m_GM[1] = gmset->get_GM(GammaMatrixSet::GAMMA2);
  m_GM[2] = gmset->get_GM(GammaMatrixSet::GAMMA3);
  m_GM[3] = gmset->get_GM(GammaMatrixSet::GAMMA4);
  m_GM[4] = gmset->get_GM(GammaMatrixSet::GAMMA5);

  m_mult     = &Fopr_Wilson::Fopr_Wilson_impl::mult_undef;
  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::mult_undef;

  delete gmset;
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::set_mode(string mode)
{
  m_mode = mode;

  if (m_mode == "D") {
    m_mult     = &Fopr_Wilson::Fopr_Wilson_impl::D;
    m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::Ddag;
  } else if (m_mode == "Ddag") {
    m_mult     = &Fopr_Wilson::Fopr_Wilson_impl::Ddag;
    m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::D;
  } else if (m_mode == "DdagD") {
    m_mult     = &Fopr_Wilson::Fopr_Wilson_impl::DdagD;
    m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::DdagD;
  } else if (m_mode == "H") {
    m_mult     = &Fopr_Wilson::Fopr_Wilson_impl::H;
    m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::H;
  } else {
    vout.crucial(m_vl, "Fopr_Wilson_impl: input mode is undefined.\n");
    abort();
  }
}


//====================================================================
string Fopr_Wilson::Fopr_Wilson_impl::get_mode() const
{
  return m_mode;
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::set_parameters(const double kappa, const valarray<int> bc)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of Wilson fermion operator:\n");
  vout.general(m_vl, "  kappa = %8.4f\n", kappa);
  for (int mu = 0; mu < m_Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

  //- range check
  // NB. kappa = 0 is allowed.
  assert(bc.size() == m_Ndim);

  //- store values
  m_kappa = kappa;

  // m_boundary.resize(m_Ndim);  // already resized in init.
  for (int mu = 0; mu < m_Ndim; ++mu) {
    m_boundary[mu] = bc[mu];
  }
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::D(Field& v, const Field& f)
{
  Field_F w(f.nvol(), f.nex());

  v = 0.0;
  w = 0.0;
  for (int mu = 0; mu < m_Ndim; ++mu) {
    mult_up(mu, w, f);
    mult_dn(mu, w, f);
  }

  v = (Field)w;

  v *= -m_kappa;
  v += f;
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::mult_gm5(Field& v, const Field& f)
{
  assert(v.nvol() == f.nvol());
  assert(v.nex() == f.nex());
  assert(v.nin() == f.nin());

  Field_F vt(f.nvol(), f.nex());

  vt.mult_GM(m_GM[4], (Field_F)f);
  v = (Field)vt;
}


//====================================================================
const Field_F Fopr_Wilson::Fopr_Wilson_impl::mult_gm5p(int mu, const Field_F& w)
{
  Field_F vt, v;

  vt = 0.0;

  assert(mu >= 0);
  assert(mu < m_Ndim);

  mult_up(mu, vt, w);

  // v = (Field_F) mult_gm5(vt);
  mult_gm5(v, vt);

  return v;
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::mult_up(int mu, Field& w, const Field& f)
{
  Field_F vt(f.nvol(), 1);

  for (int ex = 0; ex < f.nex(); ++ex) {
    vt.setpart_ex(0, f, ex);
    shift.backward(trf, f, m_boundary[mu], mu);
    trf2.mult_Field_Gn(0, *m_U, mu, trf, 0);
    vt.mult_GMproj2(-1, m_GM[mu], trf2);
    w.addpart_ex(ex, vt, 0);
  }
}


//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::mult_dn(int mu, Field& w, const Field& f)
{
  Field_F vt(f.nvol(), 1);

  for (int ex = 0; ex < f.nex(); ++ex) {
    trf.mult_Field_Gd(0, *m_U, mu, (Field_F)f, ex);
    shift.forward(trf2, trf, m_boundary[mu], mu);
    vt.mult_GMproj2(1, m_GM[mu], trf2);
    w.addpart_ex(ex, vt, 0);
  }
}


//====================================================================

/*
void Fopr_Wilson::Fopr_Wilson_impl::mult_p(int mu, Field_F& w, const Field_F& f){

  Field_F vt(f.nvol(),1);

  for(int ex = 0; ex < f.nex(); ++ex){
    vt.setpart_ex(0,f,ex);
    shift.backward(trf,f,m_boundary[mu],mu);
    trf2.mult_Field_Gn(0,*m_U,mu,trf,0);
    vt.mult_GMproj2(-1,m_GM[mu],trf2);
    w.addpart_ex(ex,vt,0);
  }

}
//====================================================================
void Fopr_Wilson::Fopr_Wilson_impl::mult_m(int mu, Field_F& w, const Field_F& f){

  Field_F vt(f.nvol(),1);

  for(int ex = 0; ex < f.nex(); ++ex){
    trf.mult_Field_Gd(0,*m_U,mu,f,ex);
    shift.forward(trf2,trf,m_boundary[mu],mu);
    vt.mult_GMproj2( 1,m_GM[mu],trf2);
    w.addpart_ex(ex,vt,0);
  }

}
*/
//====================================================================
//============================================================END=====