Bridge++  Ver. 2.0.4
asolver_BiCGStab_Cmplx-tmpl.h
Go to the documentation of this file.
1 
10 template<typename AFIELD>
11 const std::string ASolver_BiCGStab_Cmplx<AFIELD>::class_name
12  = "ASolver_BiCGStab_Cmplx";
13 
14 namespace {
15 #ifndef AFIELD_HAS_SUB
16  template<typename AFIELD>
17  void sub(AFIELD& v, const AFIELD& w)
18  {
19  axpy(v, -1.0, w);
20  }
21 #endif
22 
23 #ifndef AFIELD_HAS_DOTC_AND_NORM2
24  template<typename AFIELD>
25  void dotc_and_norm2(typename AFIELD::complex_t& dotc_vw,
26  typename AFIELD::real_t& v_norm2,
27  typename AFIELD::real_t& w_norm2,
28  const AFIELD& v, const AFIELD& w)
29  {
30  dotc_vw = dotc(v, w);
31  v_norm2 = v.norm2();
32  w_norm2 = w.norm2();
33  }
34 #endif
35 }
36 
37 //====================================================================
38 template<typename AFIELD>
39 void ASolver_BiCGStab_Cmplx<AFIELD>::init(void)
40 {
41  ThreadManager::assert_single_thread(class_name);
42 
43  vout.detailed(m_vl, "%s: being setup.\n", class_name.c_str());
44 
45  int nin = m_fopr->field_nin();
46  int nvol = m_fopr->field_nvol();
47  int nex = m_fopr->field_nex();
48 
49  m_r.reset(nin, nvol, nex);
50  m_p.reset(nin, nvol, nex);
51  m_v.reset(nin, nvol, nex);
52  m_t.reset(nin, nvol, nex);
53  m_rh.reset(nin, nvol, nex);
54 
55  m_ecrit = 1.e-32;
56  m_nconv = -1;
57  m_Omega_tolerance = 0.7;
58 
59  vout.detailed(m_vl, "%s: setup finished.\n", class_name.c_str());
60 }
61 
62 
63 //====================================================================
64 template<typename AFIELD>
65 void ASolver_BiCGStab_Cmplx<AFIELD>::tidyup(void)
66 {
67  // ThreadManager::assert_single_thread(class_name);
68  // nothing is to be deleted.
69 }
70 
71 
72 //====================================================================
73 template<typename AFIELD>
74 void ASolver_BiCGStab_Cmplx<AFIELD>::set_parameters(const Parameters& params)
75 {
76  const string str_vlevel = params.get_string("verbose_level");
77 
78  m_vl = vout.set_verbose_level(str_vlevel);
79 
80  //- fetch and check input parameters
81  int Niter, Nrestart;
82  double Stop_cond;
83  double Omega_tolerance;
84  InitialGuess initial_guess;
85  int err = 0;
86  err += params.fetch_int("maximum_number_of_iteration", Niter);
87  err += params.fetch_int("maximum_number_of_restart", Nrestart);
88  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
89 
90  if (err) {
91  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
92  class_name.c_str());
93  exit(EXIT_FAILURE);
94  }
95  err = params.fetch_double("Omega_tolerance", Omega_tolerance);
96  if (err) {
97  Omega_tolerance = 0.7;
98  }
99  std::string initial_guess_str;
100  err = params.fetch_string("initial_guess_mode", initial_guess_str);
101  initial_guess = InitialGuess::RHS;
102  if (!err) {
103  vout.detailed(m_vl, "%s: initial_guess_str=%s\n",
104  class_name.c_str(), initial_guess_str.c_str());
105 
106  if (initial_guess_str == "ZERO") {
107  initial_guess = InitialGuess::ZERO;
108  } else if (initial_guess_str == "GIVEN") {
109  initial_guess = InitialGuess::GIVEN;
110  } else {
111  vout.crucial(m_vl, "%s: unknown initial guess mode string: %d\n",
112  class_name.c_str(), initial_guess_str.c_str());
113  exit(EXIT_FAILURE);
114  }
115  }
116 
117  // "restart" is not yet supported
118  int Niter2 = Niter * Nrestart;
119 
120  // N.B.
121  // - default initial mode is RHS
122  // - negative Stop_cond is for fixed nubmer of iterations
123  set_parameters(Niter2, Stop_cond, initial_guess);
124  set_parameters_BiCGStab_series(Omega_tolerance);
125 }
126 
127 
128 //====================================================================
129 template<typename AFIELD>
130 void ASolver_BiCGStab_Cmplx<AFIELD>::set_parameters(const int Niter,
131  const real_t Stop_cond)
132 {
133  set_parameters(Niter, Stop_cond, InitialGuess::RHS);
134 }
135 
136 
137 //====================================================================
138 template<typename AFIELD>
139 void ASolver_BiCGStab_Cmplx<AFIELD>::set_parameters(const int Niter,
140  const real_t Stop_cond,
141  const InitialGuess init_guess_mode)
142 {
143  ThreadManager::assert_single_thread(class_name);
144 
145  m_Niter = Niter;
146  m_Stop_cond = Stop_cond;
147  std::string prec = "double";
148  if (sizeof(real_t) == 4) prec = "float";
149  m_initial_mode = init_guess_mode;
150 
151  vout.general(m_vl, "%s:\n", class_name.c_str());
152  vout.general(m_vl, " Precision: %s\n", prec.c_str());
153  vout.general(m_vl, " Niter = %d\n", m_Niter);
154  vout.general(m_vl, " Stop_cond = %16.8e\n", m_Stop_cond);
155  vout.general(m_vl, " init_guess_mode: %d\n", m_initial_mode);
156 }
157 
158 
159 //====================================================================
160 template<typename AFIELD>
161 void ASolver_BiCGStab_Cmplx<AFIELD>::set_parameters_BiCGStab_series(const real_t Omega_tolerance)
162 {
163  ThreadManager::assert_single_thread(class_name);
164 
165  //- print input parameters
166  vout.general(m_vl, " Omega_tolerance = %8.2e\n", Omega_tolerance);
167 
168  //- range check
169  // NB. Omega_tolerance == 0.0 is allowed.
170 
171  //- store values
172  m_Omega_tolerance = Omega_tolerance;
173 }
174 
175 
176 //====================================================================
177 template<typename AFIELD>
178 void ASolver_BiCGStab_Cmplx<AFIELD>::solve(AFIELD& xq, const AFIELD& b,
179  int& Nconv, real_t& diff)
180 {
181  vout.paranoiac(m_vl, "%s: solver start:\n", class_name.c_str());
182 
183 #pragma omp master
184  {
185  m_init_mult = 0;
186  }
187 
188  copy(m_t, b);
189  real_t snorm2 = norm2(b);
190  real_t snorm = sqrt(snorm2);
191  real_t snorm_inv = real_t(1.0) / snorm;
192  vout.detailed(m_vl, " |b| = snorm = %22.15e\n", snorm);
193  scal(m_t, snorm_inv);
194 
195  // solutoin vector
196  m_px = &xq;
197  AFIELD& m_x = *m_px;
198 
199  // coefficents
200  coeff_t coeff;
201 
202  // initialication
203  int nconv = -1;
204  real_t rr;
205  solve_init(m_t, rr, coeff, snorm_inv, m_initial_mode);
206  vout.detailed(m_vl, " init: %22.15e\n", rr);
207 
208  for (int iter = 0; iter < m_Niter; ++iter) {
209  int iflg = 0;
210  solve_step(rr, iflg, coeff);
211  vout.detailed(m_vl, "%6d %22.15e\n", iter, rr);
212  if (iflg != 0) { // restart
213  copy(m_t, b);
214  scal(m_t, snorm_inv);
215  solve_init_GIVEN(m_t, rr, coeff, real_t(1.0));
216  vout.detailed(m_vl, "%s: restarted.\n", class_name.c_str());
217  vout.detailed(m_vl, "%6d %22.15e\n", iter, rr);
218  }
219 
220  if (rr < m_Stop_cond) {
221  nconv = 2 * (iter + 1); // counting fermion mult
222  break;
223  }
224  }
225  if (nconv == -1) {
226  if ((m_Stop_cond < 0) || (Nconv < 0)) {
227  vout.detailed(m_vl, " truncating with iter= %d\n", 2 * m_Niter);
228  nconv = 2 * m_Niter;
229  } else {
230  vout.crucial(m_vl, "Error at %s: not converged\n",
231  class_name.c_str());
232  vout.crucial(m_vl, " iter(final): %8d %22.15e\n",
233  m_Niter, rr * snorm2);
234  exit(EXIT_FAILURE);
235  }
236  } else {
237  vout.detailed(m_vl, "converged:\n");
238  }
239  vout.detailed(m_vl, " nconv = %d\n", nconv);
240 
241  scal(xq, snorm);
242 
243  m_fopr->mult(m_t, xq);
244  sub(m_t, b);
245  real_t diff2 = norm2(m_t);
246 
247 #pragma omp master
248  {
249  m_nconv = nconv; // include solve_init();
250  }
251  Nconv = nconv + m_init_mult; // include solve_init();
252  diff = sqrt(diff2 / snorm2);
253 
254 #pragma omp barrier
255 }
256 
257 
258 //====================================================================
259 template<typename AFIELD>
260 inline void ASolver_BiCGStab_Cmplx<AFIELD>::solve_init(const AFIELD& b,
261  real_t& rr,
262  coeff_t& prev,
263  const real_t scale,
264  const InitialGuess init_mode)
265 {
266  if (init_mode == InitialGuess::RHS) {
267  solve_init_RHS(b, rr, prev);
268  } else if (init_mode == InitialGuess::GIVEN) {
269  solve_init_GIVEN(b, rr, prev, scale);
270  } else if (init_mode == InitialGuess::ZERO) {
271  solve_init_ZERO(b, rr, prev);
272  } else {
273  vout.crucial("%s: unkown init guess mode\n", class_name.c_str());
274  exit(EXIT_FAILURE);
275  }
276 }
277 
278 
279 //====================================================================
280 template<typename AFIELD>
281 void ASolver_BiCGStab_Cmplx<AFIELD>::solve_init_RHS(const AFIELD& b,
282  real_t& rr,
283  coeff_t& prev)
284 {
285  AFIELD& m_x = *m_px;
286  copy(m_x, b);
287  copy(m_r, b);
288 
289  {
290  real_t temp = norm2(m_x);
291  vout.detailed(m_vl, " |m_x|^2 = %23.16e [init_RHS]\n", temp);
292  }
293 
294  m_fopr->mult(m_v, m_x);
295  axpy(m_r, real_t(-1.0), m_v);
296 
297  {
298  real_t temp = norm2(m_v);
299  vout.detailed(m_vl, " |m_v|^2 = %23.16e [init_RHS]\n", temp);
300  }
301 
302  copy(m_rh, m_r);
303  rr = norm2(m_r);
304  m_p.set(real_t(0.0));
305  m_v.set(real_t(0.0));
306 
307  prev.rho = real_t(1.0);
308  prev.alpha = real_t(1.0);
309  prev.omega = real_t(1.0);
310 
311  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx rr = %23.16e [init_RHS]\n", rr);
312 
313 
314 #pragma omp master
315  {
316  m_init_mult++;
317  }
318 }
319 
320 
321 //====================================================================
322 template<typename AFIELD>
323 void ASolver_BiCGStab_Cmplx<AFIELD>::solve_init_GIVEN(const AFIELD& b,
324  real_t& rr,
325  coeff_t& prev,
326  const real_t scale)
327 {
328  AFIELD& m_x = *m_px;
329  if (scale != real_t(1.0)) {
330  scal(m_x, scale);
331  }
332  copy(m_r, b);
333  m_fopr->mult(m_v, m_x);
334  axpy(m_r, real_t(-1.0), m_v);
335 
336  copy(m_rh, m_r);
337  rr = norm2(m_r);
338  m_p.set(real_t(0.0));
339  m_v.set(real_t(0.0));
340 
341  prev.rho = real_t(1.0);
342  prev.alpha = real_t(1.0);
343  prev.omega = real_t(1.0);
344 
345 #pragma omp master
346  {
347  m_init_mult++;
348  }
349  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx rr = %23.16e [init_GIVEN]\n", rr);
350 }
351 
352 
353 //====================================================================
354 template<typename AFIELD>
355 void ASolver_BiCGStab_Cmplx<AFIELD>::solve_init_ZERO(const AFIELD& b,
356  real_t& rr,
357  coeff_t& prev)
358 {
359  // x0 = 0
360  AFIELD& m_x = *m_px;
361  m_x.set(0.0);
362  copy(m_r, b);
363  copy(m_rh, m_r);
364  rr = norm2(m_r);
365  m_p.set(real_t(0.0));
366  m_v.set(real_t(0.0));
367 
368  prev.rho = real_t(1.0);
369  prev.alpha = real_t(1.0);
370  prev.omega = real_t(1.0);
371  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx rr = %23.16e [init_ZERO]\n", rr);
372 }
373 
374 
375 //====================================================================
376 template<typename AFIELD>
377 void ASolver_BiCGStab_Cmplx<AFIELD>::solve_step(real_t& rr, int& iflg, coeff_t& prev)
378 {
379  using complex_t = typename AFIELD::complex_t;
380  const real_t Omega_tolerance = m_Omega_tolerance;
381  AFIELD& m_x = *m_px;
382 
383  complex_t rho = dotc(m_rh, m_r);
384 #ifdef DEBUG_BICGSTAB_CMPLX
385  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx rho = %23.16e %23.16e\n", real(rho), imag(rho));
386 #endif
387  if (abs(rho) < m_ecrit) {
388  iflg = 1;
389  vout.detailed(m_vl, "too small value of rho = %23.16e %23.16e\n", real(rho), imag(rho));
390 #ifdef DEBUG
391  double r2 = norm2(m_r);
392  double rh2 = norm2(m_rh);
393  vout.detailed(m_vl, " |m_r|^2 = %23.16e\n", r2);
394  vout.detailed(m_vl, " |m_rh|^2 = %23.16e\n", rh2);
395 #endif
396  return;
397  }
398  complex_t bet = rho * prev.alpha / (prev.rho * prev.omega);
399 
400  // p = r + bet * (p - omega_prev * v);
401  axpy(m_p, -prev.omega, m_v);
402  aypx(bet, m_p, m_r);
403 
404  m_fopr->mult(m_v, m_p);
405 
406  complex_t aden = dotc(m_rh, m_v);
407 #ifdef DEBUG_BICGSTAB_CMPLX
408  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx bet = %23.16e %23.16e\n", real(bet), imag(bet));
409  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx aden = %23.16e %23.16e\n", real(aden), imag(aden));
410 #endif
411  if (abs(aden) < m_ecrit) {
412  iflg = 1;
413  vout.detailed(m_vl, "too small denominator aden = %23.16e %23.16e\n", real(aden), imag(aden));
414  return;
415  }
416  complex_t alpha = rho / aden;
417 #ifdef DEBUG_BICGSTAB_CMPLX
418  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx alpha = %23.16e %23.16e\n", real(alpha), imag(alpha));
419 #endif
420 
421  // "r" from here is an alias of "s"
422  axpy(m_r, -alpha, m_v);
423  m_fopr->mult(m_t, m_r);
424 
425  // todo: use something like dotc_and_norm2()
426  // complex_t omega_n = dotc(m_t, m_r);
427  // real_t omega_d = m_t.norm2();
428  complex_t omega_n;
429  real_t omega_d, s_norm2;
430  dotc_and_norm2(omega_n, omega_d, s_norm2, m_t, m_r);
431  complex_t omega = omega_n / omega_d;
432 
433  // for omega prescription
434  // const real_t s_norm2=m_r.norm2();
435  const double abs_rho = abs(omega_n) / sqrt(omega_d * s_norm2);
436 
437 #ifdef DEBUG_BICGSTAB_CMPLX
438  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx omega_n = %23.16e %23.16e\n", real(omega_n), imag(omega_n));
439  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx omega_d = %23.16e\n", omega_d);
440  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx omega = %23.16e %23.16e\n", real(omega), imag(omega));
441  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx s_norm2 = %23.16e\n", s_norm2);
442  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx abs_rho = %23.16e\n", abs_rho);
443 #endif
444 
445 
446  //- a prescription to improve stability of BiCGStab_Cmplx
447  if (abs_rho < Omega_tolerance) {
448  vout.detailed(m_vl, "ASolver_BiCGStab_Cmplx abs_rho = %23.16e, omega = %23.16e %23.16e [Omega presciption applied]\n", abs_rho, real(omega), imag(omega));
449  omega *= Omega_tolerance / abs_rho;
450  }
451 
452  // x = x + omega "s" + alpha p
453  axpy(m_x, omega, m_r);
454  axpy(m_x, alpha, m_p);
455 
456  // "r" below is the residual vector
457  axpy(m_r, -omega, m_t);
458 
459  prev.rho = rho;
460  prev.alpha = alpha;
461  prev.omega = omega;
462  rr = norm2(m_r);
463 }
464 
465 
466 //====================================================================
467 template<typename AFIELD>
468 double ASolver_BiCGStab_Cmplx<AFIELD>::flop_count()
469 {
470  int Nin = m_fopr->field_nin();
471  int Nvol = m_fopr->field_nvol();
472  int Nex = m_fopr->field_nex();
473  int NPE = CommonParameters::NPE();
474 
475  int ninit = 1;
476  int iter = m_nconv;
477  // if(m_nconv<0){ iter = m_Niter;}
478  double flop_field = static_cast<double>(Nin * Nvol * Nex) * NPE;
479  double flop_vector = (6 + 4 * ninit + iter * 20) * flop_field;
480 
481  double flop_fopr = (1 + ninit + iter) * m_fopr->flop_count();
482 
483  if (m_initial_mode == InitialGuess::ZERO) {
484  flop_vector -= (2 * ninit) * flop_field;
485  flop_fopr -= m_fopr->flop_count();
486  }
487 
488  double flop = flop_vector + flop_fopr;
489 
490  return flop;
491 }
492 
493 
494 //====================================================================
495 //============================================================END=====
ASolver_BiCGStab_Cmplx::solve_init_ZERO
void solve_init_ZERO(const AFIELD &b, real_t &rr, coeff_t &)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:355
ASolver_BiCGStab_Cmplx::coeff_t::omega
std::complex< real_t > omega
Definition: asolver_BiCGStab_Cmplx.h:79
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Parameters
Class for parameters.
Definition: parameters.h:46
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:281
aypx
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
Definition: field.cpp:509
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
dotc_and_norm2
void dotc_and_norm2(dcomplex &yx, double &y2, double &x2, const Field &y, const Field &x)
Definition: field.cpp:806
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
norm2
REALTYPE norm2(const AField< REALTYPE, QXS > &v)
Definition: afield-inc.h:453
ASolver_BiCGStab_Cmplx::solve_init_RHS
void solve_init_RHS(const AFIELD &b, real_t &rr, coeff_t &)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:281
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
Field::real_t
double real_t
Definition: field.h:51
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:300
ASolver_BiCGStab_Cmplx::solve
void solve(AFIELD &xq, const AFIELD &b, int &nconv, real_t &diff)
solver main.
Definition: asolver_BiCGStab_Cmplx-tmpl.h:178
Field::norm2
double norm2() const
Definition: field.cpp:113
ASolver_BiCGStab_Cmplx::solve_step
void solve_step(real_t &rr, int &iflg, coeff_t &)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:377
ASolver_BiCGStab_Cmplx::coeff_t
Definition: asolver_BiCGStab_Cmplx.h:75
ASolver_BiCGStab_Cmplx::set_parameters_BiCGStab_series
void set_parameters_BiCGStab_series(const real_t Omega_tolerance)
setting BiCGStab specific parameters.
Definition: asolver_BiCGStab_Cmplx-tmpl.h:161
ASolver_BiCGStab_Cmplx::init
void init(void)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:39
ASolver_BiCGStab_Cmplx
Definition: asolver_BiCGStab_Cmplx.h:52
ASolver_BiCGStab_Cmplx::tidyup
void tidyup(void)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:65
ASolver_BiCGStab_Cmplx::solve_init_GIVEN
void solve_init_GIVEN(const AFIELD &b, real_t &rr, coeff_t &, const real_t)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:323
dotc
dcomplex dotc(const Field &y, const Field &x)
Definition: field.cpp:712
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
ASolver_BiCGStab_Cmplx::solve_init
void solve_init(const AFIELD &b, real_t &rr, coeff_t &, const real_t scale, const InitialGuess)
Definition: asolver_BiCGStab_Cmplx-tmpl.h:260
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
ASolver::InitialGuess
InitialGuess
Definition: asolver.h:31
complex_t
ComplexTraits< double >::complex_t complex_t
Definition: afopr_Clover_coarse_double.cpp:23
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
Field
Container of Field-type object.
Definition: field.h:46
ASolver_BiCGStab_Cmplx::coeff_t::rho
std::complex< real_t > rho
Definition: asolver_BiCGStab_Cmplx.h:77
ASolver_BiCGStab_Cmplx::coeff_t::alpha
std::complex< real_t > alpha
Definition: asolver_BiCGStab_Cmplx.h:78
ASolver_BiCGStab_Cmplx::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
Definition: asolver_BiCGStab_Cmplx-tmpl.h:74
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
ASolver::real_t
AFIELD::real_t real_t
Definition: asolver.h:29
ASolver_BiCGStab_Cmplx::flop_count
double flop_count()
returns the floating point operation count.
Definition: asolver_BiCGStab_Cmplx-tmpl.h:468
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:569