Bridge++  Ver. 2.0.2
fprop_alt_Standard_eo-tmpl.h
Go to the documentation of this file.
1 
10 template<typename AFIELD>
11 const std::string Fprop_alt_Standard_eo<AFIELD>::class_name
12  = "Fprop_alt_Standard_eo";
13 //====================================================================
14 template<typename AFIELD>
15 void Fprop_alt_Standard_eo<AFIELD>::init(const Parameters& params_fopr,
16  const Parameters& params_solver)
17 {
18  // this constructor assumes that the factories are available.
19  vout.general(m_vl, "\n");
20  vout.general(m_vl, "%s: being setup (without link smearing).\n",
21  class_name.c_str());
22 
23  typedef AFopr<AFIELD> AltFopr;
24  typedef ASolver<AFIELD> AltSolver;
25 
26  string fopr_type = params_fopr.get_string("fermion_type");
27  fopr_type += "_eo";
28 
29  m_fopr = AltFopr::New(fopr_type, params_fopr);
30 
31  m_dr_smear = 0;
32 
33  m_kernel = 0;
34 
35  string solver_type = params_solver.get_string("solver_type");
36  m_solver = AltSolver::New(solver_type, m_fopr);
37  m_solver->set_parameters(params_solver);
38 
39  reset_performance();
40 
41  vout.general(m_vl, "%s: setup finished.\n", class_name.c_str());
42 }
43 
44 
45 //====================================================================
46 template<typename AFIELD>
47 void Fprop_alt_Standard_eo<AFIELD>::init(const Parameters& params_fopr,
48  const Parameters& params_solver,
49  Director_Smear *dr_smear)
50 {
51  // this constructor assumes that the factories are available.
52  // vout.general(m_vl, "\n");
53  vout.general(m_vl, "%s: being setup (with link smearing).\n",
54  class_name.c_str());
55 
56  typedef AFopr<AFIELD> AltFopr;
57  typedef ASolver<AFIELD> AltSolver;
58 
59  m_dr_smear = dr_smear;
60 
61  string fopr_type = params_fopr.get_string("fermion_type");
62  fopr_type += "_eo";
63  m_kernel = AltFopr::New(fopr_type, params_fopr);
64 
65  // m_fopr = AltFopr::New("Smeared", m_kernel, m_dr_smear);
66  m_fopr = new AFopr_Smeared<AFIELD>(m_kernel, m_dr_smear);
67 
68  string solver_type = params_solver.get_string("solver_type");
69  m_solver = AltSolver::New(solver_type, m_fopr);
70  m_solver->set_parameters(params_solver);
71 
72  reset_performance();
73 
74  vout.general(m_vl, "%s: setup finished.\n", class_name.c_str());
75 }
76 
77 
78 //====================================================================
79 template<typename AFIELD>
80 void Fprop_alt_Standard_eo<AFIELD>::tidyup()
81 {
82  delete m_solver;
83  delete m_fopr;
84  if (m_kernel != 0) delete m_kernel;
85 }
86 
87 
88 //====================================================================
89 template<typename AFIELD>
90 void Fprop_alt_Standard_eo<AFIELD>::set_config(Field *U)
91 {
92  m_fopr->set_config(U);
93 }
94 
95 
96 //====================================================================
97 template<typename AFIELD>
98 void Fprop_alt_Standard_eo<AFIELD>::invert(Field& xq, const Field& b,
99  int& nconv, double& diff)
100 {
101  vout.paranoiac(m_vl, "%s: invert is called.\n", class_name.c_str());
102  vout.paranoiac(m_vl, "mode = %s.\n", m_mode.c_str());
103 
104  if ((m_mode == "D") || (m_mode == "D_prec")) {
105  invert_D(xq, b, nconv, diff);
106  } else if ((m_mode == "DdagD") || (m_mode == "DdagD_prec")) {
107  invert_DdagD(xq, b, nconv, diff);
108  } else {
109  vout.crucial(m_vl, "%s: unsupported mode: %s\n",
110  class_name.c_str(), m_mode.c_str());
111  exit(EXIT_FAILURE);
112  }
113 }
114 
115 
116 //====================================================================
117 template<typename AFIELD>
118 void Fprop_alt_Standard_eo<AFIELD>::invert_D(Field& xq, const Field& b,
119  int& nconv, double& diff)
120 {
121  ThreadManager::assert_single_thread(class_name);
122  vout.paranoiac("invert_D start.\n");
123 
124  m_timer.reset();
125  m_timer.start();
126 
127  // xq.set(0.0);
128 
129  int nin = m_fopr->field_nin();
130  int nvol2 = m_fopr->field_nvol();
131  int nex = m_fopr->field_nex();
132  int nvol = 2 * nvol2;
133 
134  AFIELD axq(nin, nvol, nex);
135  AFIELD abq(nin, nvol, nex);
136 
137  AFIELD be(nin, nvol2, nex), bo(nin, nvol2, nex);
138  AFIELD xe(nin, nvol2, nex), xo(nin, nvol2, nex);
139 
140  AIndex_lex<real_t, AFIELD::IMPL> index_alt;
141  AIndex_eo<real_t, AFIELD::IMPL> index_eo;
142 
143 #pragma omp parallel
144  {
145  if (m_fopr->needs_convert()) {
146  vout.detailed(m_vl, "convert required.\n");
147  m_fopr->convert(abq, b);
148  } else {
149  vout.detailed(m_vl, "convert not required.\n");
150  convert(index_alt, abq, b);
151  }
152 #pragma omp barrier
153  index_eo.split(be, bo, abq);
154  }
155 
156  invert_De(xe, xo, be, bo, nconv, diff);
157 
158 #pragma omp parallel
159  {
160  index_eo.merge(axq, xe, xo);
161 #pragma omp barrier
162 
163  if (m_fopr->needs_convert()) {
164  m_fopr->reverse(xq, axq);
165  } else {
166  reverse(index_alt, xq, axq);
167  }
168  }
169 
170  m_timer.stop();
171  m_elapsed_time += m_timer.elapsed_sec();
172  m_flop_count += m_solver->flop_count();
173 }
174 
175 
176 //====================================================================
177 template<typename AFIELD>
178 void Fprop_alt_Standard_eo<AFIELD>::invert_DdagD(Field& xq, const Field& b,
179  int& nconv, double& diff)
180 {
181  ThreadManager::assert_single_thread(class_name);
182  vout.paranoiac("invert_DdagD start.\n");
183 
184  m_timer.reset();
185  m_timer.start();
186 
187  // xq.set(0.0);
188 
189  int nin = m_fopr->field_nin();
190  int nvol2 = m_fopr->field_nvol();
191  int nex = m_fopr->field_nex();
192  int nvol = 2 * nvol2;
193 
194  AFIELD axq(nin, nvol, nex);
195  AFIELD abq(nin, nvol, nex);
196 
197  AFIELD be(nin, nvol2, nex), bo(nin, nvol2, nex);
198  AFIELD xe(nin, nvol2, nex), xo(nin, nvol2, nex);
199  AFIELD y1(nin, nvol2, nex), y2(nin, nvol2, nex);
200 
201  AIndex_lex<real_t, AFIELD::IMPL> index_alt;
202  AIndex_eo<real_t, AFIELD::IMPL> index_eo;
203 
204 #pragma omp parallel
205  {
206  if (m_fopr->needs_convert()) {
207  vout.detailed(m_vl, "convert required.\n");
208  m_fopr->convert(abq, b);
209  } else {
210  vout.detailed(m_vl, "convert not required.\n");
211  convert(index_alt, abq, b);
212  }
213 #pragma omp barrier
214 
215  index_eo.split(be, bo, abq);
216 #pragma omp barrier
217 
218  m_fopr->mult_gm5(y1, be);
219  m_fopr->mult_gm5(y2, bo);
220  }
221 
222  int nconv1;
223  double diff1;
224  invert_De(xe, xo, y1, y2, nconv1, diff1);
225 
226  nconv = nconv1;
227  diff = diff1;
228 
229 #pragma omp parallel
230  {
231  m_fopr->mult_gm5(y1, xe);
232  m_fopr->mult_gm5(y2, xo);
233  }
234 
235  invert_De(xe, xo, y1, y2, nconv1, diff1);
236 
237  nconv += nconv1;
238  diff += diff1;
239 
240 #pragma omp parallel
241  {
242  index_eo.merge(axq, xe, xo);
243 #pragma omp barrier
244 
245  if (m_fopr->needs_convert()) {
246  m_fopr->reverse(xq, axq);
247  } else {
248  reverse(index_alt, xq, axq);
249  }
250  }
251 }
252 
253 
254 //====================================================================
255 template<typename AFIELD>
256 void Fprop_alt_Standard_eo<AFIELD>::invert(AFIELD& xq, const AFIELD& b,
257  int& nconv, double& diff)
258 {
259  vout.paranoiac(m_vl, "%s: invert is called.\n", class_name.c_str());
260  vout.paranoiac(m_vl, "mode = %s.\n", m_mode.c_str());
261 
262  if (m_mode == "D") {
263  invert_D(xq, b, nconv, diff);
264  } else if (m_mode == "DdagD") {
265  invert_DdagD(xq, b, nconv, diff);
266  } else {
267  vout.crucial(m_vl, "%s: unsupported mode: %s\n",
268  class_name.c_str(), m_mode.c_str());
269  exit(EXIT_FAILURE);
270  }
271 }
272 
273 
274 //====================================================================
275 template<typename AFIELD>
276 void Fprop_alt_Standard_eo<AFIELD>::invert_D(AFIELD& xq, const AFIELD& b,
277  int& nconv, double& diff)
278 {
279  ThreadManager::assert_single_thread(class_name);
280  vout.paranoiac("invert_D start.\n");
281 
282  m_timer.reset();
283  m_timer.start();
284 
285  // xq.set(0.0);
286 
287  int nin = m_fopr->field_nin();
288  int nvol2 = m_fopr->field_nvol();
289  int nex = m_fopr->field_nex();
290 
291  AFIELD be(nin, nvol2, nex), bo(nin, nvol2, nex);
292  AFIELD xe(nin, nvol2, nex), xo(nin, nvol2, nex);
293 
294  AIndex_eo<real_t, AFIELD::IMPL> index_eo;
295 
296 #pragma omp parallel
297  {
298  index_eo.split(be, bo, b);
299  }
300 
301  invert_De(xe, xo, be, bo, nconv, diff);
302 
303 #pragma omp parallel
304  {
305  index_eo.merge(xq, xe, xo);
306  }
307 
308  m_timer.stop();
309  m_elapsed_time += m_timer.elapsed_sec();
310  m_flop_count += m_solver->flop_count();
311 }
312 
313 
314 //====================================================================
315 template<typename AFIELD>
316 void Fprop_alt_Standard_eo<AFIELD>::invert_DdagD(AFIELD& xq, const AFIELD& b,
317  int& nconv, double& diff)
318 {
319  ThreadManager::assert_single_thread(class_name);
320  vout.paranoiac("invert_DdagD start.\n");
321 
322  m_timer.reset();
323  m_timer.start();
324 
325  // xq.set(0.0);
326 
327  int nin = m_fopr->field_nin();
328  int nvol2 = m_fopr->field_nvol();
329  int nex = m_fopr->field_nex();
330 
331  AFIELD be(nin, nvol2, nex), bo(nin, nvol2, nex);
332  AFIELD xe(nin, nvol2, nex), xo(nin, nvol2, nex);
333  AFIELD y1(nin, nvol2, nex), y2(nin, nvol2, nex);
334 
335  AIndex_eo<real_t, AFIELD::IMPL> index_eo;
336 
337 #pragma omp parallel
338  {
339  index_eo.split(be, bo, b);
340 #pragma omp barrier
341 
342  m_fopr->mult_gm5(y1, be);
343  m_fopr->mult_gm5(y2, bo);
344  }
345 
346  int nconv1;
347  double diff1;
348  invert_De(xe, xo, y1, y2, nconv1, diff1);
349 
350  nconv = nconv1;
351  diff = diff1;
352 
353 #pragma omp parallel
354  {
355  m_fopr->mult_gm5(y1, xe);
356  m_fopr->mult_gm5(y2, xo);
357  }
358 
359  invert_De(xe, xo, y1, y2, nconv1, diff1);
360 
361  nconv += nconv1;
362  diff += diff1;
363 
364 #pragma omp parallel
365  {
366  index_eo.merge(xq, xe, xo);
367  }
368 
369  m_timer.stop();
370  m_elapsed_time += m_timer.elapsed_sec();
371  m_flop_count += m_solver->flop_count();
372 }
373 
374 
375 //====================================================================
376 template<typename AFIELD>
377 void Fprop_alt_Standard_eo<AFIELD>::invert_De(AFIELD& xe, AFIELD& xo,
378  AFIELD& be, AFIELD& bo,
379  int& nconv, double& diff)
380 {
381  int nin = m_fopr->field_nin();
382  int nvol2 = m_fopr->field_nvol();
383  int nex = m_fopr->field_nex();
384  int nvol = 2 * nvol2;
385 
386  AFIELD y1(nin, nvol2, nex), y2(nin, nvol2, nex);
387 
388 #pragma omp parallel
389  {
390  // set even source vector.
391  m_fopr->mult(y1, bo, "Doo_inv");
392 
393 #pragma omp barrier
394  m_fopr->mult(y2, y1, "Deo");
395 
396 #pragma omp barrier
397  axpy(be, real_t(-1.0), y2);
398 
399 #pragma omp barrier
400  m_fopr->mult(y1, be, "Dee_inv");
401  }
402 
403  m_fopr->set_mode("D");
404 
405  real_t diff2;
406 
407 #pragma omp parallel
408  {
409  m_solver->solve(xe, y1, nconv, diff2);
410 #pragma omp barrier
411 
412  m_fopr->normalize_fprop(xe);
413 
414  m_fopr->mult(y1, xe, "Doe");
415 #pragma omp barrier
416 
417  aypx(real_t(-1.0), y1, bo);
418 #pragma omp barrier
419 
420  m_fopr->mult(xo, y1, "Doo_inv");
421 #pragma omp barrier
422  }
423 
424  diff = double(diff2);
425 }
426 
427 
428 //====================================================================
429 template<typename AFIELD>
430 double Fprop_alt_Standard_eo<AFIELD>::flop_count()
431 {
432  return m_solver->flop_count();
433 }
434 
435 
436 //====================================================================
437 template<typename AFIELD>
438 void Fprop_alt_Standard_eo<AFIELD>::reset_performance()
439 {
440  m_flop_count = 0.0;
441  m_elapsed_time = 0.0;
442 }
443 
444 
445 //====================================================================
446 template<typename AFIELD>
447 void Fprop_alt_Standard_eo<AFIELD>::report_performance()
448 {
449  double flops = m_flop_count / m_elapsed_time;
450  double gflops = flops * 1.0e-9;
451 
452  vout.general(m_vl, "\n");
453  vout.general(m_vl, "%s: solver performance:\n", class_name.c_str());
454  vout.general(m_vl, " Elapsed time = %14.6f sec\n", m_elapsed_time);
455  vout.general(m_vl, " Flop(total) = %18.0f\n", m_flop_count);
456  vout.general(m_vl, " Performance = %11.3f GFlops\n", gflops);
457 }
458 
459 
460 //====================================================================
461 template<typename AFIELD>
462 void Fprop_alt_Standard_eo<AFIELD>::mult_performance(
463  const std::string mode,
464  const int Nrepeat)
465 {
466  int nin = m_fopr->field_nin();
467  int nvol = m_fopr->field_nvol();
468  int nex = m_fopr->field_nex();
469 
470  AFIELD axq(nin, nvol, nex), abq(nin, nvol, nex);
471  abq.set(0.0);
472  abq.set(0, 1.0);
473 
474  unique_ptr<Timer> timer(new Timer);
475 
476  std::string mode_prev = m_fopr->get_mode();
477  m_fopr->set_mode(mode);
478 
479  timer->start();
480 
481 #pragma omp parallel
482  {
483  for (int i = 0; i < Nrepeat; ++i) {
484  m_fopr->mult(axq, abq);
485  m_fopr->mult(abq, axq);
486  }
487  }
488 
489  timer->stop();
490 
491  double flop_fopr = m_fopr->flop_count();
492  double flop_total = flop_fopr * double(2 * Nrepeat);
493 
494  double elapsed_time = timer->elapsed_sec();
495  double flops = flop_total / elapsed_time;
496  double gflops = flops * 1.0e-9;
497 
498  vout.general(m_vl, "\n");
499  vout.general(m_vl, "%s: mult performance:\n", class_name.c_str());
500  vout.general(m_vl, " mult mode = %s\n", mode.c_str());
501  vout.general(m_vl, " Elapsed time = %14.6f sec\n", elapsed_time);
502  vout.general(m_vl, " Flop(Fopr) = %18.0f\n", flop_fopr);
503  vout.general(m_vl, " Flop(total) = %18.0f\n", flop_total);
504  vout.general(m_vl, " Performance = %11.3f GFlops\n", gflops);
505 
506  m_fopr->set_mode(mode_prev);
507 }
508 
509 
510 //====================================================================
511 //============================================================END=====
Fprop_alt_Standard_eo::report_performance
void report_performance()
Definition: fprop_alt_Standard_eo-tmpl.h:447
AFopr
Definition: afopr.h:48
ASolver
Definition: asolver.h:23
Fprop_alt_Standard_eo::invert
void invert(Field &, const Field &, int &, double &)
invert accordingly to the mode. [22 Sep 2018 H.Matsufuru]
Definition: fprop_alt_Standard_eo-tmpl.h:98
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
AIndex_lex
Definition: aindex_lex_base.h:17
convert
void convert(INDEX &index, AFIELD &v, const Field &w)
Definition: afield-inc.h:129
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
Fprop_alt_Standard_eo::init
void init(const Parameters &params_fopr, const Parameters &params_solver)
Definition: fprop_alt_Standard_eo-tmpl.h:15
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
Fprop_alt_Standard_eo::mult_performance
void mult_performance(const std::string mode, const int Nrepeat)
Definition: fprop_alt_Standard_eo-tmpl.h:462
Fprop_alt_Standard_eo::set_config
void set_config(Field *)
Definition: fprop_alt_Standard_eo-tmpl.h:90
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
Fprop_alt_Standard_eo::flop_count
double flop_count()
Definition: fprop_alt_Standard_eo-tmpl.h:430
Timer
Definition: timer.h:31
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:238
Fprop_alt_Standard_eo
Get quark propagator for Fopr with lexical site index: alternative version.
Definition: fprop_alt_Standard_eo.h:30
Fprop_alt_Standard_eo::invert_De
void invert_De(AFIELD &, AFIELD &, AFIELD &, AFIELD &, int &, double &)
Definition: fprop_alt_Standard_eo-tmpl.h:377
Fprop_alt_Standard_eo::reset_performance
void reset_performance()
Definition: fprop_alt_Standard_eo-tmpl.h:438
AFopr_Smeared
smeared fermion operator: alternative version.
Definition: afopr_Smeared.h:41
Fprop_alt::real_t
AFIELD::real_t real_t
Definition: fprop_alt.h:32
Fprop_alt_Standard_eo::invert_DdagD
void invert_DdagD(Field &, const Field &, int &, double &)
Definition: fprop_alt_Standard_eo-tmpl.h:178
Fprop_alt_Standard_eo::invert_D
void invert_D(Field &, const Field &, int &, double &)
Definition: fprop_alt_Standard_eo-tmpl.h:118
Director_Smear
Manager of smeared configurations.
Definition: director_Smear.h:39
Fprop_alt_Standard_eo::tidyup
void tidyup()
Definition: fprop_alt_Standard_eo-tmpl.h:80
reverse
void reverse(INDEX &index, Field &v, const AFIELD &w)
Definition: afield-inc.h:159
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
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
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
AIndex_eo< real_t, AFIELD::IMPL >