shiftsolver_CG.cpp

Go to the documentation of this file.

#include "shiftsolver_CG.h"

//- parameter entries
namespace {
  void append_entry(Parameters& param)
  {
    param.Register_int("maximum_number_of_iteration", 0);
    param.Register_double("convergence_criterion_squared", 0.0);

    param.Register_string("verbose_level", "NULL");
  }
}
//- end

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

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

//====================================================================
void Shiftsolver_CG::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
  int    Niter;
  double Stop_cond;

  int err = 0;
  err += params.fetch_int("maximum_number_of_iteration", Niter);
  err += params.fetch_double("convergence_criterion_squared", Stop_cond);

  if (err) {
    vout.crucial(m_vl, "%s: fetch error, input parameter not found.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }


  set_parameters(Niter, Stop_cond);
}


//====================================================================
void Shiftsolver_CG::set_parameters(const int Niter, const double Stop_cond)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
  vout.general(m_vl, "  Niter     = %d\n", Niter);
  vout.general(m_vl, "  Stop_cond = %16.8e\n", Stop_cond);

  //- range check
  int err = 0;
  err += ParameterCheck::non_negative(Niter);
  err += ParameterCheck::square_non_zero(Stop_cond);

  if (err) {
    vout.crucial(m_vl, "%s: parameter range check failed.\n", class_name.c_str());
    exit(EXIT_FAILURE);
  }

  //- store values
  m_Niter     = Niter;
  m_Stop_cond = Stop_cond;
}


//====================================================================
void Shiftsolver_CG::solve(std::vector<Field>& xq,
                           const std::vector<double>& sigma,
                           const Field& b,
                           int& Nconv, double& diff)
{
  int Nshift = sigma.size();

  vout.paranoiac(m_vl, "  Shift CG solver start.\n");
  vout.paranoiac(m_vl, "    number of shift = %d\n", Nshift);
  vout.paranoiac(m_vl, "    values of shift:\n");
  for (int i = 0; i < Nshift; ++i) {
    vout.paranoiac(m_vl, "    %d  %12.8f\n", i, sigma[i]);
  }

  snorm = 1.0 / b.norm2();

  int Nconv2 = -1;

  reset_field(b, sigma, Nshift);

  copy(s, b);
  copy(r, b);

  double rr = 0.0;

  solve_init(rr);

  vout.detailed(m_vl, "    iter: %8d  %22.15e\n", 0, rr * snorm);

  bool is_converged = false;

  for (int iter = 0; iter < m_Niter; iter++) {
    solve_step(rr);

    vout.detailed(m_vl, "    iter: %8d  %22.15e  %4d\n", (iter + 1), rr * snorm, Nshift2);

    if (rr * snorm < m_Stop_cond) {
      Nconv2       = iter;
      is_converged = true;
      break;
    }
  }

  if (!is_converged) {
    vout.crucial(m_vl, "Shiftsolver_CG not converged.\n");
    exit(EXIT_FAILURE);
  }


  std::vector<double> diffs(Nshift);
  for (int i = 0; i < Nshift; ++i) {
    diffs[i] = 0.0;
  }

  for (int i = 0; i < Nshift; ++i) {
    m_fopr->mult(s, x[i]);
    axpy(s, sigma[i], x[i]);
    axpy(s, -1.0, b);

    double diff1 = sqrt(s.norm2() * snorm);

    vout.paranoiac(m_vl, "    %4d  %22.15e\n", i, diff1);

    // if (diff1 > diff2) diff2 = diff1;
    diffs[i] = diff1;
  }

#pragma omp barrier
#pragma omp master
  {
    double diff2 = -1.0;

    for (int i = 0; i < Nshift; ++i) {
      if (diffs[i] > diff2) diff2 = diffs[i];
    }

    diff = diff2;

    Nconv = Nconv2;
  }
#pragma omp barrier

  for (int i = 0; i < Nshift; ++i) {
    copy(xq[i], x[i]);
  }

  vout.paranoiac(m_vl, "   diff(max) = %22.15e  \n", diff);
}


//====================================================================
void Shiftsolver_CG::solve_init(double& rr)
{
  int Nshift = p.size();

  vout.paranoiac(m_vl, "number of shift = %d\n", Nshift);

  for (int i = 0; i < Nshift; ++i) {
    copy(p[i], s);
    scal(x[i], 0.0);
  }

  copy(r, s);
  rr = r.norm2();

#pragma omp barrier
#pragma omp master
  {
    alpha_p = 0.0;
    beta_p  = 1.0;
  }
#pragma omp barrier
}


//====================================================================
void Shiftsolver_CG::solve_step(double& rr)
{
  m_fopr->mult(s, p[0]);
  axpy(s, m_sigma0, p[0]);

  double rr_p = rr;
  double pa_p = dot(s, p[0]);
  double beta = -rr_p / pa_p;

  axpy(x[0], -beta, p[0]);
  axpy(r, beta, s);
  rr = r.norm2();

  double alpha = rr / rr_p;

  aypx(alpha, p[0], r);

#pragma omp barrier
#pragma omp master
  {
    pp[0] = rr;
  }
#pragma omp barrier

  double alpha_h = 1.0 + alpha_p * beta / beta_p;

  for (int ish = 1; ish < Nshift2; ++ish) {
    double zeta = (alpha_h - csh2[ish] * beta) / zeta1[ish]
                  + (1.0 - alpha_h) / zeta2[ish];
    zeta = 1.0 / zeta;
    double zr      = zeta / zeta1[ish];
    double beta_s  = beta * zr;
    double alpha_s = alpha * zr * zr;

    axpy(x[ish], -beta_s, p[ish]);
    scal(p[ish], alpha_s);
    axpy(p[ish], zeta, r);

    double ppr = p[ish].norm2();

#pragma omp barrier
#pragma omp master
    {
      pp[ish] = ppr * snorm;

      zeta2[ish] = zeta1[ish];
      zeta1[ish] = zeta;
    }
#pragma omp barrier
  }

  int ish1 = Nshift2;

  for (int ish = Nshift2 - 1; ish >= 0; --ish) {
    vout.paranoiac(m_vl, "%4d %16.8e\n", ish, pp[ish]);
    if (pp[ish] > m_Stop_cond) {
      ish1 = ish + 1;
      break;
    }
  }

#pragma omp barrier
#pragma omp master
  {
    Nshift2 = ish1;

    alpha_p = alpha;
    beta_p  = beta;
  }
#pragma omp barrier
}


//====================================================================
void Shiftsolver_CG::reset_field(const Field& b, const std::vector<double>& sigma, const int Nshift)
{
#pragma omp barrier
#pragma omp master
  {
    int Nin  = b.nin();
    int Nvol = b.nvol();
    int Nex  = b.nex();

    p.resize(Nshift);
    x.resize(Nshift);
    zeta1.resize(Nshift);
    zeta2.resize(Nshift);
    csh2.resize(Nshift);
    pp.resize(Nshift);

    for (int i = 0; i < Nshift; ++i) {
      p[i].reset(Nin, Nvol, Nex);
      x[i].reset(Nin, Nvol, Nex);
      zeta1[i] = 1.0;
      zeta2[i] = 1.0;
      csh2[i]  = sigma[i] - sigma[0];

      pp[i] = 0.0;
    }

    s.reset(Nin, Nvol, Nex);
    r.reset(Nin, Nvol, Nex);

    m_sigma0 = sigma[0];

    Nshift2 = Nshift;
  }
#pragma omp barrier
}


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