solver_BiCGStab_DS_L_Cmplx.cpp

Go to the documentation of this file.

#include "solver_BiCGStab_DS_L_Cmplx.h"

using std::valarray;


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


  bool init = Solver::Factory::Register("BiCGStab_DS_L_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_double("tolerance_for_DynamicSelection_of_L", 0);

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


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

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

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

//====================================================================
void Solver_BiCGStab_DS_L_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_L;
  double Tol_L;

  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_L);
  err += params.fetch_double("tolerance_for_DynamicSelection_of_L", Tol_L);

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


  set_parameters(Niter, Stop_cond);
  set_parameters_DS_L(N_L, Tol_L);
}


//====================================================================
void Solver_BiCGStab_DS_L_Cmplx::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_DS_L_Cmplx::set_parameters_DS_L(const int N_L, const double Tol_L)
{
  //- print input parameters
  vout.general(m_vl, "  N_L   = %d\n", N_L);
  vout.general(m_vl, "  Tol_L = %16.8e\n", Tol_L);

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

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

  //- store values
  m_N_L   = N_L;
  m_Tol_L = Tol_L;
}


//====================================================================
void Solver_BiCGStab_DS_L_Cmplx::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();
    double rr;

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

#pragma omp master
    {
      vout.paranoiac(m_vl, "%s: 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);


    // Nconv = -1;
    int Nconv2 = 0;
    copy(s, b); // s = b;

    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);

        Nconv2 += 2 * N_L_prev;

#pragma omp master
        {
          vout.detailed(m_vl, "    iter,N_L: %8d  %8d\n", Nconv2, N_L_prev);
          vout.detailed(m_vl, "    iter: %8d  %22.15e\n", Nconv2, 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) {
            // NB. Nconv is calculated above.
            // 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);
        }
      }
    }

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

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

    double diff2 = s.norm2();

#pragma omp master
    {
      diff  = sqrt(diff2);
      Nconv = Nconv2;
    }
#pragma omp barrier

    if (diff2 * snorm > m_Stop_cond) {
#pragma omp master
      vout.crucial(m_vl, "%s: not converged.\n", class_name.c_str());
#pragma omp barrier
      abort();
    }
  } // end of parallel region
}


//====================================================================
void Solver_BiCGStab_DS_L_Cmplx::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);
      x.reset(Nin, Nvol, Nex);

      r_init.reset(Nin, Nvol, Nex);

      v_tmp.reset(Nin, Nvol, Nex);

      vout.paranoiac(m_vl, "%s: field size reset.\n", class_name.c_str());
    }

    u.resize(m_N_L + 1);
    r.resize(m_N_L + 1);

    for (int i = 0; i < m_N_L + 1; ++i) {
      u[i].reset(Nin, Nvol, Nex);
      r[i].reset(Nin, Nvol, Nex);
    }
  }
#pragma omp barrier
}


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

  //- r[0] = b - A x_0
  m_fopr->mult(v_tmp, s);   // v_tmp = m_fopr->mult(s);
  copy(r[0], b);            // r[0]  = b;
  axpy(r[0], -1.0, v_tmp);  // r[0] -= v_tmp;

  copy(r_init, r[0]);       // r_init = r[0];
  rr = r[0].norm2();        // rr     = r[0] * r[0];

  u[0].set(0.0);            // u[0] = 0.0;

  // NB. alpha_prev = 0.0 \neq 1.0
#pragma omp master
  {
    rho_prev   = cmplx(-1.0, 0.0);
    alpha_prev = cmplx(0.0, 0.0);

    N_L_prev = m_N_L;
  }
#pragma omp barrier
}


//====================================================================
void Solver_BiCGStab_DS_L_Cmplx::solve_step(double& rr)
{
  dcomplex rho_prev2   = rho_prev;
  dcomplex alpha_prev2 = alpha_prev;

  int      N_L_tmp         = 0; // superficial initialization
  dcomplex c_Rayleigh_prev = cmplx(0.0, 0.0);

  bool is_converged_L = false;


  for (int j = 0; j < N_L_prev; ++j) {
    if (!is_converged_L) {
      dcomplex rho = dotc(r[j], r_init);   // rho  = r[j] * r_init;
      rho = conj(rho);

      dcomplex beta = alpha_prev2 * (rho / rho_prev2);

      rho_prev2 = rho;

      for (int i = 0; i < j + 1; ++i) {
        aypx(-beta, u[i], r[i]);   // u[i] = - beta * u[i] + r[i];
      }

      m_fopr->mult(u[j + 1], u[j]);  // u[j+1] = m_fopr->mult(u[j]);

      dcomplex gamma = dotc(u[j + 1], r_init);
      alpha_prev2 = rho_prev2 / conj(gamma);

      for (int i = 0; i < j + 1; ++i) {
        axpy(r[i], -alpha_prev2, u[i + 1]);  // r[i] -= alpha_prev * u[i+1];
      }

      m_fopr->mult(r[j + 1], r[j]);        // r[j+1] = m_fopr->mult(r[j]);

      axpy(x, alpha_prev2, u[0]);          // x += alpha_prev * u[0];


      //- calculate the Rayleigh quotient for N_L_tmp

      // dcomplex c_Rayleigh = (r[j] * r[j+1]) / (r[j] * r[j]);
      double   r_tmp      = r[j].norm2();                 // r_tmp = r[j] * r[j];
      dcomplex c_Rayleigh = dotc(r[j], r[j + 1]) / r_tmp; // c_Rayleigh = r[j] * r[j+1] / r_tmp;

      dcomplex c_E = (c_Rayleigh - c_Rayleigh_prev) / c_Rayleigh;

      c_Rayleigh_prev = c_Rayleigh;

      N_L_tmp = j + 1;

      // vout.paranoiac(m_vl, "N_L_tmp,abs(c_E),m_Tol_L = %d %f %f\n",N_L_tmp,abs(c_E),m_Tol_L);

      if (abs(c_E) < m_Tol_L) {
        // vout.paranoiac(m_vl, "N_L_tmp = %d\n",N_L_tmp);

        // break;
        is_converged_L = true;
      }
    }
  }


  valarray<double>   sigma(m_N_L + 1);
  valarray<dcomplex> gamma_prime(m_N_L + 1);

  // NB. tau(m_N_L,m_N_L+1), not (m_N_L+1,m_N_L+1)
  valarray<dcomplex> tau(m_N_L * (m_N_L + 1));
  int                ij, ji;

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

      dcomplex r_ji = dotc(r[j], r[i]);
      tau[ij] = conj(r_ji) / sigma[i];  // tau[ij]  = (r[j] * r[i]) / sigma[i];
      axpy(r[j], -tau[ij], r[i]);       // r[j]    -= tau[ij] * r[i];
    }

    sigma[j] = r[j].norm2();  // sigma[j] = r[j] * r[j];

    dcomplex r_0j = dotc(r[0], r[j]);
    gamma_prime[j] = conj(r_0j) / sigma[j]; // gamma_prime[j] = (r[0] * r[j]) / sigma[j];
  }


  valarray<dcomplex> gamma(m_N_L + 1);
  dcomplex           c_tmp;

  gamma[N_L_tmp] = gamma_prime[N_L_tmp];

  for (int j = N_L_tmp - 1; j > 0; --j) {
    c_tmp = cmplx(0.0, 0.0);

    for (int i = j + 1; i < N_L_tmp + 1; ++i) {
      ji     = index_ij(j, i);
      c_tmp += tau[ji] * gamma[i];
    }

    gamma[j] = gamma_prime[j] - c_tmp;
  }


  // NB. gamma_double_prime(m_N_L), not (m_N_L+1)
  valarray<dcomplex> gamma_double_prime(m_N_L);

  for (int j = 1; j < N_L_tmp; ++j) {
    c_tmp = cmplx(0.0, 0.0);

    for (int i = j + 1; i < N_L_tmp; ++i) {
      ji     = index_ij(j, i);
      c_tmp += tau[ji] * gamma[i + 1];
    }

    gamma_double_prime[j] = gamma[j + 1] + c_tmp;
  }


  axpy(x, gamma[1], r[0]);                        // x    += gamma[            1] * r[      0];
  axpy(r[0], -gamma_prime[N_L_tmp], r[N_L_tmp]);  // r[0] -= gamma_prime[N_L_tmp] * r[N_L_tmp];
  axpy(u[0], -gamma[N_L_tmp], u[N_L_tmp]);        // u[0] -= gamma[      N_L_tmp] * u[N_L_tmp];

  for (int j = 1; j < N_L_tmp; ++j) {
    axpy(x, gamma_double_prime[j], r[j]);      // x    += gamma_double_prime[j] * r[j];
    axpy(r[0], -gamma_prime[j], r[j]);         // r[0] -= gamma_prime[       j] * r[j];
    axpy(u[0], -gamma[j], u[j]);               // u[0] -= gamma[             j] * u[j];
  }

  rr = r[0].norm2();  // rr = r[0] * r[0];

  // #pragma omp master
  int ith = ThreadManager_OpenMP::get_thread_id();
  if (ith == 0)
  {
    rho_prev   = rho_prev2;
    alpha_prev = alpha_prev2;
    N_L_prev   = N_L_tmp;
    rho_prev  *= -gamma_prime[N_L_tmp];
  }
#pragma omp barrier
}


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