fopr_Wilson_eo.cpp

Go to the documentation of this file.

#include "fopr_Wilson_eo.h"

#ifdef USE_PARAMETERS_FACTORY
#include "parameters_factory.h"
#endif

using std::valarray;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_string("gamma_matrix_type", "NULL");

    param.Register_double("hopping_parameter", 0.0);
    param.Register_int_vector("boundary_condition", std::valarray<int>());

    param.Register_string("verbose_level", "NULL");
  }


#ifdef USE_PARAMETERS_FACTORY
  bool init_param = ParametersFactory::Register("Fopr.Wilson_eo", append_entry);
#endif
}
//- end

//===================================================================
namespace {
  void mult_Field_Gn_eo(Field_F& w, int ex,
                        const Field_G& U, const int ex1,
                        const Field_F& x, const int ex2,
                        const int ieo)
  {
    assert(ex < w.nex());
    assert(ex1 < U.nex());
    assert(ex2 < x.nex());
    assert(U.nvol() == w.nvol() * 2);
    assert(x.nvol() == w.nvol());

    Vec_SU_N vec(w.nc());

    for (int site = 0, nvol = w.nvol(); site < nvol; ++site) {
      for (int s = 0, nd = w.nd(); s < nd; ++s) {
        vec = U.mat(site + ieo * nvol, ex1) * x.vec(s, site, ex2);
        w.set_vec(s, site, ex, vec);
      }
    }
  }


  //====================================================================
  void mult_Field_Gd_eo(Field_F& w, int ex,
                        const Field_G& U, int ex1,
                        const Field_F& x, int ex2,
                        const int ieo)
  {
    assert(ex < w.nex());
    assert(ex1 < U.nex());
    assert(ex2 < x.nex());
    assert(U.nvol() == w.nvol() * 2);
    assert(x.nvol() == w.nvol());

    Vec_SU_N vec(w.nc());

    for (int site = 0, nvol = w.nvol(); site < nvol; ++site) {
      for (int s = 0, nd = w.nd(); s < nd; ++s) {
        vec = U.mat_dag(site + ieo * nvol, ex1) * x.vec(s, site, ex2);
        w.set_vec(s, site, ex, vec);
      }
    }
  }
}

//====================================================================
void Fopr_Wilson_eo::set_parameters(const Parameters& params)
{
  const string str_vlevel = params.get_string("verbose_level");

  m_vl = vout.set_verbose_level(str_vlevel);

  //- fetch and check input parameters
  double        kappa;
  valarray<int> bc;

  int err = 0;
  err += params.fetch_double("hopping_parameter", kappa);
  err += params.fetch_int_vector("boundary_condition", bc);

  if (err) {
    vout.crucial(m_vl, "Fopr_Wilson_eo: fetch error, input parameter not found.\n");
    abort();
  }


  set_parameters(kappa, bc);
}


//====================================================================
void Fopr_Wilson_eo::set_parameters(const double kappa, const std::valarray<int> bc)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of Fopr_Wilson_eo:\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_eo::set_config(Field *U)
{
  int Ndim = CommonParameters::Ndim();
  int Nvol = CommonParameters::Nvol();

  m_index.convertField(*m_Ueo, *U);
}


//====================================================================
void Fopr_Wilson_eo::init(const std::string repr)
{
  m_repr = repr;

  m_boundary.resize(m_Ndim);
  m_GM.resize(m_Ndim + 1);

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

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

  m_Ueo = new Field_G(m_Nvol, m_Ndim);

  m_mult     = &Fopr_Wilson_eo::mult_undef;
  m_mult_dag = &Fopr_Wilson_eo::mult_undef;
  m_preProp  = 0;
  m_postProp = 0;

  delete gmset;
}


//====================================================================
void Fopr_Wilson_eo::prePropD(Field& Be, Field& bo, const Field& b)
{
  int Nin = b.nin();
  int Nex = b.nex();

  Field be(Nin, m_Nvol2, Nex);

  m_index.convertField(be, b, 0);
  m_index.convertField(bo, b, 1);

  Be.reset(Nin, m_Nvol2, Nex);
  Be = be - (Field)Meo(bo, 0);
}


void Fopr_Wilson_eo::postPropD(Field& x, const Field& xe, const Field& bo)
{
  //Field bo();
  //index.convertField(bo, b, 1);
  int Nin = xe.nin();
  int Nex = xe.nex();

  Field xo(Nin, m_Nvol2, Nex);

  xo = bo - (Field)Meo(xe, 1);

  x.reset(Nin, m_Nvol2 * 2, Nex);
  m_index.reverseField(x, xe, 0);
  m_index.reverseField(x, xo, 1);
}


//====================================================================
void Fopr_Wilson_eo::prePropDag(Field& Be, Field& bo, const Field& b)
{
  int Nin = b.nin();
  int Nex = b.nex();

  Field be(Nin, m_Nvol2, Nex);

  m_index.convertField(be, b, 0);
  m_index.convertField(bo, b, 1);

  Be.reset(Nin, m_Nvol2, Nex);
  Be = be - (Field)Mdageo(bo, 0);
}


void Fopr_Wilson_eo::postPropDag(Field& x, const Field& xe, const Field& bo)
{
  //Field bo();
  //index.convertField(bo, b, 1);
  int Nin = xe.nin();
  int Nex = xe.nex();

  Field xo(Nin, m_Nvol2, Nex);

  xo = bo - (Field)Mdageo(xe, 1);

  x.reset(Nin, m_Nvol2 * 2, Nex);
  m_index.reverseField(x, xe, 0);
  m_index.reverseField(x, xo, 1);
}


//====================================================================
const Field_F Fopr_Wilson_eo::Meo(const Field_F& f, const int ieo)
{
  Field_F w(m_Nvol2, f.nex());

  w = 0.0;

  for (int mu = 0; mu < m_Ndim; ++mu) {
    mult_p(mu, w, (Field_F)f, ieo);
    mult_m(mu, w, (Field_F)f, ieo);
  }

  w *= -m_kappa;

  return w;
}


//====================================================================
const Field_F Fopr_Wilson_eo::Mdageo(const Field_F& f, const int ieo)
{
  Field_F w(m_Nvol2, f.nex());
  Field_F v(m_Nvol2, f.nex());

  v = (Field_F)mult_gm5((Field)f);

  w = Meo_gm5(v, ieo);

  return w;
}


//====================================================================
const Field_F Fopr_Wilson_eo::Meo_gm5(const Field_F& f, const int ieo)
{
  Field_F w(m_Nvol2, f.nex());

  w = Meo((Field)f, ieo);


  Field_F v(m_Nvol2, f.nex());
  v = (Field_F)mult_gm5((Field)w);

  return v;
}


//====================================================================
const Field Fopr_Wilson_eo::mult_gm5(const Field& f)
{
  Field_F w = (Field_F)f;

  w.mult_GM(m_GM[4], (Field_F)f);

  return (Field)w;
}


//====================================================================
const Field Fopr_Wilson_eo::gm5p(const int mu, const Field& f)
{
  Field_F w(f.nvol(), f.nex());
  Field_F vt1(f.nvol(), 1);
  Field_F vt2(f.nvol(), 1);

  for (int ex = 0; ex < f.nex(); ++ex) {
    vt1.setpart_ex(0, f, ex);
    shift.backward(vt2, vt1, m_boundary[mu], mu); // vt2 = vt1(x+mu)
    vt1.mult_GMproj2(-1, m_GM[mu], vt2);          // vt1 = (1 - gamma_mu) vt2
    vt2.mult_GM(m_GM[4], vt1);                    // vt2 = gamma_5 vt1
    w.addpart_ex(ex, vt2, 0);
  }

  return (Field)w;
}


//====================================================================
void Fopr_Wilson_eo::mult_p(int mu,
                            Field_F& w, const Field_F& f,
                            const int ieo)
{
  Field_F vt(f.nvol(), 1);

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


//=====================================================================
void Fopr_Wilson_eo::mult_m(int mu,
                            Field_F& w, const Field_F& f,
                            const int ieo)
{
  Field_F vt(f.nvol(), 1);

  for (int ex = 0; ex < f.nex(); ++ex) {
    mult_Field_Gd_eo(trf, 0, *m_Ueo, mu, f, ex, 1 - ieo);
    shift.forward_h(trf2, trf, m_boundary[mu], mu, ieo);
    vt.mult_GMproj2(1, m_GM[mu], trf2);
    w.addpart_ex(ex, vt, 0);
  }
}


//====================================================================
//============================================================END=====