fopr_Wilson_eo_impl.cpp

Go to the documentation of this file.

#include "fopr_Wilson_eo_impl.h"

using std::valarray;

const std::string Fopr_Wilson_eo::Fopr_Wilson_eo_impl::class_name
  = "Fopr_Wilson_eo_impl";

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

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


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::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);

  delete gmset;
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::prePropD(Field& Be, Field& bo, const Field& b)
{
#pragma omp master
  {
    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);
    Field_F v1(m_Nvol2, Nex);
    copy(Be, be);
    Meo(v1, bo, 0);
    axpy(Be, -1.0, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::postPropD(Field& x, const Field& xe, const Field& bo)
{
#pragma omp master
  {
    int Nin = xe.nin();
    int Nex = xe.nex();

    Field xo(Nin, m_Nvol2, Nex);

    //  xo = bo - (Field)Meo(xe, 1);
    Field_F v1(m_Nvol2, Nex);
    copy(xo, bo);
    Meo(v1, xe, 1);
    axpy(xo, -1.0, v1);

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


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::prePropDag(Field& Be, Field& bo, const Field& b)
{
#pragma omp master
  {
    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);
    Field_F v1(m_Nvol2, Nex);
    copy(Be, be);
    Mdageo(v1, bo, 0);
    axpy(Be, -1.0, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::postPropDag(Field& x, const Field& xe, const Field& bo)
{
#pragma omp master
  {
    //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);
    Field_F v1(m_Nvol2, Nex);
    copy(xo, bo);
    Mdageo(v1, xe, 1);
    axpy(xo, -1.0, v1);

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


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::D(Field& v, const Field& f)
{
#pragma omp master
  {
    assert(f.nex() == 1);
    Field_F v1(m_Nvol2, 1), v2(m_Nvol2, 1);

    copy(v, f); //  v  = f;

    //  v -= (Field)Meo(Meo(f, 1), 0);
    Meo(v1, f, 1);
    Meo(v2, v1, 0);
    axpy(v, -1.0, v2);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::Ddag(Field& v, const Field& f)
{
#pragma omp master
  {
    Field_F v1(m_Nvol2, 1), v2(m_Nvol2, 1);

    copy(v, f); //  v  = f;

    //  v -= (Field)Mdageo(Mdageo(f, 1), 0);
    Mdageo(v1, f, 1);
    Mdageo(v2, v1, 0);
    axpy(v, -1.0, v2);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::DdagD(Field& v, const Field& f)
{
#pragma omp master
  {
    Field w(v);
    D(w, f);
    Ddag(v, w);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::DDdag(Field& v, const Field& f)
{
#pragma omp master
  {
    Field w(v);
    Ddag(w, f);
    D(v, w);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::H(Field& v, const Field& f)
{
#pragma omp master
  {
    Field w(v);
    D(w, f);
    mult_gm5(v, w);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::MeoMoe(Field& v, const Field& f)
{
#pragma omp master
  {
    Field_F v1(m_Nvol2, 1), v2(m_Nvol2, 1);

    //  v -= (Field)Meo(Meo(f, 1), 0);
    Meo(v1, f, 1);
    Meo(v2, v1, 0);
    axpy(v, -1.0, v2);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::Meo(Field& w,
                                              const Field& f, const int ieo)
{
#pragma omp master
  {
    Field_F w2(m_Nvol2, 1), f2(m_Nvol2, 1);
    copy(f2, f);
    w2.set(0.0); // w2 = 0.0

    for (int mu = 0; mu < m_Ndim; ++mu) {
      mult_p(mu, w2, f2, ieo);
      mult_m(mu, w2, f2, ieo);
    }

    scal(w2, -m_kappa); //  w *= -m_kappa;

    copy(w, w2);
  }
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::Mdageo(Field& w,
                                                 const Field& f, const int ieo)
{
#pragma omp master
  {
    Field_F v(m_Nvol2, f.nex());

    mult_gm5(v, f);
    Meo_gm5(w, v, ieo);
  }
}


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

  Meo(w, f, ieo);
  mult_gm5(v, w);
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::mult_gm5(Field& w, const Field& f)
{
#pragma omp master
  {
    Field_F w2 = (Field_F)f;
    mult_GM(w2, m_GM[4], (Field_F)f);
    w = w2;
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::gm5p(const int mu, Field& w,
                                               const Field& f)
{
  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)
    mult_GMproj2(vt1, -1, m_GM[mu], vt2);         // vt1 = (1 - gamma_mu) vt2
    mult_GM(vt2, m_GM[4], vt1);                   // vt2 = gamma_5 vt1
    w.addpart_ex(ex, vt2, 0);
  }
}


//====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::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);
    mult_GMproj2(vt, -1, m_GM[mu], trf2);
    w.addpart_ex(ex, vt, 0);
  }
}


//=====================================================================
void Fopr_Wilson_eo::Fopr_Wilson_eo_impl::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);
    mult_GMproj2(vt, 1, m_GM[mu], trf2);
    w.addpart_ex(ex, vt, 0);
  }
}


//====================================================================
double Fopr_Wilson_eo::Fopr_Wilson_eo_impl::flop_count()
{
  // this counting is based on the Org-implementation of ver.1.2.0.
  // Flop count of mult_GMproj2 is different for upward and downward
  // directions due to the implemetation in Field_F.cpp.
  // Present counting is based on rev.1130. [10 Sep 2014 H.Matsufuru]

  int Nc   = m_Nc;
  int Nd   = m_Nd;
  int Lvol = CommonParameters::Lvol();

  int flop_Meo = Nc * Nd * 2 * 8 * (4 * Nc - 1);      // #(mult_Field_Gn/d)

  flop_Meo += Nc * Nd * 2 * (4 * 3 + 4 * 2);          // #(mult_GMproj2)
  flop_Meo += Nc * Nd * 2 * 8;                        // #(addpart_ex)
  flop_Meo += Nc * Nd * 2;                            // #(scal(kappa))

  int flop_per_site = 2 * flop_Meo + Nc * Nd * 2 * 2; // #(2*Meo + axpy)

  double flop = static_cast<double>(flop_per_site) *
                static_cast<double>(Lvol / 2);

  if ((m_mode == "DdagD") || (m_mode == "DDdag")) {
    flop *= 2.0;
  }

  return flop;
}


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