Bridge++  Ver. 2.0.4
asolver_FBiCGStab-tmpl.h
Go to the documentation of this file.
1 
14 #include "asolver_FBiCGStab.h"
15 
16 #include "lib/ResourceManager/threadManager.h"
17 
18 
19 template<typename AFIELD>
20 const std::string ASolver_FBiCGStab<AFIELD>::class_name
21  = "ASolver_FBiCGStab";
22 //====================================================================
23 template<typename AFIELD>
24 void ASolver_FBiCGStab<AFIELD>::init()
25 {
26  ThreadManager::assert_single_thread(class_name);
27 
28  int nin = m_fopr->field_nin();
29  int nvol = m_fopr->field_nvol();
30  int nex = m_fopr->field_nex();
31 
32  m_r.reset(nin, nvol, nex);
33  m_x.reset(nin, nvol, nex);
34  m_p.reset(nin, nvol, nex);
35  m_s.reset(nin, nvol, nex);
36  m_v.reset(nin, nvol, nex);
37  m_t.reset(nin, nvol, nex);
38  m_rh.reset(nin, nvol, nex);
39  m_u.reset(nin, nvol, nex);
40 
41  m_nconv = -1;
42  m_Omega_tolerance = 0.7;
43 }
44 
45 
46 //====================================================================
47 template<typename AFIELD>
48 void ASolver_FBiCGStab<AFIELD>::tidyup()
49 {
50  // ThreadManager::assert_single_thread(class_name);
51  // nothing is to be deleted.
52 }
53 
54 
55 //====================================================================
56 template<typename AFIELD>
57 void ASolver_FBiCGStab<AFIELD>::set_parameters(
58  const Parameters& params)
59 {
60  const string str_vlevel = params.get_string("verbose_level");
61 
62  m_vl = vout.set_verbose_level(str_vlevel);
63 
64  //- fetch and check input parameters
65  int Niter, Nrestart;
66  double Stop_cond;
67 
68  int err = 0;
69  err += params.fetch_int("maximum_number_of_iteration", Niter);
70  err += params.fetch_int("maximum_number_of_restart", Nrestart);
71  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
72 
73  if (err) {
74  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
75  class_name.c_str());
76  exit(EXIT_FAILURE);
77  }
78 
79  double Omega_tolerance;
80  err = params.fetch_double("Omega_tolerance", Omega_tolerance);
81  if (err) {
82  Omega_tolerance = 0.7;
83  }
84 
85  // "restart" is not yet supported
86  int Niter2 = Niter * Nrestart;
87  set_parameters(Niter2, Stop_cond);
88 
89  set_parameters_BiCGStab_series(Omega_tolerance);
90 }
91 
92 
93 //====================================================================
94 template<typename AFIELD>
95 void ASolver_FBiCGStab<AFIELD>::set_parameters(const int Niter,
96  const real_t Stop_cond,
97  const InitialGuess init_guess_mode)
98 {
99  ThreadManager::assert_single_thread(class_name);
100 
101  m_Niter = Niter;
102  m_Stop_cond = Stop_cond;
103  m_initial_mode = init_guess_mode;
104  std::string prec = "double";
105  if (sizeof(real_t) == 4) prec = "float";
106 
107  vout.general(m_vl, "%s:\n", class_name.c_str());
108  vout.general(m_vl, " Precision: %s\n", prec.c_str());
109  vout.general(m_vl, " Niter = %d\n", m_Niter);
110  vout.general(m_vl, " Stop_cond = %16.8e\n", m_Stop_cond);
111  vout.general(m_vl, " InitialGuess = %d\n", m_initial_mode);
112 }
113 
114 
115 //====================================================================
116 template<typename AFIELD>
117 void ASolver_FBiCGStab<AFIELD>::set_parameters(const int Niter,
118  const real_t Stop_cond)
119 {
120  set_parameters(Niter, Stop_cond, InitialGuess::RHS);
121 }
122 
123 
124 //====================================================================
125 template<typename AFIELD>
126 void ASolver_FBiCGStab<AFIELD>::set_parameters(const int Niter,
127  const real_t Stop_cond,
128  const bool use_init_guess)
129 {
130  // for backward compatibility
131  if (use_init_guess) {
132  set_parameters(Niter, Stop_cond, InitialGuess::GIVEN);
133  } else {
134  set_parameters(Niter, Stop_cond, InitialGuess::RHS);
135  }
136 }
137 
138 
139 //====================================================================
140 template<typename AFIELD>
141 void ASolver_FBiCGStab<AFIELD>::set_parameters_BiCGStab_series(const real_t Omega_tolerance)
142 {
143  ThreadManager::assert_single_thread(class_name);
144 
145  //- print input parameters
146  vout.general(m_vl, " Omega_tolerance = %8.2e\n", Omega_tolerance);
147 
148  //- range check
149  // NB. Omega_tolerance == 0.0 is allowed.
150 
151  //- store values
152  m_Omega_tolerance = Omega_tolerance;
153 }
154 
155 
156 //====================================================================
157 template<typename AFIELD>
158 void ASolver_FBiCGStab<AFIELD>::solve(AFIELD& xq,
159  const AFIELD& b,
160  int& Nconv, real_t& diff)
161 {
162  vout.paranoiac(m_vl, "%s: solver start.\n", class_name.c_str());
163 
164 
165  real_t bnorm2, bnorm, bnorm_inv, bnorm2_inv;
166  real_t rr;
167  coeff_t prev;
168 
169  // initial barrier
170 #pragma omp barrier
171 
172  bnorm2 = norm2(b);
173  if (!(bnorm2 > 0.0)) {
174  xq.set(0.0);
175  return;
176  }
177  bnorm = sqrt(bnorm2);
178  bnorm_inv = real_t(1.0) / bnorm;
179  bnorm2_inv = real_t(1.0) / bnorm2;
180 
181  m_prec->reset_flop_count();
182 
183  int nconv = -1;
184 
185  solve_init(b, xq, rr, prev);
186  vout.detailed(m_vl, " FBiCGStab init: %22.15e\n", rr * bnorm2_inv);
187  assert(rr > 0.0);
188 
189  int iter2 = 0;
190  for (int iter = 0; iter < m_Niter; ++iter) {
191  solve_step1(rr, prev);
192  iter2++;
193  double flop = 0.0;
194 #pragma omp master
195  {
196  flop = flop_count_intermediate(iter2);
197  }
198  vout.detailed(m_vl, "FBiCGStab %d %22.15e Flop(double+float) %e\n", iter2, rr * bnorm2_inv, flop);
199  if (rr * bnorm2_inv < m_Stop_cond) {
200  nconv = iter2; // counting fermion mult
201  break;
202  }
203  solve_step2(rr, prev);
204  iter2++;
205 #pragma omp master
206  {
207  flop = flop_count_intermediate(iter2);
208  }
209  vout.detailed(m_vl, "FBiCGStab %d %22.15e Flop(double+float) %e\n", iter2, rr * bnorm2_inv, flop);
210  if (rr * bnorm2_inv < m_Stop_cond) {
211  nconv = iter2; // counting fermion mult
212  break;
213  }
214  }
215 
216  if (nconv == -1) {
217  vout.detailed(m_vl, "FBiCGStab NOT converged:\n");
218  // vout.crucial(m_vl, "Error at %s: not converged\n",
219  // class_name.c_str());
220  // vout.crucial(m_vl, " iter(final): %8d %22.15e\n",
221  // m_Niter, rr*snorm);
222  // exit(EXIT_FAILURE);
223  }
224 
225  vout.detailed(m_vl, "FBiCGStab converged:\n");
226 
227  copy(xq, m_x);
228 
229  // check the solution
230  m_fopr->mult(m_s, xq);
231  axpy(m_s, real_t(-1.0), b);
232  real_t diff2 = norm2(m_s);
233  vout.general(m_vl, "FBiCGStab nconv = %d, diff2 = %22.15e\n", nconv, diff2 * bnorm2_inv);
234 
235 
236 #pragma omp master
237  {
238  m_nconv = nconv;
239  Nconv = m_nconv + 1; // include solve_init();
240  diff = sqrt(diff2 * bnorm2_inv);
241  }
242 #pragma omp barrier
243 }
244 
245 
246 //====================================================================
247 template<typename AFIELD>
248 void ASolver_FBiCGStab<AFIELD>::solve_init(const AFIELD& b,
249  const AFIELD& xq,
250  real_t& rr,
251  coeff_t& prev)
252 {
253  // initial guess and its residual
254  if (m_initial_mode == InitialGuess::RHS) {
255  copy(m_r, b);
256  copy(m_x, b);
257  m_fopr->mult(m_v, m_x);
258  axpy(m_r, real_t(-1.0), m_v);
259  } else if (m_initial_mode == InitialGuess::GIVEN) {
260  copy(m_r, b);
261  copy(m_x, xq);
262  m_fopr->mult(m_v, m_x);
263  axpy(m_r, real_t(-1.0), m_v);
264  } else if (m_initial_mode == InitialGuess::ZERO) {
265  copy(m_r, b);
266  m_x.set(0.0);
267  } else {
268  vout.crucial("%s: unkown init guess mode\n", class_name.c_str());
269  exit(EXIT_FAILURE);
270  }
271 
272  copy(m_rh, m_r);
273 
274  rr = norm2(m_r);
275  assert(rr > 0.0);
276 
277  copy(m_p, m_r);
278 
279  // initial parameters
280  prev.rho = rr;
281  prev.alpha = 0.0;
282 
283  vout.detailed(m_vl, "ASolver_FBiCGStab rr = %23.16e [init]\n", rr);
284 }
285 
286 
287 //====================================================================
288 template<typename AFIELD>
289 void ASolver_FBiCGStab<AFIELD>::solve_step1(real_t& rr, coeff_t& prev)
290 {
291  using complex_t = std::complex<real_t>;
292 
293  m_prec->mult(m_u, m_p);
294  m_fopr->mult(m_v, m_u); // m_v = D M p
295 
296  complex_t alpha_d = dotc(m_rh, m_v);
297  complex_t alpha = prev.rho / alpha_d;
298  prev.alpha = alpha;
299 
300  axpy(m_r, -alpha, m_v);
301  axpy(m_x, alpha, m_u);
302  rr = norm2(m_r);
303 
304 #ifdef DEBUG_FBICGSTAB
305  vout.detailed(m_vl, "ASolver_FBiCGStab alpha_d = %23.16e %23.16e\n", real(alpha_d), std::imag(alpha_d));
306  vout.detailed(m_vl, "ASolver_FBiCGStab alpha = %23.16e %23.16e\n", real(alpha), std::imag(alpha));
307  vout.detailed(m_vl, "ASolver_FBiCGStab rr = %23.16e\n", rr);
308 #endif
309 }
310 
311 
312 //====================================================================
313 template<typename AFIELD>
314 void ASolver_FBiCGStab<AFIELD>::solve_step2(real_t& rr, coeff_t& prev)
315 {
316  using complex_t = std::complex<real_t>;
317  m_prec->mult(m_u, m_r);
318  m_fopr->mult(m_t, m_u);
319 
320  complex_t omega_n = dotc(m_t, m_r); // omega_n = t * s;
321  real_t omega_d = norm2(m_t); // omega_d = t * t;
322  real_t r_norm2 = norm2(m_r);
323  complex_t omega = omega_n / omega_d;
324 
325  // for omega prescrition
326  const real_t Omega_tolerance = m_Omega_tolerance;
327  const real_t abs_rho = abs(omega_n) / sqrt(omega_d * r_norm2);
328 
329 #ifdef DEBUG_FBICGSTAB
330  vout.detailed(m_vl, "ASolver_FBiCGStab omega_n = %23.16e %23.16e\n", real(omega_n), std::imag(omega_n));
331  vout.detailed(m_vl, "ASolver_FBiCGStab omega_d = %23.16e\n", omega_d);
332  vout.detailed(m_vl, "ASolver_FBiCGStab omega = %23.16e %23.16e\n", real(omega), std::imag(omega));
333  vout.detailed(m_vl, "ASolver_FBiCGStab r_norm2 = %23.16e\n", r_norm2);
334  vout.detailed(m_vl, "ASolver_FBiCGStab abs_rho = %23.16e\n", abs_rho);
335 #endif
336 
337  //- a prescription to improve stability of BiCGStab_Cmplx
338  if (abs_rho < Omega_tolerance) {
339  omega *= Omega_tolerance / abs_rho;
340  vout.detailed(m_vl, "ASolver_FBiCGStab omega = %23.16e %23.16e [Omega presciption applied]\n", real(omega), std::imag(omega));
341  }
342 
343  axpy(m_x, omega, m_u);
344  axpy(m_r, -omega, m_t);
345 
346  complex_t rho = dotc(m_rh, m_r);
347  complex_t beta = (rho / prev.rho) * (prev.alpha * omega);
348  prev.rho = rho;
349 
350  // p = r + bet * (p - omega * v);
351  axpy(m_p, -omega, m_v); // p := p - omega v
352  aypx(beta, m_p, m_r);
353 
354  rr = norm2(m_r); // rr = r * r;
355 #ifdef DEBUG_FBICGSTAB
356  vout.detailed(m_vl, "ASolver_FBiCGStab rho = %23.16e %23.16e\n", real(rho), std::imag(rho));
357  vout.detailed(m_vl, "ASolver_FBiCGStab beta = %23.16e %23.16e\n", real(beta), std::imag(beta));
358  vout.detailed(m_vl, "ASolver_FBiCGStab rr = %23.16e\n", rr);
359 #endif
360 }
361 
362 
363 //====================================================================
364 template<typename AFIELD>
365 void ASolver_FBiCGStab<AFIELD>::prec(AFIELD& v, AFIELD& w)
366 {
367  copy(v, w);
368 }
369 
370 
371 //====================================================================
372 template<typename AFIELD>
373 double ASolver_FBiCGStab<AFIELD>::flop_count()
374 {
375  return flop_count_intermediate(m_nconv);
376 }
377 
378 
379 //====================================================================
380 template<typename AFIELD>
381 double ASolver_FBiCGStab<AFIELD>::flop_count_intermediate(const int iter)
382 {
383  int Nin = m_fopr->field_nin();
384  int Nvol = m_fopr->field_nvol();
385  int Nex = m_fopr->field_nex();
386  int NPE = CommonParameters::NPE();
387 
388  int ninit = 1;
389 
390  double flop_field = static_cast<double>(Nin * Nvol * Nex) * NPE;
391  // step1: 4+4+4+2 = 14
392  // step2: 4+2+4+4+4+4+4+2 = 4*6 + 2*2= 28
393  double flop_vector = (6 + ninit * 4 + (iter / 2) * 42) * flop_field;
394  if (iter % 2 == 1) {
395  flop_vector += 14 * flop_field;
396  }
397 
398  double flop_fopr = (1 + ninit + iter) * m_fopr->flop_count();
399  double flop_prec = m_prec->flop_count();
400 
401  double flop = flop_vector + flop_fopr + flop_prec;
402 
403  return flop;
404 }
405 
406 
407 //====================================================================
ASolver_FBiCGStab::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
Definition: asolver_FBiCGStab-tmpl.h:57
ASolver_FBiCGStab
Definition: asolver_FBiCGStab.h:32
ASolver_FBiCGStab::flop_count_intermediate
double flop_count_intermediate(const int iter)
Definition: asolver_FBiCGStab-tmpl.h:381
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
ASolver_FBiCGStab::solve_step1
void solve_step1(real_t &rr, coeff_t &prev)
Definition: asolver_FBiCGStab-tmpl.h:289
asolver_FBiCGStab.h
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
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_FBiCGStab::solve_init
void solve_init(const AFIELD &b, const AFIELD &xq, real_t &rr, coeff_t &prev)
Definition: asolver_FBiCGStab-tmpl.h:248
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:300
ASolver_FBiCGStab::init
void init(void)
initial setup.
Definition: asolver_FBiCGStab-tmpl.h:24
ASolver_FBiCGStab::flop_count
double flop_count()
returns the floating point operation count.
Definition: asolver_FBiCGStab-tmpl.h:373
ASolver_FBiCGStab::coeff_t::alpha
std::complex< real_t > alpha
Definition: asolver_FBiCGStab.h:43
ASolver_FBiCGStab::solve_step2
void solve_step2(real_t &rr, coeff_t &prev)
Definition: asolver_FBiCGStab-tmpl.h:314
threadManager.h
dotc
dcomplex dotc(const Field &y, const Field &x)
Definition: field.cpp:712
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
ASolver_FBiCGStab::tidyup
void tidyup(void)
final tidy-up.
Definition: asolver_FBiCGStab-tmpl.h:48
ASolver_FBiCGStab::coeff_t::rho
std::complex< real_t > rho
Definition: asolver_FBiCGStab.h:42
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
ASolver_FBiCGStab::prec
void prec(AFIELD &, AFIELD &)
Definition: asolver_FBiCGStab-tmpl.h:365
ASolver_FBiCGStab::solve
void solve(AFIELD &xq, const AFIELD &b, int &nconv, real_t &diff)
solver main.
Definition: asolver_FBiCGStab-tmpl.h:158
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
ASolver_FBiCGStab::set_parameters_BiCGStab_series
void set_parameters_BiCGStab_series(const real_t Omega_tolerance)
setting BiCGStab specific parameters.
Definition: asolver_FBiCGStab-tmpl.h:141
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
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_FBiCGStab::coeff_t
Definition: asolver_FBiCGStab.h:40
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