field_F_imp.cpp

Go to the documentation of this file.

#include "field_F.h"

using std::valarray;

//#include "threadManager_OpenMP.h"

#if defined USE_GROUP_SU3
#include "field_F_imp_SU3.inc"
#elif defined USE_GROUP_SU2
#include "field_F_imp_SU2.inc"
#elif defined USE_GROUP_SU_N
#include "field_F_imp_SU_N.inc"
#endif

//====================================================================
void Field_F::check()
{
  //  assert(NC == CommonParameters::Nc());
  //  assert(ND == CommonParameters::Nd());

  check_Nc();
}


//====================================================================
void mult_GM(Field_F& y, const GammaMatrix& gm, const Field_F& x)
{
  assert(x.nex() == y.nex());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;

  int    id[Nd];
  int    idc_r[Nd];
  int    idc_i[Nd];
  double gv_r[Nd];
  double gv_i[Nd];

  for (int s = 0; s < Nd; ++s) {
    id[s]    = gm.index(s);
    gv_r[s]  = gm.value_r(s);
    gv_i[s]  = gm.value_i(s);
    idc_r[s] = gm.index_c(s);
    idc_i[s] = 1 - idc_r[s];
  }

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  const int Nex = y.nex();

  for (int ex = 0; ex < Nex; ++ex) {
    double *v = y.ptr(0, is, ex);
    double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex);

    for (int site = 0; site < ns; ++site) {
      int iv = Ncd2 * site;

      for (int s = 0; s < Nd; ++s) {
        int iv2 = s * Nc2 + iv;
        int iw2 = id[s] * Nc2 + iv;
        for (int ic = 0; ic < Nc; ++ic) {
          v[2 * ic + iv2]     = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
          v[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
        }
      }
    }
  }
}


//====================================================================
void mult_iGM(Field_F& y, const GammaMatrix& gm, const Field_F& x)
{
  assert(x.nex() == y.nex());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;

  int    id[Nd];
  int    idc_r[Nd];
  int    idc_i[Nd];
  double gv_r[Nd];
  double gv_i[Nd];

  for (int s = 0; s < Nd; ++s) {
    id[s]    = gm.index(s);
    gv_r[s]  = -gm.value_i(s);
    gv_i[s]  = gm.value_r(s);
    idc_r[s] = 1 - gm.index_c(s);
    idc_i[s] = 1 - idc_r[s];
  }

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  const int Nex = y.nex();

  for (int ex = 0; ex < Nex; ++ex) {
    double *v = y.ptr(0, is, ex);
    double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex);

    for (int site = 0; site < ns; ++site) {
      int iv = Ncd2 * site;

      for (int s = 0; s < Nd; ++s) {
        int iv2 = s * Nc2 + iv;
        int iw2 = id[s] * Nc2 + iv;
        for (int ic = 0; ic < Nc; ++ic) {
          v[2 * ic + iv2]     = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
          v[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
        }
      }
    }
  }
}


//====================================================================
void mult_GMproj(Field_F& y,
                 const int pm, const GammaMatrix& gm,
                 const Field_F& w)
{
  assert(w.nvol() == y.nvol());
  assert(w.nex() == y.nex());

  mult_GM(y, gm, w);

  if (pm == 1) {
    y += w;
    y *= 0.5;
  } else if (pm == -1) {
    y -= w;
    y *= -0.5;
  } else {
    vout.general("%s: wrong pm.\n", __func__);
    abort();
  }
}


//====================================================================
void mult_GMproj2(Field_F& y,
                  const int pm, const GammaMatrix& gm,
                  const Field_F& w)
{
  assert(w.nvol() == y.nvol());
  assert(w.nex() == y.nex());

  mult_GM(y, gm, w);

  if (pm == 1) {
    y += w;
  } else if (pm == -1) {
    y -= w;
    y *= -1;
  } else {
    vout.general("%s: wrong pm.\n", __func__);
    abort();
  }
}


//====================================================================
void mult_Field_Gn(Field_F& y, const int ex,
                   const Field_G& U, const int ex1,
                   const Field_F& x, const int ex2)
{
  assert(ex < y.nex());
  assert(ex1 < U.nex());
  assert(ex2 < x.nex());
  assert(U.nvol() == y.nvol());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;
  const int Ndf  = Nc2 * Nc;

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  double *v = y.ptr(0, is, ex);
  double *g = const_cast<Field_G *>(&U)->ptr(0, is, ex1);
  double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex2);

  for (int site = 0; site < ns; ++site) {
    int ig = Ndf * site;
    int iv = Ncd2 * site;
    for (int s = 0; s < Nd; ++s) {
      for (int ic = 0; ic < Nc; ++ic) {
        int ig2 = ic * Nc2 + ig;
        int iv2 = s * Nc2 + iv;
        v[2 * ic + iv2]     = mult_Gn_r(&g[ig2], &w[iv2], Nc);
        v[2 * ic + 1 + iv2] = mult_Gn_i(&g[ig2], &w[iv2], Nc);
      }
    }
  }
}


//====================================================================
void mult_Field_Gd(Field_F& y, const int ex,
                   const Field_G& U, const int ex1,
                   const Field_F& x, const int ex2)
{
  assert(ex < y.nex());
  assert(ex1 < U.nex());
  assert(ex2 < x.nex());
  assert(U.nvol() == y.nvol());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;
  const int Ndf  = Nc2 * Nc;

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  double *v = y.ptr(0, is, ex);
  double *g = const_cast<Field_G *>(&U)->ptr(0, is, ex1);
  double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex2);

  for (int site = 0; site < ns; ++site) {
    int ig = Ndf * site;
    int iv = Ncd2 * site;
    for (int s = 0; s < Nd; ++s) {
      for (int ic = 0; ic < Nc; ++ic) {
        int ig2 = ic * 2 + ig;
        int iv2 = s * Nc2 + iv;
        v[2 * ic + iv2]     = mult_Gd_r(&g[ig2], &w[iv2], Nc);
        v[2 * ic + 1 + iv2] = mult_Gd_i(&g[ig2], &w[iv2], Nc);
      }
    }
  }
}


//====================================================================
void multadd_Field_Gn(Field_F& y, const int ex,
                      const Field_G& U, const int ex1,
                      const Field_F& x, const int ex2,
                      const double a)
{
  assert(ex < y.nex());
  assert(ex1 < U.nex());
  assert(ex2 < x.nex());
  assert(U.nvol() == y.nvol());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;
  const int Ndf  = Nc2 * Nc;

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  double *v = y.ptr(0, is, ex);
  double *g = const_cast<Field_G *>(&U)->ptr(0, is, ex1);
  double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex2);

  for (int site = 0; site < ns; ++site) {
    int ig = Ndf * site;
    int iv = Ncd2 * site;
    for (int s = 0; s < Nd; ++s) {
      for (int ic = 0; ic < Nc; ++ic) {
        int ig2 = ic * Nc2 + ig;
        int iv2 = s * Nc2 + iv;
        v[2 * ic + iv2]     += a * mult_Gn_r(&g[ig2], &w[iv2], Nc);
        v[2 * ic + 1 + iv2] += a * mult_Gn_i(&g[ig2], &w[iv2], Nc);
      }
    }
  }
}


//====================================================================
void multadd_Field_Gd(Field_F& y, const int ex,
                      const Field_G& U, const int ex1,
                      const Field_F& x, const int ex2,
                      const double a)
{
  assert(ex < y.nex());
  assert(ex1 < U.nex());
  assert(ex2 < x.nex());
  assert(U.nvol() == y.nvol());
  assert(x.nvol() == y.nvol());

  const int Nd   = x.nd();
  const int Nc   = x.nc();
  const int Nc2  = x.nc2();
  const int Ncd2 = Nc2 * Nd;
  const int Ndf  = Nc2 * Nc;

  //  int nth = ThreadManager_OpenMP::get_num_threads();
  //  int ith = ThreadManager_OpenMP::get_thread_id();
  int nth = 1;
  int ith = 0;

  int is = y.nvol() * ith / nth;
  int ns = y.nvol() * (ith + 1) / nth - is;

  double *v = y.ptr(0, is, ex);
  double *g = const_cast<Field_G *>(&U)->ptr(0, is, ex1);
  double *w = const_cast<Field_F *>(&x)->ptr(0, is, ex2);

  for (int site = 0; site < ns; ++site) {
    int ig = Ndf * site;
    int iv = Ncd2 * site;
    for (int s = 0; s < Nd; ++s) {
      for (int ic = 0; ic < Nc; ++ic) {
        int ig2 = ic * 2 + ig;
        int iv2 = s * Nc2 + iv;
        v[2 * ic + iv2]     += a * mult_Gd_r(&g[ig2], &w[iv2], Nc);
        v[2 * ic + 1 + iv2] += a * mult_Gd_i(&g[ig2], &w[iv2], Nc);
      }
    }
  }
}


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