projection_Stout_SU3.cpp

Go to the documentation of this file.

#include "projection_Stout_SU3.h"

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

#ifdef USE_FACTORY
namespace {
  Projection *create_object()
  {
    return new Projection_Stout_SU3();
  }


  bool init = Projection::Factory::Register("Stout_SU3", create_object);
}
#endif

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


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

//- parameters class
Parameters_Projection_Stout_SU3::Parameters_Projection_Stout_SU3() { append_entry(*this); }
//- end

const std::string Projection_Stout_SU3::class_name = "Projection_Stout_SU3";

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

  m_vl = vout.set_verbose_level(str_vlevel);
}


//====================================================================
void Projection_Stout_SU3::setup()
{
  m_flop = 0;
  m_time = 0.0;

  assert(CommonParameters::Nc() == NC);
}


//====================================================================
void Projection_Stout_SU3::print_stat()
{
  double gflops = 1.e-9 * double(m_flop) / m_time;

  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
  vout.general(m_vl, "  total time: %f\n", m_time);
  vout.general(m_vl, "  total flop: %d\n", m_flop);
  vout.general(m_vl, "  GFlops    : %f\n", gflops);
}


//====================================================================
void Projection_Stout_SU3::project(Field_G& U,
                                   double alpha,
                                   const Field_G& Cst, const Field_G& Uorg)
{
  //  int id = 31;
  //  KEK_FopCountStart(id);
  double time0 = Communicator::get_time();

  // in stout projection, parameter alpha is dummy.

  int Nex  = Uorg.nex();
  int Nvol = Uorg.nvol();

  assert(Cst.nex() == Nex);
  assert(Cst.nvol() == Nvol);
  assert(U.nex() == Nex);
  assert(U.nvol() == Nvol);

  int NinG = Uorg.nin();

  Mat_SU_N iQ0(NC), iQ1(NC), iQ2(NC), iQ3(NC);
  Mat_SU_N ct(NC), ut(NC), ut2(NC), e_iQ(NC);
  iQ0.unit();

  double   u, w;
  dcomplex f0, f1, f2, qt;

  for (int mu = 0; mu < Nex; ++mu) {
    for (int site = 0; site < Nvol; ++site) {
      Uorg.mat(ut, site, mu);
      Cst.mat(ct, site, mu);
      iQ1.mult_nd(ct, ut);
      iQ1.at();
      iQ2.mult_nn(iQ1, iQ1);
      iQ3.mult_nn(iQ1, iQ2);

      double norm = iQ1.norm2();
      if (norm > 1.0e-10) {
        set_uw(u, w, iQ2, iQ3);
        set_fj(f0, f1, f2, u, w);

        for (int cc = 0; cc < NC * NC; ++cc) {
          qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
               + f1 *cmplx(iQ1.i(cc), -iQ1.r(cc))
               - f2 *cmplx(iQ2.r(cc), iQ2.i(cc));
          e_iQ.setr(cc, real(qt));
          e_iQ.seti(cc, imag(qt));
        }
      } else {
        //  vout.general(m_vl,"project: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
        e_iQ.unit();
      }

      ut2.mult_nn(e_iQ, ut);
      U.set_mat(site, mu, ut2);
    }
  }

  /*
  unsigned long count;
  double time;
  KEK_FopCountFinish(id,&count,&time);
  m_time += time;
  m_flop += count;
  */
  double time1 = Communicator::get_time();
  m_time += time1 - time0;
}


//====================================================================
void Projection_Stout_SU3::exp_iQ(Field_G& e_iQ, const Field_G& iQ)
{
  int Nvol = iQ.nvol();
  int Nex  = iQ.nex();

  Mat_SU_N iQ0(NC), iQ1(NC), iQ2(NC), iQ3(NC);

  iQ0.unit();

  double   u, w;
  dcomplex f0, f1, f2, qt;

  for (int mu = 0; mu < Nex; ++mu) {
    for (int site = 0; site < Nvol; ++site) {
      iQ1 = iQ.mat(site, mu);
      iQ2 = iQ1 * iQ1;
      iQ3 = iQ1 * iQ2;

      double norm = iQ1.norm2();
      if (norm > 1.0e-10) {
        set_uw(u, w, iQ2, iQ3);
        set_fj(f0, f1, f2, u, w);

        for (int cc = 0; cc < NC * NC; ++cc) {
          qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
               + f1 *cmplx(iQ1.i(cc), -iQ1.r(cc))
               - f2 *cmplx(iQ2.r(cc), iQ2.i(cc));
          e_iQ.set_ri(cc, site, mu, real(qt), imag(qt));
        }
      } else {
        //      vout.general(m_vl,"exp_iQ: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
        e_iQ.set_mat(site, mu, iQ0);
      }
    }
  }
}


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

//====================================================================
void Projection_Stout_SU3::force_recursive(Field_G& Xi, Field_G& iTheta,
                                           double alpha, const Field_G& Sigmap,
                                           const Field_G& Cst, const Field_G& Uorg)
{
  // in stout projection, parameter alpha is dummy.

  //  int id = 31;
  //  KEK_FopCountStart(id);
  double time0 = Communicator::get_time();

  int Nvol = CommonParameters::Nvol();
  // int NinG = 2 * NC * NC;
  int Nex  = Xi.nex();

  assert(Xi.nvol() == Nvol);
  assert(iTheta.nvol() == Nvol);
  assert(Sigmap.nvol() == Nvol);
  assert(Cst.nvol() == Nvol);
  assert(Uorg.nvol() == Nvol);
  assert(iTheta.nex() == Nex);
  assert(Sigmap.nex() == Nex);
  assert(Cst.nex() == Nex);
  assert(Uorg.nex() == Nex);

  Mat_SU_N C_tmp(NC), U_tmp(NC), Sigmap_tmp(NC);
  Mat_SU_N iQ0(NC), iQ1(NC), iQ2(NC), iQ3(NC), e_iQ(NC);
  Mat_SU_N B1(NC), B2(NC);
  Mat_SU_N USigmap(NC), iQUS(NC), iUSQ(NC), iGamma(NC);
  Mat_SU_N Xi_tmp(NC), iTheta_tmp(NC);
  Mat_SU_N tmp1(NC), tmp2(NC);
  iQ0.unit();

  double   u, w, u2, w2, cos_w, xi0, xi1, fden;
  dcomplex f0, f1, f2, h0, h1, h2, e2iu, emiu, qt;
  dcomplex r01, r11, r21, r02, r12, r22, tr1, tr2;
  dcomplex b10, b11, b12, b20, b21, b22;

  for (int mu = 0; mu < Nex; ++mu) {
    for (int site = 0; site < Nvol; ++site) {
      Cst.mat(C_tmp, site, mu);
      Uorg.mat(U_tmp, site, mu);
      // Sigmap_tmp \f$=\Sigma_\mu'(x)\f$
      Sigmap.mat(Sigmap_tmp, site, mu);

      iQ1.mult_nd(C_tmp, U_tmp);
      iQ1.at();
      iQ2.mult_nn(iQ1, iQ1);
      iQ3.mult_nn(iQ1, iQ2);

      // In order to aviod 1Q1=0
      double norm = iQ1.norm2();
      if (norm > 1.0e-10) {
        set_uw(u, w, iQ2, iQ3);
        set_fj(f0, f1, f2, u, w);

        for (int cc = 0; cc < NC * NC; ++cc) {
          qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
               + f1 *cmplx(iQ1.i(cc), -iQ1.r(cc))
               - f2 *cmplx(iQ2.r(cc), iQ2.i(cc));
          e_iQ.set(cc, real(qt), imag(qt));
        }

        xi0   = func_xi0(w);
        xi1   = func_xi1(w);
        u2    = u * u;
        w2    = w * w;
        cos_w = cos(w);

        emiu = cmplx(cos(u), -sin(u));
        e2iu = cmplx(cos(2.0 * u), sin(2.0 * u));

        r01 = cmplx(2.0 * u, 2.0 * (u2 - w2)) * e2iu
              + emiu *cmplx(16.0 * u * cos_w + 2.0 * u * (3.0 * u2 + w2) * xi0,
                            -8.0 * u2 * cos_w + 2.0 * (9.0 * u2 + w2) * xi0);

        r11 = cmplx(2.0, 4.0 * u) * e2iu
              + emiu *cmplx(-2.0 * cos_w + (3.0 * u2 - w2) * xi0,
                            2.0 * u * cos_w + 6.0 * u * xi0);

        r21 = cmplx(0.0, 2.0) * e2iu
              + emiu *cmplx(-3.0 * u * xi0, cos_w - 3.0 * xi0);

        r02 = cmplx(-2.0, 0.0) * e2iu
              + emiu *cmplx(-8.0 * u2 * xi0,
                            2.0 * u * (cos_w + xi0 + 3.0 * u2 * xi1));

        r12 = emiu * cmplx(2.0 * u * xi0,
                           -cos_w - xi0 + 3.0 * u2 * xi1);

        r22 = emiu * cmplx(xi0, -3.0 * u * xi1);

        fden = 1.0 / (2 * (9.0 * u2 - w2) * (9.0 * u2 - w2));

        b10 = cmplx(2.0 * u, 0.0) * r01 + cmplx(3.0 * u2 - w2, 0.0) * r02
              - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f0;

        b11 = cmplx(2.0 * u, 0.0) * r11 + cmplx(3.0 * u2 - w2, 0.0) * r12
              - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f1;

        b12 = cmplx(2.0 * u, 0.0) * r21 + cmplx(3.0 * u2 - w2, 0.0) * r22
              - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f2;

        b20 = r01 - cmplx(3.0 * u, 0.0) * r02 - cmplx(24.0 * u, 0.0) * f0;

        b21 = r11 - cmplx(3.0 * u, 0.0) * r12 - cmplx(24.0 * u, 0.0) * f1;

        b22 = r21 - cmplx(3.0 * u, 0.0) * r22 - cmplx(24.0 * u, 0.0) * f2;

        b10 *= cmplx(fden, 0.0);
        b11 *= cmplx(fden, 0.0);
        b12 *= cmplx(fden, 0.0);
        b20 *= cmplx(fden, 0.0);
        b21 *= cmplx(fden, 0.0);
        b22 *= cmplx(fden, 0.0);

        for (int cc = 0; cc < NC * NC; ++cc) {
          qt = b10 * cmplx(iQ0.r(cc), iQ0.i(cc))
               + b11 *cmplx(iQ1.i(cc), -iQ1.r(cc))
               - b12 *cmplx(iQ2.r(cc), iQ2.i(cc));
          B1.set(cc, real(qt), imag(qt));
          qt = b20 * cmplx(iQ0.r(cc), iQ0.i(cc))
               + b21 *cmplx(iQ1.i(cc), -iQ1.r(cc))
               - b22 *cmplx(iQ2.r(cc), iQ2.i(cc));
          B2.set(cc, real(qt), imag(qt));
        }

        USigmap.mult_nn(U_tmp, Sigmap_tmp);

        tmp1.mult_nn(USigmap, B1);
        tmp2.mult_nn(USigmap, B2);
        tr1 = cmplx(tmp1.r(0) + tmp1.r(4) + tmp1.r(8),
                    tmp1.i(0) + tmp1.i(4) + tmp1.i(8));
        tr2 = cmplx(tmp2.r(0) + tmp2.r(4) + tmp2.r(8),
                    tmp2.i(0) + tmp2.i(4) + tmp2.i(8));

        iQUS.mult_nn(iQ1, USigmap);
        iUSQ.mult_nn(USigmap, iQ1);

        for (int cc = 0; cc < NC * NC; ++cc) {
          qt = tr1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
               - tr2 *cmplx(iQ2.r(cc), iQ2.i(cc))
               + f1  *cmplx(USigmap.r(cc), USigmap.i(cc))
               + f2  *cmplx(iQUS.i(cc), -iQUS.r(cc))
               + f2  *cmplx(iUSQ.i(cc), -iUSQ.r(cc));
          iGamma.set(cc, -imag(qt), real(qt));
        }
      } else {
        // vout.general(m_vl,"force_recursive: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
        iGamma.zero();
        e_iQ.unit();
      }

      iGamma.at();
      iTheta_tmp.mult_nn(iGamma, U_tmp);
      iTheta.set_mat(site, mu, iTheta_tmp);

      Xi_tmp.mult_nn(Sigmap_tmp, e_iQ);
      Xi_tmp.multadd_dn(C_tmp, iGamma);
      Xi.set_mat(site, mu, Xi_tmp);
    }
  }

  /*
  unsigned long count;
  double time;
  KEK_FopCountFinish(id,&count,&time);
  m_time += time;
  m_flop += count;
  */
  double time1 = Communicator::get_time();
  m_time += time1 - time0;
}


//====================================================================
void Projection_Stout_SU3::set_fj(dcomplex& f0, dcomplex& f1, dcomplex& f2,
                                  const double& u, const double& w)
{
  dcomplex h0, h1, h2, e2iu, emiu, ixi0, qt;
  double   xi0, u2, w2, cos_w, fden;
  // double   c0, c1, c0max;

  xi0   = func_xi0(w);
  u2    = u * u;
  w2    = w * w;
  cos_w = cos(w);

  double cos_u = cos(u);
  double sin_u = sin(u);
  emiu = cmplx(cos_u, -sin_u);
  e2iu = cmplx(cos_u * cos_u - sin_u * sin_u, 2.0 * sin_u * cos_u);

  h0 = e2iu * cmplx(u2 - w2, 0.0)
       + emiu *cmplx(8.0 * u2 * cos_w, 2.0 * u * (3.0 * u2 + w2) * xi0);
  h1 = cmplx(2 * u, 0.0) * e2iu
       - emiu      *cmplx(2.0 * u * cos_w, -(3.0 * u2 - w2) * xi0);
  h2 = e2iu - emiu *cmplx(cos_w, 3.0 * u * xi0);

  fden = 1.0 / (9.0 * u2 - w2);
  f0   = h0 * fden;
  f1   = h1 * fden;
  f2   = h2 * fden;
}


//====================================================================
void Projection_Stout_SU3::set_uw(double& u, double& w,
                                  const Mat_SU_N& iQ2, const Mat_SU_N& iQ3)
{
  double c0, c1, c0max, theta;

  c0 = -(iQ3.i(0, 0) + iQ3.i(1, 1) + iQ3.i(2, 2)) / 3.0;
  c1 = -0.5 * (iQ2.r(0, 0) + iQ2.r(1, 1) + iQ2.r(2, 2));
  double c13r = sqrt(c1 / 3.0);
  c0max = 2.0 * c13r * c13r * c13r;

  theta = acos(c0 / c0max);
  u     = c13r * cos(theta / 3.0);
  w     = sqrt(c1) * sin(theta / 3.0);
}


//====================================================================
double Projection_Stout_SU3::func_xi0(double w)
{
  double epsilon_criterion = CommonParameters::epsilon_criterion();

  double xi0 = sin(w) / w;

  if (w < epsilon_criterion) {
    vout.general(m_vl, "%s: too small w = %18.6e\n", class_name.c_str(), w);
  }

  return xi0;
}


//====================================================================
double Projection_Stout_SU3::func_xi1(double w)
{
  double epsilon_criterion = CommonParameters::epsilon_criterion();

  double xi1 = cos(w) / (w * w) - sin(w) / (w * w * w);

  if (w < epsilon_criterion) {
    vout.general(m_vl, "%s: too small w = %18.6e\n", class_name.c_str(), w);
  }

  return xi1;
}


//====================================================================
void Projection_Stout_SU3::exp_iQ_bf(Field_G& e_iQ, const Field_G& iQ)
{
  // brute force version of exponentiation: for check

  int Nprec = 32;

  int Nvol = iQ.nvol();
  int Nex  = iQ.nex();

  Mat_SU_N u0(NC), u1(NC), u2(NC);
  Mat_SU_N h1(NC);

  for (int ex = 0; ex < Nex; ++ex) {
    for (int site = 0; site < Nvol; ++site) {
      u0.unit();
      u1.unit();
      h1 = iQ.mat(site, ex);

      for (int iprec = 0; iprec < Nprec; ++iprec) {
        double exf = 1.0 / (Nprec - iprec);
        u2  = h1 * u1;
        u2 *= exf;
        u1  = u2;
        u1 += u0;
      }

      u1.reunit();
      e_iQ.set_mat(site, ex, u1);
    }
  }
}


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