fopr_Clover_eo.cpp

Go to the documentation of this file.

#include "fopr_Clover_eo.h"

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

using std::valarray;

//====================================================================
//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("hopping_parameter", 0.0);
    param.Register_double("clover_coefficient", 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.Clover_eo", append_entry);
#endif
}
//- end

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

//====================================================================
void Fopr_Clover_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, cSW;
  valarray<int> bc;

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

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


  set_parameters(kappa, cSW, bc);
}


//====================================================================
void Fopr_Clover_eo::set_parameters(const double kappa, const double cSW,
                                    const std::valarray<int> bc)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of Fopr_Clover_eo:\n");
  vout.general(m_vl, "  kappa = %8.4f\n", kappa);
  vout.general(m_vl, "  cSW   = %8.4f\n", cSW);
  for (int mu = 0; mu < m_Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

  //- range check
  // NB. kappa,cSW == 0 is allowed.
  assert(bc.size() == m_Ndim);

  //- store values
  m_kappa = kappa;
  m_cSW   = cSW;

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

  //- propagate parameters to components
  m_fopr_w->set_parameters(m_kappa, m_boundary);
  m_fopr_csw->set_parameters(m_kappa, m_cSW, m_boundary);
}


//====================================================================
void Fopr_Clover_eo::set_config(Field *U)
{
  int Ndim = CommonParameters::Ndim();
  int Nvol = CommonParameters::Nvol();

  // if (m_Ueo == 0) {
  //   Field_G* Ueo = new Field_G(Nvol,Ndim);
  //   m_idx.convertField(*Ueo, *U);
  //   m_Ueo = Ueo;
  // }
  m_idx.convertField(*m_Ueo, *U);

  // m_fopr_w->set_config(  m_Ueo);
  m_fopr_w->set_config(U);
  m_fopr_csw->set_config(m_Ueo);

  solve_csw_inv();
}


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

  m_idx.convertField(tmp, b, 0);
  be = (Field)mult_csw_inv((Field_F)tmp, 0);

  m_idx.convertField(tmp, b, 1);
  bo = (Field)mult_csw_inv((Field_F)tmp, 1);

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


void Fopr_Clover_eo::postPropD(Field& x,
                               const Field& xe, const Field& bo)
{
  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_idx.reverseField(x, xe, 0);
  m_idx.reverseField(x, xo, 1);
}


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

  m_idx.convertField(tmp, b, 0);
  be = (Field)mult_csw_inv((Field_F)tmp, 0);
  m_idx.convertField(tmp, b, 1);
  bo = (Field)mult_csw_inv((Field_F)tmp, 1);

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


void Fopr_Clover_eo::postPropDag(Field& x,
                                 const Field& xe, const Field& bo)
{
  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_idx.reverseField(x, xe, 0);
  m_idx.reverseField(x, xo, 1);
}


//====================================================================
void Fopr_Clover_eo::solve_csw_inv()
{
  double eps2 = CommonParameters::epsilon_criterion2();

  Parameters *params_solver = ParametersFactory::New("Solver");

  params_solver->set_string("solver_type", "CG");
  params_solver->set_int("maximum_number_of_iteration", 1000);
  params_solver->set_double("convergence_criterion_squared", 1.0e-30);
  //- NB. set VerboseLevel to CRUCIAL to suppress frequent messages.
  params_solver->set_string("verbose_level", "Crucial");

  int    Nconv;
  double diff;

  Solver *solver = new Solver_CG(m_fopr_csw);
  solver->set_parameters(*params_solver);

  Field_F w(m_Nvol2);
  Field_F w2(m_Nvol2);

  for (int ispin = 0; ispin < m_Nd; ++ispin) {
    for (int icolor = 0; icolor < m_Nc; ++icolor) {
      int spin_color = icolor + m_Nc * ispin;
      w = 0.0;
      for (int isite = 0; isite < m_Nvol2; ++isite) {
        w.set_ri(icolor, ispin, isite, 0, 1, 0);
      }

      if (m_cSW * m_cSW < eps2) {
        m_fee_inv->setpart_ex(spin_color, w, 0);
        m_foo_inv->setpart_ex(spin_color, w, 0);
      } else {
        m_fopr_csw->set_mode("even");
        solver->solve(w2, w, Nconv, diff);
        m_fee_inv->setpart_ex(spin_color, w2, 0);

        vout.detailed(m_vl, "  Nconv,diff = %d %12.4e\n", Nconv, diff);

        m_fopr_csw->set_mode("odd");
        solver->solve(w2, w, Nconv, diff);
        m_foo_inv->setpart_ex(spin_color, w2, 0);

        vout.detailed(m_vl, "  Nconv,diff = %d %12.4e\n", Nconv, diff);
      }
    }
  }

  delete params_solver;
  delete solver;
}


//====================================================================
const Field_F Fopr_Clover_eo::mult_csw_inv(const Field_F& f, const int ieo)
{
  int     nex = f.nex();
  Field_F w(m_Nvol2, nex);
  Field_F *csw_inv;

  if (ieo == 0) {
    csw_inv = m_fee_inv;
  } else if (ieo == 1) {
    csw_inv = m_foo_inv;
  } else {
    vout.crucial(m_vl, "Fopr_clover_eo: wrong parameter, ieo=%d.\n", ieo);
    abort();
  }

  for (int iex = 0; iex < nex; ++iex) {
    for (int isite = 0; isite < m_Nvol2; ++isite) {
      for (int ispin = 0; ispin < m_Nd; ++ispin) {
        for (int icolor = 0; icolor < m_Nc; ++icolor) {
          double re = 0;
          double im = 0;

          for (int jspin = 0; jspin < m_Nd; ++jspin) {
            for (int jcolor = 0; jcolor < m_Nc; ++jcolor) {
              int spin_color = jcolor + m_Nc * jspin;

              re += csw_inv->cmp_r(icolor, ispin, isite, spin_color) *
                    f.cmp_r(jcolor, jspin, isite, iex);
              re -= csw_inv->cmp_i(icolor, ispin, isite, spin_color) *
                    f.cmp_i(jcolor, jspin, isite, iex);

              im += csw_inv->cmp_r(icolor, ispin, isite, spin_color) *
                    f.cmp_i(jcolor, jspin, isite, iex);
              im += csw_inv->cmp_i(icolor, ispin, isite, spin_color) *
                    f.cmp_r(jcolor, jspin, isite, iex);
            }
          }

          w.set_ri(icolor, ispin, isite, iex, re, im);
        }
      }
    }
  }

  return w;
}


//====================================================================
const Field_G Fopr_Clover_eo::trSigmaInv(const int mu, const int nu)
{
  int      nex_finv = m_fee_inv->nex();
  Vec_SU_N v;
  Field_F  sigma_inv(m_Nvol, nex_finv);
  Field_G  tr_sigma_inv(m_Nvol, 1);
  {
    Field_F sigma_eo_inv(m_Nvol2, nex_finv);
    mult_isigma(sigma_eo_inv, *m_fee_inv, mu, nu);
    m_idx.reverseField(sigma_inv, sigma_eo_inv, 0);
    mult_isigma(sigma_eo_inv, *m_foo_inv, mu, nu);
    m_idx.reverseField(sigma_inv, sigma_eo_inv, 1);
  }

  for (int isite = 0; isite < m_Nvol; ++isite) {
    for (int ispin = 0; ispin < m_Nd; ++ispin) {
      for (int icolor = 0; icolor < m_Nc; ++icolor) {
        v = sigma_inv.vec(ispin, isite, icolor + m_Nc * ispin);
        for (int jcolor = 0; jcolor < m_Nc; ++jcolor) {
          int cc = icolor + m_Nc * jcolor;
          tr_sigma_inv.set_r(cc, isite, 0, v.r(jcolor));
          tr_sigma_inv.set_i(cc, isite, 0, v.i(jcolor));
        }
      }
    }
  }
  return tr_sigma_inv;
}


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