solver_GMRES_m_Cmplx.cpp

Go to the documentation of this file.

#include "solver_GMRES_m_Cmplx.h"

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

#ifdef USE_FACTORY
namespace {
  Solver *create_object(Fopr *fopr)
  {
    return new Solver_GMRES_m_Cmplx(fopr);
  }


  bool init = Solver::Factory::Register("GMRES_m_Cmplx", create_object);
}
#endif

//- 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_int("number_of_orthonormal_vectors", 0);

    param.Register_string("verbose_level", "NULL");
  }


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

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

//====================================================================
void Solver_GMRES_m_Cmplx::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    N_M;

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

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


  set_parameters(Niter, Stop_cond);
  set_parameters_GMRES_m(N_M);
}


//====================================================================
void Solver_GMRES_m_Cmplx::set_parameters(const int Niter, const double Stop_cond)
{
  //- print input parameters
  vout.general(m_vl, "Parameters of Solver_GMRES_m_Cmplx:\n");
  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, "Solver_GMRES_m_Cmplx: parameter range check failed.\n");
    abort();
  }

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


//====================================================================
void Solver_GMRES_m_Cmplx::set_parameters_GMRES_m(const int N_M)
{
  //- print input parameters
  vout.general(m_vl, "  N_M   = %d\n", N_M);

  //- range check
  int err = 0;
  err += ParameterCheck::non_negative(N_M);

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

  //- store values
  m_N_M = N_M;
}


//====================================================================
void Solver_GMRES_m_Cmplx::solve(Field& xq, const Field& b,
                                 int& Nconv, double& diff)
{
  vout.detailed(m_vl, "  GMRES_m_Cmplx solver starts\n");

  reset_field(b);

  vout.paranoiac(m_vl, "    norm of b = %16.8e\n", b.norm2());
  vout.paranoiac(m_vl, "    size of b = %d\n", b.size());

  double snorm = 1.0 / b.norm2();
  double rr;

  Nconv = -1;
  s     = b;

  solve_init(b, rr);

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

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

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

    if (rr * snorm < m_Stop_cond) {
      s    = m_fopr->mult(x);
      s   -= b;
      diff = s.norm();

      if (diff * diff * snorm < m_Stop_cond) {
        Nconv = m_N_M * (iter + 1);
        break;
      }

      s = x;
      solve_init(b, rr);
    }
  }
  if (Nconv == -1) {
    vout.crucial(m_vl, "GMRES_m_Cmplx solver not converged.\n");
    abort();
  }

  s    = m_fopr->mult(x);
  s   -= b;
  diff = s.norm();

  xq = x;
}


//====================================================================
void Solver_GMRES_m_Cmplx::reset_field(const Field& b)
{
  int Nin  = b.nin();
  int Nvol = b.nvol();
  int Nex  = b.nex();

  if ((s.nin() != Nin) || (s.nvol() != Nvol) || (s.nex() != Nex)) {
    s.reset(Nin, Nvol, Nex);
    r.reset(Nin, Nvol, Nex);
    x.reset(Nin, Nvol, Nex);

    v_tmp.reset(Nin, Nvol, Nex);

    v.resize(m_N_M + 1);

    for (int i = 0; i < m_N_M + 1; ++i) {
      v[i].reset(Nin, Nvol, Nex);
    }

    vout.paranoiac(m_vl, "  Solver_GMRES_m_Cmplx: field size reset.\n");
  }
}


//====================================================================
void Solver_GMRES_m_Cmplx::solve_init(const Field& b, double& rr)
{
  x = s;

  v_tmp = m_fopr->mult(x);
  r     = b;
  r    -= v_tmp;

  rr     = r * r;
  beta_p = sqrt(rr);

  v[0] = (1.0 / beta_p) * r;
}


//====================================================================
void Solver_GMRES_m_Cmplx::solve_step(const Field& b, double& rr)
{
  int    ij;
  double const_r, const_i;
  std::valarray<dcomplex> h((m_N_M + 1) * m_N_M), y(m_N_M);

  h = cmplx(0.0, 0.0);
  y = cmplx(0.0, 0.0);


  for (int j = 0; j < m_N_M; ++j) {
    v_tmp = m_fopr->mult(v[j]);

    for (int i = 0; i < (j + 1); ++i) {
      // h[i,j] = (A v[j], v[i]);
      ij = index_ij(i, j);
      innerprod_c(const_r, const_i, v_tmp, v[i]);
      h[ij] = cmplx(const_r, const_i);
    }

    // v[j+1] = A v[j] - \Sum_{i=0}^{j-1} h[i,j] * v[i]
    v[j + 1] = v_tmp;

    for (int i = 0; i < (j + 1); ++i) {
      // v[j+1] -= h[i,j] * v[i];
      ij = index_ij(i, j);
      mult_c(v_tmp, v[i], real(h[ij]), imag(h[ij]));
      v[j + 1] -= v_tmp;
    }

    const_r = v[j + 1] * v[j + 1];

    ij        = index_ij(j + 1, j);
    h[ij]     = sqrt(const_r);
    v[j + 1] /= sqrt(const_r);
  }


  // Compute y, which minimizes J := |r_new| = |beta_p - h * y|
  min_J(y, h);


  // x += Sum_{i=0}^{N_M-1} y[i] * v[i];
  for (int i = 0; i < m_N_M; ++i) {
    mult_c(v_tmp, v[i], real(y[i]), imag(y[i]));
    x += v_tmp;
  }


  // r = b - m_fopr->mult(x);
  s = x;
  solve_init(b, rr);
}


//====================================================================
void Solver_GMRES_m_Cmplx::innerprod_c(double& prod_r, double& prod_i,
                                       const Field& v, const Field& w)
{
  // prod = (v,w);

  int size = w.size();

  assert(v.size() == size);

  prod_r = 0.0;
  prod_i = 0.0;

  for (int i = 0; i < size; i += 2) {
    prod_r += v.cmp(i) * w.cmp(i) + v.cmp(i + 1) * w.cmp(i + 1);
    prod_i += v.cmp(i) * w.cmp(i + 1) - v.cmp(i + 1) * w.cmp(i);
  }

  prod_r = Communicator::reduce_sum(prod_r);
  prod_i = Communicator::reduce_sum(prod_i);
}


//====================================================================
void Solver_GMRES_m_Cmplx::min_J(std::valarray<dcomplex>& y,
                                 std::valarray<dcomplex>& h)
{
  // Compute y, which minimizes J := |r_new| = |beta_p - h * y|

  int      ii, i1i, ij, i1j;
  double   const_r, const_1_r, const_2_r;
  dcomplex cs, sn, const_1_c, const_2_c;

  std::valarray<dcomplex> g(m_N_M + 1);

  g    = dcomplex(0.0);
  g[0] = beta_p;


  for (int i = 0; i < m_N_M; ++i) {
    ii        = index_ij(i, i);
    const_1_r = abs(h[ii]);

    i1i       = index_ij(i + 1, i);
    const_2_r = abs(h[i1i]);

    const_r = sqrt(const_1_r * const_1_r
                   + const_2_r * const_2_r);

    cs = h[ii] / const_r;
    sn = h[i1i] / const_r;

    for (int j = i; j < m_N_M; ++j) {
      ij  = index_ij(i, j);
      i1j = index_ij(i + 1, j);

      const_1_c = conj(cs) * h[ij] + sn * h[i1j];
      const_2_c = -sn * h[ij] + cs * h[i1j];

      h[ij]  = const_1_c;
      h[i1j] = const_2_c;
    }

    const_1_c = conj(cs) * g[i] + sn * g[i + 1];
    const_2_c = -sn * g[i] + cs * g[i + 1];

    g[i]     = const_1_c;
    g[i + 1] = const_2_c;
  }


  for (int i = m_N_M - 1; i > -1; --i) {
    for (int j = i + 1; j < m_N_M; ++j) {
      ij    = index_ij(i, j);
      g[i] -= h[ij] * y[j];
    }

    ii   = index_ij(i, i);
    y[i] = g[i] / h[ii];
  }
}


//====================================================================
void Solver_GMRES_m_Cmplx::mult_c(Field& v,
                                  const Field& w,
                                  const double& prod_r, const double& prod_i)
{
  // v = dcomplex(prod_r,prod_i) * w;

  int size = w.size();

  assert(v.size() == size);

  double vr, vi;
  for (int i = 0; i < size; i += 2) {
    vr = prod_r * w.cmp(i) - prod_i * w.cmp(i + 1);
    vi = prod_r * w.cmp(i + 1) + prod_i * w.cmp(i);

    v.set(i, vr);
    v.set(i + 1, vi);
  }
}


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