solver_BiCGStab.cpp

Go to the documentation of this file.

#include "solver_BiCGStab.h"


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


  bool init = Solver::Factory::Register("BiCGStab", 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_string("verbose_level", "NULL");
  }


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

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

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

//====================================================================
void Solver_BiCGStab::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());
    abort();
  }


  set_parameters(Niter, Stop_cond);
}


//====================================================================
void Solver_BiCGStab::set_parameters(const int Niter, const double Stop_cond)
{
  ThreadManager_OpenMP::assert_single_thread(class_name);

  //- print input parameters
  vout.general(m_vl, "%s: input parameters\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());
    abort();
  }

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


//====================================================================
void Solver_BiCGStab::solve(Field& xq, const Field& b,
                            int& Nconv, double& diff)
{
  //#pragma omp parallel
  {
    double bnorm2 = b.norm2();
    double snorm  = 1.0 / bnorm2;
    int    bsize  = b.size();

    int nth = ThreadManager_OpenMP::get_num_threads();
    int ith = ThreadManager_OpenMP::get_thread_id();

#pragma omp master
    {
      vout.paranoiac(m_vl, "%s: solver starts\n", class_name.c_str());
      vout.paranoiac(m_vl, "  norm of b = %16.8e\n", bnorm2);
      vout.paranoiac(m_vl, "  size of b = %d\n", bsize);
    }
#pragma omp barrier


    reset_field(b);

    copy(s, b); //  s = b;

    double rr;
    int    Nconv2 = -1;

    solve_init(b, rr);

    bool is_converged = false;

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


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

#pragma omp master
        vout.detailed(m_vl, "    iter: %8d  %22.15e\n", 2 * (iter + 1), rr * snorm);
#pragma omp barrier

        if (rr * snorm < m_Stop_cond) {
          m_fopr->mult(s, x); // s  = m_fopr->mult(x);
          axpy(s, -1.0, b);   // s -= b;

          double diff2 = s.norm2();
          ThreadManager_OpenMP::sync_barrier_all();

          nth = ThreadManager_OpenMP::get_num_threads();
          if (ith == 0) vout.detailed(m_vl, "    iter0: %8d  %22.15e\n", nth, diff2 * snorm);

          if (diff2 * snorm < m_Stop_cond) {
            Nconv2 = 2 * (iter + 1);

            // break;
            is_converged = true;
          }

          copy(s, x); //  s = x;
          solve_init(b, rr);

          nth = ThreadManager_OpenMP::get_num_threads();
          if (ith == 0) vout.detailed(m_vl, "    iter1: %8d  %22.15e\n", nth, rr * snorm);
        }
      }
    }
    if (Nconv2 == -1) {
#pragma omp master
      vout.crucial(m_vl, "%s: not converged.\n", class_name.c_str());
#pragma omp barrier
      abort();
    }

    m_fopr->mult(p, x); // p  = m_fopr->mult(x);
    axpy(p, -1.0, b);   // p -= b;

    copy(xq, x);        // xq = x;

    double diff2 = p.norm2();

#pragma omp master
    {
      diff  = sqrt(diff2);
      Nconv = Nconv2;
    }
#pragma omp barrier
  } // end of parallel region
}


//====================================================================
void Solver_BiCGStab::reset_field(const Field& b)
{
#pragma omp master
  {
    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);
      p.reset(Nin, Nvol, Nex);
      v.reset(Nin, Nvol, Nex);
      t.reset(Nin, Nvol, Nex);
      rh.reset(Nin, Nvol, Nex);

      vout.detailed(m_vl, "    %s: field size reset.\n", class_name.c_str());
    }
  }
#pragma omp barrier
}


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

  //- r = b - A x_0
  m_fopr->mult(v, s); // v = m_fopr->mult(s);
  copy(r, b);         // r  = b;
  axpy(r, -1.0, v);   // r -= v;
  copy(rh, r);        // rh = r;

  rr = r.norm2();     // rr = r * r;

  p.set(0.0);         // p = 0.0
  v.set(0.0);         // v = 0.0

#pragma omp master
  {
    rho_prev   = 1.0;
    alpha_prev = 1.0;
    omega_prev = 1.0;
  }
#pragma omp barrier
}


//====================================================================
void Solver_BiCGStab::solve_step(double& rr)
{
  double rho = dot(rh, r);  // double rho = rh * r;
  double bet = rho * alpha_prev / (rho_prev * omega_prev);

  // p = r + bet * (p - omega_prev * v);
  axpy(p, -omega_prev, v);   // p += - omega_prev * v;
  aypx(bet, p, r);           // p  = bet * p + r;

  m_fopr->mult(v, p);        // v  = m_fopr->mult(p);

  double aden  = dot(rh, v); // dcomplex aden = rh * v;
  double alpha = rho / aden;

  copy(s, r);                 // s  = r
  axpy(s, -alpha, v);         // s += - alpha * v;

  m_fopr->mult(t, s);         // t  = m_fopr->mult(s);

  double omega_d = dot(t, t); // omega_d = t * t;
  double omega_n = dot(t, s); // omega_n = t * s;
  double omega   = omega_n / omega_d;

  axpy(x, omega, s);   // x += omega * s;
  axpy(x, alpha, p);   // x += alpha * p;

  copy(r, s);          // r  = s
  axpy(r, -omega, t);  // r += - omega * t;

  rr = r.norm2();      // rr = r * r;

#pragma omp master
  {
    rho_prev   = rho;
    alpha_prev = alpha;
    omega_prev = omega;
  }
#pragma omp barrier
}


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