fopr_Wilson_General_impl.cpp

Go to the documentation of this file.

#include "fopr_Wilson_General_impl.h"

#if defined USE_GROUP_SU3
#include "fopr_Wilson_impl_SU3.inc"
#elif defined USE_GROUP_SU2
#include "fopr_Wilson_impl_SU2.inc"
#elif defined USE_GROUP_SU_N
#include "fopr_Wilson_impl_SU_N.inc"
#endif

const std::string Fopr_Wilson_General::Fopr_Wilson_General_impl::class_name = "Fopr_Wilson_General_impl";

//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::init(std::string repr)
{
  m_vl = CommonParameters::Vlevel();

  vout.general(m_vl, "%s: Construction of fermion operator(imp).\n", class_name.c_str());

  check_Nc();

  m_Nc  = CommonParameters::Nc();
  m_Nd  = CommonParameters::Nd();
  m_Nvc = 2 * m_Nc;
  m_Ndf = 2 * m_Nc * m_Nc;

  m_Nx = CommonParameters::Nx();
  m_Ny = CommonParameters::Ny();
  m_Nz = CommonParameters::Nz();
  m_Nt = CommonParameters::Nt();

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

  m_repr = repr;

  m_GM.resize(m_Ndim + 1);

  unique_ptr<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_U = 0;

  m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_undef;
  m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_undef;

  if (m_repr == "Dirac") {
    m_D            = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D_dirac;
    m_gm5          = &Fopr_Wilson_General::Fopr_Wilson_General_impl::gm5_dirac;
    m_mult_t_plus  = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_plus_dirac;
    m_mult_t_minus = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_minus_dirac;
    m_D_ex         = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D_ex_dirac;
  } else if (m_repr == "Chiral") {
    m_D            = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D_chiral;
    m_gm5          = &Fopr_Wilson_General::Fopr_Wilson_General_impl::gm5_chiral;
    m_mult_t_plus  = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_plus_chiral;
    m_mult_t_minus = &Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_minus_chiral;
    m_D_ex         = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D_ex_chiral;
  } else {
    vout.crucial(m_vl, "%s: input repr is undefined.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  m_w1.reset(m_Nvc * m_Nd, m_Nvol, 1);
  m_w2.reset(m_Nvc * m_Nd, m_Nvol, 1);

  int Nvx = 2 * m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
  vcp1_x_plus  = new double[Nvx];
  vcp2_x_plus  = new double[Nvx];
  vcp1_x_minus = new double[Nvx];
  vcp2_x_minus = new double[Nvx];

  int Nvy = 2 * m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
  vcp1_y_plus  = new double[Nvy];
  vcp2_y_plus  = new double[Nvy];
  vcp1_y_minus = new double[Nvy];
  vcp2_y_minus = new double[Nvy];

  int Nvz = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
  vcp1_z_plus  = new double[Nvz];
  vcp2_z_plus  = new double[Nvz];
  vcp1_z_minus = new double[Nvz];
  vcp2_z_minus = new double[Nvz];

  // int Nvt = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
  int Nvt = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
  vcp1_t_plus  = new double[Nvt];
  vcp2_t_plus  = new double[Nvt];
  vcp1_t_minus = new double[Nvt];
  vcp2_t_minus = new double[Nvt];

  setup_thread();
}


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

  if (m_mode == "D") {
    m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D;
    m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::Ddag;
  } else if (m_mode == "Ddag") {
    m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::Ddag;
    m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::D;
  } else if (m_mode == "DdagD") {
    m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::DdagD;
    m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::DdagD;
  } else if (m_mode == "DDdag") {
    m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::DDdag;
    m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::DDdag;
  } else if (m_mode == "H") {
    m_mult     = &Fopr_Wilson_General::Fopr_Wilson_General_impl::H;
    m_mult_dag = &Fopr_Wilson_General::Fopr_Wilson_General_impl::H;
  } else {
    vout.crucial(m_vl, "%s: input mode is undefined.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::tidyup()
{
  delete[]  vcp1_x_plus;
  delete[]  vcp2_x_plus;
  delete[]  vcp1_x_minus;
  delete[]  vcp2_x_minus;

  delete[]  vcp1_y_plus;
  delete[]  vcp2_y_plus;
  delete[]  vcp1_y_minus;
  delete[]  vcp2_y_minus;

  delete[]  vcp1_z_plus;
  delete[]  vcp2_z_plus;
  delete[]  vcp1_z_minus;
  delete[]  vcp2_z_minus;

  delete[]  vcp1_t_plus;
  delete[]  vcp2_t_plus;
  delete[]  vcp1_t_minus;
  delete[]  vcp2_t_minus;
}


//====================================================================
std::string Fopr_Wilson_General::Fopr_Wilson_General_impl::get_mode() const
{
  return m_mode;
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::set_parameters(const double kappa_s,
                                                                   const double kappa_t,
                                                                   const double nu_s,
                                                                   const double r_s,
                                                                   const std::vector<int> bc)
{
  assert(bc.size() == m_Ndim);

  //- print input parameters
  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
  vout.general(m_vl, "  kappa_s  = %12.8f\n", kappa_s);
  vout.general(m_vl, "  kappa_t  = %12.8f\n", kappa_t);
  vout.general(m_vl, "  nu_s     = %12.8f\n", nu_s);
  vout.general(m_vl, "  r_s      = %12.8f\n", r_s);
  for (int mu = 0; mu < m_Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

  //- store values
  m_kappa_s = kappa_s;
  m_kappa_t = kappa_t;
  m_nu_s    = nu_s;
  m_r_s     = r_s;

  for (int mu = 0; mu < m_Ndim; ++mu) {
    m_boundary[mu] = bc[mu];
  }

  //- additional boundary condition for each node:
  for (int idir = 0; idir < m_Ndim; ++idir) {
    m_boundary2[idir] = 1.0;
    if (Communicator::ipe(idir) == 0) m_boundary2[idir] = m_boundary[idir];
  }
}


//====================================================================
double Fopr_Wilson_General::Fopr_Wilson_General_impl::flop_count()
{
  // The following counting explicitly depends on the implementation.
  // It will be recalculated when the code is modified.
  // The present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]

  int    Lvol = CommonParameters::Lvol();
  double flop_site, flop;

  if (m_repr == "Dirac") {
    flop_site = static_cast<double>(
      m_Nc * m_Nd * (4 + 6 * (4 * m_Nc + 2) + 2 * (4 * m_Nc + 1)));
  } else if (m_repr == "Chiral") {
    flop_site = static_cast<double>(
      m_Nc * m_Nd * (4 + 8 * (4 * m_Nc + 2)));
  } else {
    vout.crucial(m_vl, "%s: input repr is undefined.\n",
                 class_name.c_str());
    exit(EXIT_FAILURE);
  }

  flop = flop_site * static_cast<double>(Lvol);

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

  return flop;
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::D_dirac(Field& w, const Field& f)
{
  // D_ex_dirac(w, 0, f, 0);
  // daypx(w, -m_kappa_s, f);

  clear(w);

  clear(m_w1);
  for (int mu = 0; mu < (m_Ndim - 1); ++mu) {
    mult_up(mu, m_w1, f);
    mult_dn(mu, m_w1, f);
  }
  daxpy(w, -m_kappa_s, m_w1);

  clear(m_w1);
  int mu = (m_Ndim - 1);
  {
    mult_up(mu, m_w1, f);
    mult_dn(mu, m_w1, f);
  }
  daxpy(w, -m_kappa_t, m_w1);


  daxpy(w, 1.0, f); // w += f;
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::D_chiral(Field& w, const Field& f)
{
  // D_ex_chiral(w, 0, f, 0);
  // daypx(w, -m_kappa_s, f);

  clear(w);

  clear(m_w1);
  for (int mu = 0; mu < (m_Ndim - 1); ++mu) {
    mult_up(mu, m_w1, f);
    mult_dn(mu, m_w1, f);
  }
  daxpy(w, -m_kappa_s, m_w1);

  clear(m_w1);
  int mu = (m_Ndim - 1);
  {
    mult_up(mu, m_w1, f);
    mult_dn(mu, m_w1, f);
  }
  daxpy(w, -m_kappa_t, m_w1);


  daxpy(w, 1.0, f); // w += f;
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_up(int mu,
                                                            Field& w, const Field& f)
{
  if (mu == 0) {
    mult_x_plus(w, f);
  } else if (mu == 1) {
    mult_y_plus(w, f);
  } else if (mu == 2) {
    mult_z_plus(w, f);
  } else if (mu == 3) {
    (this->*m_mult_t_plus)(w, f);
  } else {
    vout.crucial(m_vl, "%s: illegal direction of mu in mult_up.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_dn(int mu, Field& w, const Field& f)
{
  if (mu == 0) {
    mult_x_minus(w, f);
  } else if (mu == 1) {
    mult_y_minus(w, f);
  } else if (mu == 2) {
    mult_z_minus(w, f);
  } else if (mu == 3) {
    (this->*m_mult_t_minus)(w, f);
  } else {
    vout.crucial(m_vl, "%s: illegal direction of mu in mult_dn.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::D_ex_dirac(Field& w, const int ex1,
                                                               const Field& f, const int ex2)
{
  int          Ninvol = m_Nvc * m_Nd * m_Nvol;
  const double *v1    = f.ptr(Ninvol * ex2);
  double       *v2    = w.ptr(Ninvol * ex1);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus1_thread(i, vcp1_x_plus, v1);
    mult_x_minus1_thread(i, vcp1_x_minus, v1);
    mult_y_plus1_thread(i, vcp1_y_plus, v1);
    mult_y_minus1_thread(i, vcp1_y_minus, v1);
    mult_z_plus1_thread(i, vcp1_z_plus, v1);
    mult_z_minus1_thread(i, vcp1_z_minus, v1);
    mult_t_plus1_dirac_thread(i, vcp1_t_plus, v1);
    mult_t_minus1_dirac_thread(i, vcp1_t_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_x = 0;
    int Nvx    = 2 * m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
    Communicator::exchange(Nvx, vcp2_x_plus, vcp1_x_plus, idir_x, 1, 1);
    Communicator::exchange(Nvx, vcp2_x_minus, vcp1_x_minus, idir_x, -1, 2);

    int idir_y = 1;
    int Nvy    = 2 * m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
    Communicator::exchange(Nvy, vcp2_y_plus, vcp1_y_plus, idir_y, 1, 3);
    Communicator::exchange(Nvy, vcp2_y_minus, vcp1_y_minus, idir_y, -1, 4);

    int idir_z = 2;
    int Nvz    = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
    Communicator::exchange(Nvz, vcp2_z_plus, vcp1_z_plus, idir_z, 1, 5);
    Communicator::exchange(Nvz, vcp2_z_minus, vcp1_z_minus, idir_z, -1, 6);

    int idir_t = 3;
    // int Nvt = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nvt = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nvt, vcp2_t_plus, vcp1_t_plus, idir_t, 1, 7);
    Communicator::exchange(Nvt, vcp2_t_minus, vcp1_t_minus, idir_t, -1, 8);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    clear_thread(i, v2);
    mult_x_plus_bulk_thread(i, v2, v1);
    mult_x_minus_bulk_thread(i, v2, v1);
    mult_y_plus_bulk_thread(i, v2, v1);
    mult_y_minus_bulk_thread(i, v2, v1);
    mult_z_plus_bulk_thread(i, v2, v1);
    mult_z_minus_bulk_thread(i, v2, v1);
    mult_t_plus_bulk_dirac_thread(i, v2, v1);
    mult_t_minus_bulk_dirac_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus2_thread(i, v2, vcp2_x_plus);
    mult_x_minus2_thread(i, v2, vcp2_x_minus);
    mult_y_plus2_thread(i, v2, vcp2_y_plus);
    mult_y_minus2_thread(i, v2, vcp2_y_minus);
    mult_z_plus2_thread(i, v2, vcp2_z_plus);
    mult_z_minus2_thread(i, v2, vcp2_z_minus);
    mult_t_plus2_dirac_thread(i, v2, vcp2_t_plus);
    mult_t_minus2_dirac_thread(i, v2, vcp2_t_minus);
  }

  // for (int i = is; i < is + ns; ++i) {
  // daypx_thread(i, v2, -m_kappa, v1);   // w = -m_kappa * w + f.
  //   daypx_thread(i, v2, -m_kappa_s, v1);
  // }


  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::D_ex_chiral(Field& w, const int ex1,
                                                                const Field& f, const int ex2)
{
  int          Ninvol = m_Nvc * m_Nd * m_Nvol;
  const double *v1    = f.ptr(Ninvol * ex2);
  double       *v2    = w.ptr(Ninvol * ex1);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus1_thread(i, vcp1_x_plus, v1);
    mult_x_minus1_thread(i, vcp1_x_minus, v1);
    mult_y_plus1_thread(i, vcp1_y_plus, v1);
    mult_y_minus1_thread(i, vcp1_y_minus, v1);
    mult_z_plus1_thread(i, vcp1_z_plus, v1);
    mult_z_minus1_thread(i, vcp1_z_minus, v1);
    mult_t_plus1_chiral_thread(i, vcp1_t_plus, v1);
    mult_t_minus1_chiral_thread(i, vcp1_t_minus, v1);
  }
#pragma omp barrier

#pragma omp master
  {
    int idir_x = 0;
    int Nvx    = 2 * m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
    Communicator::exchange(Nvx, vcp2_x_plus, vcp1_x_plus, idir_x, 1, 1);
    Communicator::exchange(Nvx, vcp2_x_minus, vcp1_x_minus, idir_x, -1, 2);

    int idir_y = 1;
    int Nvy    = 2 * m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
    Communicator::exchange(Nvy, vcp2_y_plus, vcp1_y_plus, idir_y, 1, 3);
    Communicator::exchange(Nvy, vcp2_y_minus, vcp1_y_minus, idir_y, -1, 4);

    int idir_z = 2;
    int Nvz    = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
    Communicator::exchange(Nvz, vcp2_z_plus, vcp1_z_plus, idir_z, 1, 5);
    Communicator::exchange(Nvz, vcp2_z_minus, vcp1_z_minus, idir_z, -1, 6);

    int idir_t = 3;
    // int Nvt = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nvt = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nvt, vcp2_t_plus, vcp1_t_plus, idir_t, 1, 7);
    Communicator::exchange(Nvt, vcp2_t_minus, vcp1_t_minus, idir_t, -1, 8);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    clear_thread(i, v2);
    mult_x_plus_bulk_thread(i, v2, v1);
    mult_x_minus_bulk_thread(i, v2, v1);
    mult_y_plus_bulk_thread(i, v2, v1);
    mult_y_minus_bulk_thread(i, v2, v1);
    mult_z_plus_bulk_thread(i, v2, v1);
    mult_z_minus_bulk_thread(i, v2, v1);
    mult_t_plus_bulk_chiral_thread(i, v2, v1);
    mult_t_minus_bulk_chiral_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus2_thread(i, v2, vcp2_x_plus);
    mult_x_minus2_thread(i, v2, vcp2_x_minus);
    mult_y_plus2_thread(i, v2, vcp2_y_plus);
    mult_y_minus2_thread(i, v2, vcp2_y_minus);
    mult_z_plus2_thread(i, v2, vcp2_z_plus);
    mult_z_minus2_thread(i, v2, vcp2_z_minus);
    mult_t_plus2_chiral_thread(i, v2, vcp2_t_plus);
    mult_t_minus2_chiral_thread(i, v2, vcp2_t_minus);
  }

  // for (int i = is; i < is + ns; ++i) {
  // daypx_thread(i, v2, -m_kappa, v1);   // w = -m_kappa * w + f.
  //  daypx_thread(i, v2, -m_kappa_s, v1);   // w = -m_kappa * w + f.
  // }


  ThreadManager_OpenMP::sync_barrier_all();

  // clear(w);
  // mult_x_plus(w,f);
  // mult_x_minus(w,f);
  // mult_y_plus(w,f);
  // mult_y_minus(w,f);
  // mult_z_plus(w,f);
  // mult_z_minus(w,f);
  // mult_t_plus_chiral(w,f);
  // mult_t_minus_chiral(w,f);
  // daypx(w, -m_kappa, f);   // w = -m_kappa * w + f.
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_gm5p(int mu,
                                                              Field& v, const Field& w)
{
  clear(m_w2);

  mult_up(mu, m_w2, w);
  mult_gm5(v, m_w2);
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::proj_chiral(Field& w, const int ex1,
                                                                const Field& v, const int ex2,
                                                                const int ipm)
{
  double fpm = 0.0;

  if (ipm == 1) {
    fpm = 1.0;
  } else if (ipm == -1) {
    fpm = -1.0;
  } else {
    vout.crucial(m_vl, "%s: illegal chirality = %d.\n", class_name.c_str(), ipm);
    exit(EXIT_FAILURE);
  }

  m_w1.setpart_ex(0, v, ex2);
  mult_gm5(m_w2, m_w1);
  m_w1.addpart_ex(0, m_w2, 0, fpm);
  w.setpart_ex(ex1, m_w1, 0);

#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::daxpy(Field& w,
                                                          double fac, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    daxpy_thread(i, v2, fac, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::daypx(Field& w,
                                                          double fac, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    daypx_thread(i, v2, fac, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::scal(Field& w,
                                                         double fac)
{
  double *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    scal_thread(i, v2, fac);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::clear(Field& w)
{
  double *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    clear_thread(i, v2);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::gm5_dirac(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    gm5_dirac_thread(i, v2, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::gm5_chiral(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

  for (int i = is; i < is + ns; ++i) {
    gm5_chiral_thread(i, v2, v1);
  }
#pragma omp barrier
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_x_plus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus1_thread(i, vcp1_x_plus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_x = 0;
    int Nv     = 2 * m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
    Communicator::exchange(Nv, vcp2_x_plus, vcp1_x_plus, idir_x, 1, 1);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_x_plus2_thread(i, v2, vcp2_x_plus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_x_minus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_x_minus1_thread(i, vcp1_x_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_x = 0;
    int Nv     = 2 * m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
    Communicator::exchange(Nv, vcp2_x_minus, vcp1_x_minus, idir_x, -1, 2);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_x_minus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_x_minus2_thread(i, v2, vcp2_x_minus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_y_plus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_y_plus1_thread(i, vcp1_y_plus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_y = 1;
    int Nv     = 2 * m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
    Communicator::exchange(Nv, vcp2_y_plus, vcp1_y_plus, idir_y, 1, 3);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_y_plus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_y_plus2_thread(i, v2, vcp2_y_plus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_y_minus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_y_minus1_thread(i, vcp1_y_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_y = 1;
    int Nv     = 2 * m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
    Communicator::exchange(Nv, vcp2_y_minus, vcp1_y_minus, idir_y, -1, 4);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_y_minus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_y_minus2_thread(i, v2, vcp2_y_minus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_z_plus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_z_plus1_thread(i, vcp1_z_plus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_z = 2;
    int Nv     = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
    Communicator::exchange(Nv, vcp2_z_plus, vcp1_z_plus, idir_z, 1, 5);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_z_plus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_z_plus2_thread(i, v2, vcp2_z_plus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_z_minus(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_z_minus1_thread(i, vcp1_z_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_z = 2;
    int Nv     = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
    Communicator::exchange(Nv, vcp2_z_minus, vcp1_z_minus, idir_z, -1, 6);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_z_minus_bulk_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_z_minus2_thread(i, v2, vcp2_z_minus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_plus_dirac(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus1_dirac_thread(i, vcp1_t_plus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_t = 3;
    // int Nv = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nv, vcp2_t_plus, vcp1_t_plus, idir_t, 1, 7);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus_bulk_dirac_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus2_dirac_thread(i, v2, vcp2_t_plus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_minus_dirac(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus1_dirac_thread(i, vcp1_t_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_t = 3;
    // int Nv = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nv, vcp2_t_minus, vcp1_t_minus, idir_t, -1, 8);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus_bulk_dirac_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus2_dirac_thread(i, v2, vcp2_t_minus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_plus_chiral(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus1_chiral_thread(i, vcp1_t_plus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_t = 3;
    // int Nv = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nv, vcp2_t_plus, vcp1_t_plus, idir_t, 1, 7);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus_bulk_chiral_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_t_plus2_chiral_thread(i, v2, vcp2_t_plus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


//====================================================================
void Fopr_Wilson_General::Fopr_Wilson_General_impl::mult_t_minus_chiral(Field& w, const Field& f)
{
  const double *v1 = f.ptr(0);
  double       *v2 = w.ptr(0);

  int Nthread  = ThreadManager_OpenMP::get_num_threads();
  int i_thread = ThreadManager_OpenMP::get_thread_id();

  int is = m_Ntask * i_thread / Nthread;
  int ns = m_Ntask * (i_thread + 1) / Nthread - is;

#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus1_chiral_thread(i, vcp1_t_minus, v1);
  }

#pragma omp barrier
#pragma omp master
  {
    int idir_t = 3;
    // int Nv = 2 * m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nz;
    Communicator::exchange(Nv, vcp2_t_minus, vcp1_t_minus, idir_t, -1, 8);
  }
#pragma omp barrier

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus_bulk_chiral_thread(i, v2, v1);
  }

  for (int i = is; i < is + ns; ++i) {
    mult_t_minus2_chiral_thread(i, v2, vcp2_t_minus);
  }

  ThreadManager_OpenMP::sync_barrier_all();
}


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