Bridge++  Ver. 2.0.2
asolver_BiCGStab-tmpl.h
Go to the documentation of this file.
1 /*
2  @file asolver_BiCGStab-tmpl.h
3  @brief
4  @author Hideo Matsufuru (matufuru)
5  $LastChangedBy: matufuru $
6  @date $LastChangedDate:: 2#$
7  @version $LastChangedRevision: 2492 $
8 */
9 
10 template<typename AFIELD>
11 const std::string ASolver_BiCGStab<AFIELD>::class_name
12  = "ASolver_BiCGStab";
13 //#define DEBUG_BICGSTAB
14 //====================================================================
15 template<typename AFIELD>
16 void ASolver_BiCGStab<AFIELD>::init(void)
17 {
18  ThreadManager::assert_single_thread(class_name);
19 
20  vout.detailed(m_vl, "%s: being setup.\n", class_name.c_str());
21 
22  int nin = m_fopr->field_nin();
23  int nvol = m_fopr->field_nvol();
24  int nex = m_fopr->field_nex();
25  vout.detailed(m_vl, " Field size: Nin = %d Nvol = %d Nex = %d\n",
26  nin, nvol, nex);
27 
28  m_r.reset(nin, nvol, nex);
29  m_x.reset(nin, nvol, nex);
30  m_p.reset(nin, nvol, nex);
31  m_s.reset(nin, nvol, nex);
32  m_v.reset(nin, nvol, nex);
33  m_t.reset(nin, nvol, nex);
34  m_rh.reset(nin, nvol, nex);
35 
36  m_ecrit = 1.e-32;
37  m_nconv = -1;
38  m_Omega_tolerance = 0.7;
39 
40  m_initial_mode = InitialGuess::RHS;
41 
42  vout.detailed(m_vl, "%s: setup finished.\n", class_name.c_str());
43 }
44 
45 
46 //====================================================================
47 template<typename AFIELD>
48 void ASolver_BiCGStab<AFIELD>::tidyup(void)
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_BiCGStab<AFIELD>::set_parameters(const Parameters& params)
58 {
59  const string str_vlevel = params.get_string("verbose_level");
60 
61  m_vl = vout.set_verbose_level(str_vlevel);
62 
63  //- fetch and check input parameters
64  int Niter, Nrestart;
65  double Stop_cond;
66  double Omega_tolerance;
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  err = params.fetch_double("Omega_tolerance", Omega_tolerance);
79  if (err) {
80  Omega_tolerance = 0.7;
81  }
82 
83  // - default initial mode is RHS
84  // - negative Stop_cond is for fixed nubmer of iterations
85  int Niter2 = Niter * Nrestart;
86  set_parameters(Niter2, Stop_cond, InitialGuess::RHS);
87  set_parameters_BiCGStab_series(Omega_tolerance);
88 }
89 
90 
91 //====================================================================
92 template<typename AFIELD>
93 void ASolver_BiCGStab<AFIELD>::set_parameters(const int Niter,
94  const real_t Stop_cond)
95 {
96  set_parameters(Niter, Stop_cond, InitialGuess::RHS);
97 }
98 
99 
100 //====================================================================
101 template<typename AFIELD>
102 void ASolver_BiCGStab<AFIELD>::set_parameters(const int Niter,
103  const real_t Stop_cond,
104  const InitialGuess init_guess_mode)
105 
106 {
107  ThreadManager::assert_single_thread(class_name);
108 
109  m_Niter = Niter;
110  m_Stop_cond = Stop_cond;
111  std::string prec = "double";
112  if (sizeof(real_t) == 4) prec = "float";
113  m_initial_mode = init_guess_mode;
114 
115  vout.general(m_vl, "%s:\n", class_name.c_str());
116  vout.general(m_vl, " Precision: %s\n", prec.c_str());
117  vout.general(m_vl, " Niter = %d\n", m_Niter);
118  vout.general(m_vl, " Stop_cond = %16.8e\n", m_Stop_cond);
119  vout.general(m_vl, " init_guess_mode: %d\n", m_initial_mode);
120 }
121 
122 
123 //====================================================================
124 template<typename AFIELD>
125 void ASolver_BiCGStab<AFIELD>::set_parameters_BiCGStab_series(const real_t Omega_tolerance)
126 {
127  ThreadManager::assert_single_thread(class_name);
128 
129  //- print input parameters
130  vout.general(m_vl, " Omega_tolerance = %8.2e\n", Omega_tolerance);
131 
132  //- range check
133  // NB. Omega_tolerance == 0.0 is allowed.
134 
135  //- store values
136  m_Omega_tolerance = Omega_tolerance;
137 }
138 
139 
140 //====================================================================
141 template<typename AFIELD>
142 void ASolver_BiCGStab<AFIELD>::solve(AFIELD& xq, const AFIELD& b,
143  int& Nconv, real_t& diff)
144 {
145  vout.paranoiac(m_vl, "%s: solver start:\n", class_name.c_str());
146 #pragma omp master
147  {
148  m_init_mult = 0;
149  }
150  copy(m_s, b);
151  real_t sr = norm2(m_s);
152  real_t snorm = real_t(1.0) / sr;
153  real_t snorm2 = sqrt(snorm);
154  vout.detailed(m_vl, " sr = %22.15e\n", sr);
155  vout.detailed(m_vl, " snorm = %22.15e\n", snorm);
156 
157  real_t alpha_prev, rho_prev, omega_prev;
158  scal(m_s, snorm2);
159 
160  if (m_initial_mode == InitialGuess::RHS) {
161  {
162  vout.detailed(m_vl, "%s: initial_mode = RHS\n", class_name.c_str());
163 #ifdef DEBUG_BICGSTAB
164  double s2 = m_s.norm2();
165 
166  vout.detailed(m_vl, "%s: m_x: %d %d %d\n", class_name.c_str(), m_x.nin(), m_x.nvol(), m_x.nex());
167  vout.detailed(m_vl, "%s: m_s: %d %d %d: nrom2=%23.16e\n", class_name.c_str(), m_s.nin(), m_s.nvol(), m_s.nex(), s2);
168 #endif
169  }
170  copy(m_x, m_s);
171  } else if (m_initial_mode == InitialGuess::GIVEN) {
172  vout.detailed(m_vl, "%s: initial_mode = GIVEN\n", class_name.c_str());
173  copy(m_x, xq);
174  scal(m_x, snorm2);
175  } else if (m_initial_mode == InitialGuess::ZERO) {
176  vout.detailed(m_vl, "%s: initial_mode = ZERO\n", class_name.c_str());
177  } else {
178  vout.crucial("%s: unkown init guess mode\n", class_name.c_str());
179  exit(EXIT_FAILURE);
180  }
181 
182  int nconv = -1;
183 
184  real_t rr;
185  solve_init(b, rr, alpha_prev, rho_prev, omega_prev, m_initial_mode);
186  vout.detailed(m_vl, " init: %22.15e\n", rr);
187 
188  for (int iter = 0; iter < m_Niter; ++iter) {
189  int iflg = 0;
190  solve_step(rr, iflg, alpha_prev, rho_prev, omega_prev);
191 #ifdef DEBUG_BICGSTAB
192  {
193  m_fopr->mult(xq, m_x);
194  copy(m_t, b);
195  m_t.scal(snorm2);
196  axpy(xq, -1.0, m_t);
197  double r2 = xq.norm2();
198  vout.detailed(m_vl, " iter=%d, alg_rr=%22.15e, actual r2=%22.15e\n", iter, rr, r2);
199  }
200 #endif
201  if (iflg != 0) {
202  copy(m_s, b);
203  scal(m_s, snorm2);
204  solve_init(b, rr, alpha_prev, rho_prev, omega_prev, InitialGuess::GIVEN);
205  }
206  vout.detailed(m_vl, "%6d %22.15e\n", iter, rr);
207 
208  // if(rr*snorm < m_Stop_cond)
209  if (rr < m_Stop_cond) {
210  nconv = 2 * (iter + 1); // counting fermion mult
211  break;
212  }
213  }
214 
215  if (nconv == -1) {
216  if ((m_Stop_cond < 0) || (Nconv < 0)) {
217  vout.detailed(m_vl, " truncating with iter= %d\n", 2 * m_Niter);
218  nconv = 2 * m_Niter;
219  } else {
220  vout.crucial(m_vl, "Error at %s: not converged\n",
221  class_name.c_str());
222  vout.crucial(m_vl, " iter(final): %8d %22.15e\n",
223  m_Niter, rr);
224  exit(EXIT_FAILURE);
225  }
226  } else {
227  vout.detailed(m_vl, "converged:\n");
228  }
229  vout.detailed(m_vl, " nconv = %d\n", nconv);
230 
231  copy(xq, m_x);
232  real_t sr2 = sqrt(sr);
233  scal(xq, sr2);
234 
235  m_fopr->mult(m_s, xq);
236 
237  axpy(m_s, real_t(-1.0), b);
238  real_t diff2 = norm2(m_s);
239 
240 #pragma omp master
241  {
242  m_nconv = nconv;
243  Nconv = nconv + m_init_mult; // include solve_init();
244  diff = sqrt(diff2 / sr);
245  }
246 #pragma omp barrier
247 }
248 
249 
250 //====================================================================
251 template<typename AFIELD>
252 void ASolver_BiCGStab<AFIELD>::solve_init(const AFIELD& b,
253  real_t& rr,
254  real_t& alpha_prev, real_t& rho_prev, real_t& omega_prev,
255  const InitialGuess init_mode)
256 {
257  if (init_mode == InitialGuess::ZERO) {
258  m_x.set(0.0);
259  m_v.set(0.0);
260  copy(m_r, m_s);
261  } else {
262  m_fopr->mult(m_v, m_x);
263  copy(m_r, m_s);
264  axpy(m_r, real_t(-1.0), m_v);
265 #pragma omp master
266  {
267  m_init_mult++;
268  }
269  }
270 #ifdef DEBUG_BICGSTAB
271  {
272  double x2 = m_x.norm2();
273  double v2 = m_v.norm2();
274  vout.general(m_vl, "%s: x2=%23.15e\n", class_name.c_str(), x2);
275  vout.general(m_vl, "%s: v2=%23.15e\n", class_name.c_str(), v2);
276  }
277 #endif
278 
279 
280  copy(m_rh, m_r);
281  rr = norm2(m_r);
282 #pragma omp barrier
283  m_p.set(real_t(0.0));
284  m_v.set(real_t(0.0));
285 
286  rho_prev = real_t(1.0);
287  alpha_prev = real_t(1.0);
288  omega_prev = real_t(1.0);
289 }
290 
291 
292 //====================================================================
293 template<typename AFIELD>
294 void ASolver_BiCGStab<AFIELD>::solve_step(real_t& rr, int& iflg,
295  real_t& alpha_prev, real_t& rho_prev, real_t& omega_prev)
296 {
297  real_t rho = dot(m_rh, m_r);
298  if (fabs(rho) < m_ecrit) {
299  iflg = 1;
300  vout.detailed(m_vl, "too small value of rho = %e\n", rho);
301  return;
302  }
303  real_t bet = rho * alpha_prev / (rho_prev * omega_prev);
304 #ifdef DEBUG_BICGSTAB
305  vout.general(m_vl, "%s: rho = %23.15e\n", class_name.c_str(), rho);
306  vout.general(m_vl, "%s: alpha_prev = %23.15e\n", class_name.c_str(), alpha_prev);
307  vout.general(m_vl, "%s: rho_prev = %23.15e\n", class_name.c_str(), rho_prev);
308  vout.general(m_vl, "%s: omega_prev = %23.15e\n", class_name.c_str(), omega_prev);
309  vout.general(m_vl, "%s: bet = %23.15e\n", class_name.c_str(), bet);
310 #endif
311 
312  // p = r + bet * (p - omega_prev * v);
313  axpy(m_p, -omega_prev, m_v);
314  aypx(bet, m_p, m_r);
315 
316  m_fopr->mult(m_v, m_p);
317  real_t aden = dot(m_rh, m_v);
318 
319 #ifdef DEBUG_BICGSTAB
320  {
321  double p2 = m_p.norm2();
322  double v2 = m_v.norm2();
323  vout.general(m_vl, "%s: p2=%23.15e\n", class_name.c_str(), p2);
324  vout.general(m_vl, "%s: v2=%23.15e\n", class_name.c_str(), v2);
325  vout.general(m_vl, "%s: aden=%23.15e\n", class_name.c_str(), aden);
326  }
327 #endif
328 
329  if (fabs(aden) < m_ecrit) {
330  iflg = 1;
331  vout.detailed(m_vl, "too small denominator aden = %e\n", aden);
332  return;
333  }
334  real_t alpha = rho / aden;
335  // if(fabs(alpha) >100.0) {
336  // iflg = 1;
337  // vout.detailed(m_vl, "too large alpha = %e\n", alpha);
338  // return;
339  // }
340  copy(m_s, m_r);
341  axpy(m_s, -alpha, m_v);
342 
343  m_fopr->mult(m_t, m_s);
344 
345  real_t omega_n = dot(m_t, m_s);
346  real_t omega_d = m_t.norm2();
347  real_t s_norm2 = m_s.norm2();
348 #ifdef DEBUG_BICGSTAB
349  vout.general(m_vl, "%s: alpha = %23.15e\n", class_name.c_str(), alpha);
350  vout.general(m_vl, "%s: omega_n = %23.15e\n", class_name.c_str(), omega_n);
351  vout.general(m_vl, "%s: omega_d = %23.15e\n", class_name.c_str(), omega_d);
352  vout.general(m_vl, "%s: s_norm2 = %23.15e\n", class_name.c_str(), s_norm2);
353 #endif
354  if (fabs(omega_d) < m_ecrit) {
355  iflg = 1;
356  vout.detailed(m_vl, "too small denominator omega_d = %e\n", omega_d);
357  return;
358  }
359  real_t omega = omega_n / omega_d;
360  real_t abs_rho = abs(omega_n) / sqrt(omega_d * s_norm2);
361 #ifdef DEBUG_BICGSTAB
362  vout.general(m_vl, "%s: omega =%23.15e\n", class_name.c_str(), omega);
363  vout.general(m_vl, "%s: abs_rho=%23.15e\n", class_name.c_str(), abs_rho);
364 #endif
365  if (abs_rho < m_Omega_tolerance) {
366  vout.detailed(m_vl, "Omega_prescription: abs_rho = %23.16e, omega= %23.16e\n", abs_rho, omega);
367  omega *= m_Omega_tolerance / abs_rho;
368  }
369  axpy(m_x, omega, m_s);
370  axpy(m_x, alpha, m_p);
371 
372  copy(m_r, m_s);
373  axpy(m_r, -omega, m_t);
374 
375  rr = norm2(m_r);
376 #ifdef DEBUG_BICGSTAB
377  vout.general(m_vl, "%s: rr = %23.15e\n", class_name.c_str(), rr);
378 #endif
379  rho_prev = rho;
380  alpha_prev = alpha;
381  omega_prev = omega;
382 }
383 
384 
385 //====================================================================
386 template<typename AFIELD>
387 double ASolver_BiCGStab<AFIELD>::flop_count()
388 {
389  int Nin = m_fopr->field_nin();
390  int Nvol = m_fopr->field_nvol();
391  int Nex = m_fopr->field_nex();
392  int NPE = CommonParameters::NPE();
393 
394  int ninit = m_init_mult;
395 
396  double flop_field = static_cast<double>(Nin * Nvol * Nex) * NPE;
397  double flop_vector = (6 + ninit * 4 + m_nconv * 11) * flop_field;
398  double flop_fopr = (1 + ninit + m_nconv) * m_fopr->flop_count();
399 
400  if (m_initial_mode == InitialGuess::ZERO) {
401  flop_vector -= m_init_mult * 2 * flop_field;
402  }
403  double flop = flop_vector + flop_fopr;
404 
405  return flop;
406 }
407 
408 
409 //====================================================================
410 //============================================================END=====
ASolver_BiCGStab::tidyup
void tidyup(void)
Definition: asolver_BiCGStab-tmpl.h:48
Parameters
Class for parameters.
Definition: parameters.h:46
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
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
dot
double dot(const Field &y, const Field &x)
Definition: field.cpp:576
ASolver_BiCGStab
Definition: asolver_BiCGStab.h:15
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::set_parameters_BiCGStab_series
void set_parameters_BiCGStab_series(const real_t Omega_tolerance)
setting BiCGStab specific parameters.
Definition: asolver_BiCGStab-tmpl.h:125
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:238
Field::norm2
double norm2() const
Definition: field.cpp:113
ASolver_BiCGStab::solve
void solve(AFIELD &xq, const AFIELD &b, int &nconv, real_t &diff)
solver main.
Definition: asolver_BiCGStab-tmpl.h:142
ASolver_BiCGStab::init
void init(void)
Definition: asolver_BiCGStab-tmpl.h:16
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
ASolver_BiCGStab::solve_init
void solve_init(const AFIELD &b, real_t &rr, real_t &alpha_prev, real_t &rho_prev, real_t &omega_prev, const InitialGuess init_mode)
Definition: asolver_BiCGStab-tmpl.h:252
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
ASolver_BiCGStab::flop_count
double flop_count()
returns the floating point operation count.
Definition: asolver_BiCGStab-tmpl.h:387
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
ASolver_BiCGStab::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
Definition: asolver_BiCGStab-tmpl.h:57
ASolver::InitialGuess
InitialGuess
Definition: asolver.h:31
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:180
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:200
ASolver_BiCGStab::solve_step
void solve_step(real_t &rr, int &iflg, real_t &alpha_prev, real_t &rho_prev, real_t &omega_prev)
Definition: asolver_BiCGStab-tmpl.h:294
ASolver::real_t
AFIELD::real_t real_t
Definition: asolver.h:29
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:512