fopr_Overlap.cpp

Go to the documentation of this file.

#include "fopr_Overlap.h"

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

using std::valarray;

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_double("quark_mass", 0.0);
    param.Register_double("domain_wall_height", 0.0);
    param.Register_int("number_of_poles", 0);
    param.Register_double("lower_bound", 0.0);
    param.Register_double("upper_bound", 0.0);
    param.Register_int("maximum_number_of_iteration", 0);
    param.Register_double("convergence_criterion_squared", 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.Overlap", append_entry);
#endif
}
//- end

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

//====================================================================
void Fopr_Overlap::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        mq, M0;
  int           Np;
  double        x_min, x_max;
  int           Niter;
  double        Stop_cond;
  valarray<int> bc;

  int err = 0;
  err += params.fetch_double("quark_mass", mq);
  err += params.fetch_double("domain_wall_height", M0);
  err += params.fetch_int("number_of_poles", Np);
  err += params.fetch_double("lower_bound", x_min);
  err += params.fetch_double("upper_bound", x_max);
  err += params.fetch_int("maximum_number_of_iteration", Niter);
  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
  err += params.fetch_int_vector("boundary_condition", bc);

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


  set_parameters(mq, M0, Np, x_min, x_max, Niter, Stop_cond, bc);
}


//====================================================================
void Fopr_Overlap::set_parameters(const double mq, const double M0, const int Np,
                                  const double x_min, const double x_max,
                                  const int Niter, const double Stop_cond,
                                  const std::valarray<int> bc)
{
  int Ndim = CommonParameters::Ndim();

  //- print input parameters
  vout.general(m_vl, "Fopr_Overlap paramters:\n");
  vout.general(m_vl, "  mq        = %10.6f\n", mq);
  vout.general(m_vl, "  M0        = %10.6f\n", M0);
  vout.general(m_vl, "  Np        = %4d\n", Np);
  vout.general(m_vl, "  x_min     = %12.8f\n", x_min);
  vout.general(m_vl, "  x_max     = %12.6f\n", x_max);
  vout.general(m_vl, "  Niter     = %6d\n", Niter);
  vout.general(m_vl, "  Stop_cond = %12.4e\n", Stop_cond);
  for (int mu = 0; mu < Ndim; ++mu) {
    vout.general(m_vl, "  boundary[%d] = %2d\n", mu, bc[mu]);
  }

  //- range check
  int err = 0;
  err += ParameterCheck::non_zero(mq);
  err += ParameterCheck::non_zero(M0);
  err += ParameterCheck::non_zero(Np);
  // NB. x_min,x_max == 0 is allowed.
  err += ParameterCheck::non_zero(Niter);
  err += ParameterCheck::square_non_zero(Stop_cond);

  if (err) {
    vout.crucial(m_vl, "Fopr_Overlap: parameter range check failed.\n");
    abort();
  }

  assert(bc.size() == Ndim);

  //- store values
  m_mq = mq;
  m_M0 = M0;

  m_Np    = Np;
  m_x_min = x_min;
  m_x_max = x_max;

  m_Niter     = Niter;
  m_Stop_cond = Stop_cond;

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

  //- propagate parameters
  m_sign->set_parameters(m_Np, m_x_min, m_x_max, m_Niter, m_Stop_cond);
}


//====================================================================
void Fopr_Overlap::set_lowmodes(int Nsbt, valarray<double> *ev, valarray<Field> *vk)
{
  m_Nsbt = Nsbt;
  m_ev   = ev;
  m_vk   = vk;

  m_sign->set_lowmodes(m_Nsbt, m_ev, m_vk);
}


//====================================================================
void Fopr_Overlap::DdagD(Field& v, const Field& b)
{
  Field w(b.nin(), b.nvol(), b.nex());

  H(w, b);
  H(v, w);
}


//====================================================================
void Fopr_Overlap::H(Field& v, const Field& b)
{
  Field w(b.nin(), b.nvol(), b.nex());

  D(w, b);
  m_fopr_w->mult_gm5(v, w);
}


//====================================================================
void Fopr_Overlap::Ddag(Field& v, const Field& b)
{
  Field w(b.nin(), b.nvol(), b.nex());

  m_fopr_w->mult_gm5(w, b);
  H(v, w);
}


//====================================================================
void Fopr_Overlap::D(Field& v, const Field& b)
{
  assert(b.nin() == m_fopr_w->field_nin());
  assert(b.nvol() == m_fopr_w->field_nvol());
  assert(b.nex() == m_fopr_w->field_nex());

  assert(v.nin() == m_fopr_w->field_nin());
  assert(v.nvol() == m_fopr_w->field_nvol());
  assert(v.nex() == m_fopr_w->field_nex());

  double f1 = m_M0 + 0.5 * m_mq;
  double f2 = m_M0 - 0.5 * m_mq;

  Field w(b.nin(), b.nvol(), b.nex());

  m_sign->mult(w, b);
  m_fopr_w->mult_gm5(v, w);

  v *= f2;
  v += f1 * b;
}


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