Bridge++  Ver.2.1.3
afopr_Domainwall_5din_eo-tmpl.h
Go to the documentation of this file.
1 
11 
13 
14 template<typename AFIELD>
16  = "AFopr_Domainwall_5din_eo";
17 
18 #define AFOPR_TIMER
19 
20 #ifdef AFOPR_TIMER
21 #include "lib/Tools/timer.h"
22 #define START_TIMER(var_timer) var_timer->start()
23 #define STOP_TIMER(var_timer) var_timer->stop()
24 #else
25 #define START_TIMER(var_timer)
26 #define STOP_TIMER(var_timer)
27 #endif
28 
29 //====================================================================
30 template<typename AFIELD>
32 {
34 
35  std::string vlevel;
36  if (!params.fetch_string("verbose_level", vlevel)) {
37  m_vl = vout.set_verbose_level(vlevel);
38  } else {
39  m_vl = CommonParameters::Vlevel();
40  }
41 
42  vout.general(m_vl, "%s: construction\n", class_name.c_str());
44 
45  m_repr = "Dirac"; // now only the Dirac repr is available.
46 
47  std::string repr;
48  if (!params.fetch_string("gamma_matrix_type", repr)) {
49  if (repr != "Dirac") {
50  vout.crucial("Error at %s: unsupported gamma-matrix type: %s\n",
51  class_name.c_str(), repr.c_str());
52  exit(EXIT_FAILURE);
53  }
54  }
55 
56  int Nc = CommonParameters::Nc();
57  if (Nc != 3) {
58  vout.crucial("%s: only applicable to Nc = 3\n",
59  class_name.c_str());
60  exit(EXIT_FAILURE);
61  }
62 
63  // reset the timers
64 #ifdef AFOPR_TIMER
65  timer_mult_Deo.reset( new Timer("afopr_Domainwall_5din_eo: Deo & dagger "));
66  timer_mult_Dee_inv.reset( new Timer("afopr_Domainwall_5din_eo: Dee_inv & dag "));
67  timer_pack.reset( new Timer("afopr_Domainwall_5din_eo: pack "));
68  timer_bulk.reset( new Timer("afopr_Domainwall_5din_eo: bulk "));
69  timer_boundary.reset( new Timer("afopr_Domainwall_5din_eo: boundary "));
70  timer_comm.reset( new Timer("afopr_Domainwall_5din_eo: comm "));
71  timer_comm_recv_wait.reset( new Timer("afopr_Domainwall_5din_eo: comm_recv_wait "));
72  timer_comm_send_wait.reset( new Timer("afopr_Domainwall_5din_eo: comm_send_wait "));
73  timer_comm_recv_start.reset(new Timer("afopr_Domainwall_5din_eo: comm_recv_start"));
74  timer_comm_send_start.reset(new Timer("afopr_Domainwall_5din_eo: comm_send_start"));
75 #pragma omp barrier
76  vout.detailed(m_vl, "%s: detailed timer was initialized\n", class_name.c_str());
77 #endif
78 
79  m_Nd = CommonParameters::Nd();
80  m_Nd2 = m_Nd/2;
81  m_Nvcd = 2 * Nc * m_Nd;
82  m_Ndf = 2 * Nc * Nc;
83 
84  m_Nx = CommonParameters::Nx();
85  m_Ny = CommonParameters::Ny();
86  m_Nz = CommonParameters::Nz();
87  m_Nt = CommonParameters::Nt();
88  m_Nvol = CommonParameters::Nvol();
89  m_Ndim = CommonParameters::Ndim();
90 
91  m_Nx2 = m_Nx / 2;
92  m_Nst2 = m_Nvol / 2;
93 
94  int ipe3 = Communicator::ipe(3);
95  int ipe2 = Communicator::ipe(2);
96  int ipe1 = Communicator::ipe(1);
97  m_Ieo_origin = (ipe1 * m_Ny + ipe2 * m_Nz + ipe3 * m_Nt) % 2;
98 
99  // condition check
100  if (m_Nx % 2 != 0) {
101  vout.crucial(m_vl, "%s: Nx must be even.\n",
102  class_name.c_str());
103  exit(EXIT_FAILURE);
104  }
105 
106  m_Nsize[0] = m_Nx2;
107  m_Nsize[1] = m_Ny;
108  m_Nsize[2] = m_Nz;
109  m_Nsize[3] = m_Nt;
110 
111  // switches for coomunication
112  int req_comm = 0; // 0: communication only if necessary
113  // 1: communication is enforced any time
114  if (!params.fetch_int("require_communication", req_comm)) {
115  vout.general(m_vl, "req_comm = %d (input)\n", req_comm);
116  } else {
117  vout.general(m_vl, "req_comm = %d (default)\n", req_comm);
118  }
119 
120  do_comm_any = 0;
121  for (int mu = 0; mu < m_Ndim; ++mu) {
122  do_comm[mu] = 1;
123  if ((req_comm == 0) && (Communicator::npe(mu) == 1)) do_comm[mu] = 0;
124  do_comm_any += do_comm[mu];
125  vout.general("do_comm[%d] = %d\n", mu, do_comm[mu]);
126  }
127 
128  m_Ns = 0; // temporary set
129  set_parameters(params);
130 
131  m_Nbdsize.resize(m_Ndim);
132  int Nbdin = (m_Nvcd / 2) * m_Ns;
133  m_Nbdsize[0] = Nbdin * ceil_nwp((m_Ny * m_Nz * m_Nt + 1)/2);
134  m_Nbdsize[1] = Nbdin * ceil_nwp(m_Nx2 * m_Nz * m_Nt);
135  m_Nbdsize[2] = Nbdin * ceil_nwp(m_Nx2 * m_Ny * m_Nt);
136  m_Nbdsize[3] = Nbdin * ceil_nwp(m_Nx2 * m_Ny * m_Nz);
137 
138  setup_channels();
139 
140  // gauge configuration.
141  int Nst_pad2 = 2 * ceil_nwp(m_Nst2);
142  m_Ueo.reset(m_Ndf, Nst_pad2, m_Ndim);
143 
145  vout.detailed(m_vl, "%s: initalization finished.\n",
146  class_name.c_str());
147 }
148 
149 //====================================================================
150 template<typename AFIELD>
152 {
154 
155  // openacc device memory clean up
156  for(int mu = 0; mu < m_Ndim; ++mu){
157 
158 #ifdef USE_MPI
159  real_t* buf_dn1 = (real_t*)chsend_dn[mu].ptr();
160  real_t* buf_dn2 = (real_t*)chrecv_dn[mu].ptr();
161  real_t* buf_up1 = (real_t*)chrecv_up[mu].ptr();
162  real_t* buf_up2 = (real_t*)chsend_up[mu].ptr();
163  BridgeACC::afield_tidyup(buf_dn1, m_Nbdsize[mu]);
164  BridgeACC::afield_tidyup(buf_dn2, m_Nbdsize[mu]);
165  BridgeACC::afield_tidyup(buf_up1, m_Nbdsize[mu]);
166  BridgeACC::afield_tidyup(buf_up2, m_Nbdsize[mu]);
167 #else
168  real_t* buf_up = (real_t*)chsend_up[mu].ptr();
169  real_t* buf_dn = (real_t*)chsend_dn[mu].ptr();
170  BridgeACC::afield_tidyup(buf_up, m_Nbdsize[mu]);
171  BridgeACC::afield_tidyup(buf_dn, m_Nbdsize[mu]);
172 #endif
173 
174  }
175 
176 #ifdef AFOPR_TIMER
177  timer_mult_Deo->report();
178  timer_mult_Dee_inv->report();
179  timer_pack->report();
180  timer_bulk->report();
181  timer_boundary->report();
182  timer_comm->report();
183  timer_comm_recv_wait->report();
184  timer_comm_send_wait->report();
185  timer_comm_recv_start->report();
186  timer_comm_send_start->report();
187 #endif
188 
189 }
190 
191 
192 //====================================================================
193 template<typename AFIELD>
195 {
197 
198  chsend_up.resize(m_Ndim);
199  chrecv_up.resize(m_Ndim);
200  chsend_dn.resize(m_Ndim);
201  chrecv_dn.resize(m_Ndim);
202 
203  for (int mu = 0; mu < m_Ndim; ++mu) {
204  size_t Nvsize = m_Nbdsize[mu] * sizeof(real_t);
205 
206  chsend_dn[mu].send_init(Nvsize, mu, -1);
207  chsend_up[mu].send_init(Nvsize, mu, 1);
208 #ifdef USE_MPI
209  chrecv_up[mu].recv_init(Nvsize, mu, 1);
210  chrecv_dn[mu].recv_init(Nvsize, mu, -1);
211 #else
212  void *buf_up = (void *)chsend_dn[mu].ptr();
213  chrecv_up[mu].recv_init(Nvsize, mu, 1, buf_up);
214  void *buf_dn = (void *)chsend_up[mu].ptr();
215  chrecv_dn[mu].recv_init(Nvsize, mu, -1, buf_dn);
216 #endif
217 
218  if (do_comm[mu] == 1) {
219  chset_send.append(chsend_up[mu]);
220  chset_send.append(chsend_dn[mu]);
221  chset_recv.append(chrecv_up[mu]);
222  chset_recv.append(chrecv_dn[mu]);
223  }
224 
225  // openacc device memory allocation
226 #ifdef USE_MPI
227  real_t* buf_dn1 = (real_t*)chsend_dn[mu].ptr();
228  real_t* buf_dn2 = (real_t*)chrecv_dn[mu].ptr();
229  real_t* buf_up1 = (real_t*)chsend_up[mu].ptr();
230  real_t* buf_up2 = (real_t*)chrecv_up[mu].ptr();
231  BridgeACC::afield_init(buf_dn1, m_Nbdsize[mu]);
232  BridgeACC::afield_init(buf_dn2, m_Nbdsize[mu]);
233  BridgeACC::afield_init(buf_up1, m_Nbdsize[mu]);
234  BridgeACC::afield_init(buf_up2, m_Nbdsize[mu]);
235 #else
236  BridgeACC::afield_init((real_t*)buf_up, m_Nbdsize[mu]);
237  BridgeACC::afield_init((real_t*)buf_dn, m_Nbdsize[mu]);
238 #endif
239 
240  }
241 
242 }
243 
244 //====================================================================
245 template<typename AFIELD>
247  const Parameters& params)
248 {
249  std::string vlevel;
250  if (!params.fetch_string("verbose_level", vlevel)) {
251  m_vl = vout.set_verbose_level(vlevel);
252  }
253 
254  string gmset_type;
255  double mq, M0;
256  int Ns;
257  std::vector<int> bc;
258  double b, c, alpha;
259 
260  int err_optional = 0;
261  err_optional += params.fetch_string("gamma_matrix_type", m_repr);
262  if(err_optional){
263  vout.crucial(m_vl, " gamma_matrix_type is not given\n");
264  }
265 
266  err_optional = 0;
267  err_optional += params.fetch_string("code_implementation", m_impl);
268  if(err_optional){
269  vout.crucial(m_vl, " code_implementation is not given\n");
270  m_impl = "5d";
271  }
272 
273  int err = 0;
274  err += params.fetch_double("quark_mass", mq);
275  err += params.fetch_double("domain_wall_height", M0);
276  err += params.fetch_int("extent_of_5th_dimension", Ns);
277  err += params.fetch_int_vector("boundary_condition", bc);
278 
279  if (err) {
280  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
281  class_name.c_str());
282  exit(EXIT_FAILURE);
283  }
284 
285  int err2 = 0;
286  err2 += params.fetch_double("coefficient_b", b);
287  err2 += params.fetch_double("coefficient_c", c);
288 
289  if (err2) {
290  vout.general(m_vl, " coefficients b, c are not provided:"
291  " set to Shamir's form.\n");
292  b = 1.0;
293  c = 0.0;
294  }
295 
296  int err3 = 0;
297  err3 += params.fetch_double("parameter_alpha", alpha);
298  if (err3) {
299  vout.general(m_vl, " parameter alpha is not provided: set to 1.0.\n");
300  alpha = 1.0;
301  }
302 
303  set_parameters(real_t(mq), real_t(M0), Ns, bc,
304  real_t(b), real_t(c), real_t(alpha));
305 
306 }
307 
308 //====================================================================
309 template<typename AFIELD>
311  Parameters& params) const
312 {
313  params.set_string("kernel_type", m_kernel_type);
314  params.set_string("gamma_matrix_type", m_repr);
315  params.set_string("code_implementation", m_impl);
316  params.set_double("quark_mass", double(m_mq));
317  params.set_double("domain_wall_height", double(m_M0));
318  params.set_int("extent_of_5th_dimension", m_Ns);
319  params.set_int_vector("boundary_condition", m_boundary);
320  params.set_double("coefficient_b", double(m_b[0]));
321  params.set_double("coefficient_c", double(m_c[0]));
322  params.set_double("parameter_alpha", double(m_alpha));
323  params.set_string("gamma_matrix_type", m_repr);
324 
325  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
326 }
327 
328 
329 //====================================================================
330 template<typename AFIELD>
332  const real_t mq,
333  const real_t M0,
334  const int Ns,
335  const std::vector<int> bc,
336  const real_t b,
337  const real_t c,
338  const real_t alpha)
339 {
340 #pragma omp barrier
341 
342  int ith = ThreadManager::get_thread_id();
343 
344  if (ith == 0) {
345  m_M0 = real_t(M0);
346  m_mq = real_t(mq);
347  m_Ns = Ns;
348  m_NinF = m_Nvcd * m_Ns;
349  m_alpha = alpha;
350 
351  assert(bc.size() == m_Ndim);
352  if (m_boundary.size() != m_Ndim) m_boundary.resize(m_Ndim);
353 
354  for (int mu = 0; mu < m_Ndim; ++mu) {
355  m_boundary[mu] = bc[mu];
356  m_bc[mu] = 1;
357  if(do_comm[mu] > 0){ // do communication
358  if(Communicator::ipe(mu) == 0) m_bc[mu] = m_boundary[mu];
359  m_bc2[mu] = 0;
360  }else{ // no communication
361  m_bc[mu] = 0; // for boundar part (dummy)
362  m_bc2[mu] = m_boundary[mu]; // for bulk part
363  }
364  }
365 
366  if (m_b.size() != m_Ns) {
367  m_b.resize(m_Ns);
368  m_c.resize(m_Ns);
369  }
370  for (int is = 0; is < m_Ns; ++is) {
371  m_b[is] = real_t(b);
372  m_c[is] = real_t(c);
373  }
374  }
375 
376 #pragma omp barrier
377 
378  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
379  vout.general(m_vl, " gamma matrix repr.: %s\n", m_repr.c_str());
380  vout.general(m_vl, " code implementation: %s\n", m_impl.c_str());
381  vout.general(m_vl, " mq = %8.4f\n", m_mq);
382  vout.general(m_vl, " M0 = %8.4f\n", m_M0);
383  vout.general(m_vl, " Ns = %4d\n", m_Ns);
384  for (int mu = 0; mu < m_Ndim; ++mu) {
385  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
386  }
387  vout.general(m_vl, " coefficients:\n");
388  for (int is = 0; is < m_Ns; ++is) {
389  vout.general(m_vl, " b[%2d] = %16.10f c[%2d] = %16.10f\n",
390  is, m_b[is], is, m_c[is]);
391  }
392  vout.general(m_vl, " alpha = %8.4f\n", m_alpha);
393 
394  // working 5d vectors.
395  if (m_w1.nex() != Ns) {
396  m_w1.reset(m_NinF, m_Nst2, 1);
397  m_v1.reset(m_NinF, m_Nst2, 1);
398  m_v2.reset(m_NinF, m_Nst2, 1);
399  }
400 
401  set_precond_parameters();
402 
403 #pragma omp barrier
404 }
405 
406 
407 //====================================================================
408 template<typename AFIELD>
410  const std::vector<real_t> vec_b,
411  const std::vector<real_t> vec_c)
412 {
413 #pragma omp barrier
414 
415  if ((vec_b.size() != m_Ns) || (vec_c.size() != m_Ns)) {
416  vout.crucial(m_vl, "%s: size of coefficient vectors incorrect.\n",
417  class_name.c_str());
418  }
419 
420  vout.general(m_vl, "%s: coefficient vectors are set:\n",
421  class_name.c_str());
422 
423  int ith = ThreadManager::get_thread_id();
424  if (ith == 0) {
425  for (int is = 0; is < m_Ns; ++is) {
426  m_b[is] = vec_b[is];
427  m_c[is] = vec_c[is];
428  vout.general(m_vl, "b[%2d] = %16.10f c[%2d] = %16.10f\n",
429  is, m_b[is], is, m_c[is]);
430  }
431  }
432 
433  set_precond_parameters();
434 
435 #pragma omp barrier
436 }
437 
438 //====================================================================
439 template<typename AFIELD>
441 {
442  int ith = ThreadManager::get_thread_id();
443  if (ith == 0) {
444 
445  if (m_dp.size() != m_Ns) {
446  m_dp.resize(m_Ns);
447  m_dm.resize(m_Ns);
448  m_dpinv.resize(m_Ns);
449  m_e.resize(m_Ns - 1);
450  m_f.resize(m_Ns - 1);
451  }
452 
453  for (int is = 0; is < m_Ns; ++is) {
454  m_dp[is] = m_alpha * (1.0 + m_b[is] * (4.0 - m_M0));
455  m_dm[is] = m_alpha * (1.0 - m_c[is] * (4.0 - m_M0));
456  }
457 
458  m_e[0] = m_mq * m_dm[m_Ns - 1] / m_dp[0];
459  m_f[0] = m_mq * m_dm[0]/m_alpha;
460  for (int is = 1; is < m_Ns - 1; ++is) {
461  m_e[is] = m_e[is - 1] * m_dm[is - 1] / m_dp[is];
462  m_f[is] = m_f[is - 1] * m_dm[is] / m_dp[is - 1];
463  }
464 
465  m_g = m_e[m_Ns - 2] * m_dm[m_Ns - 2];
466 
467  for (int is = 0; is < m_Ns - 1; ++is) {
468  m_dpinv[is] = 1.0 / m_dp[is];
469  }
470  m_dpinv[m_Ns - 1] = 1.0 / (m_dp[m_Ns - 1] + m_g);
471 
472  }
473 
474  // setting inverse matrix
475  set_matrix5d_inverse();
476 
477 }
478 
479 //====================================================================
480 template<typename AFIELD>
482 {
483  int Nc = CommonParameters::Nc();
484 
485  int mat_size = m_Nd2 * m_Ns;
486  m_mat_inv.resize(mat_size * mat_size);
487 
488  for(int isb = 0; isb < m_Ns; ++isb){
489  for(int idb = 0; idb < m_Nd2; ++idb){
490  m_w1.set(0.0);
491  int inb = 2 * Nc * (idb * 2 + m_Nd * isb);
492  int site = 0;
493  int idxb = m_index_eo.idxh(inb, m_NinF, site, 0);
494  m_w1.set(idxb, 1.0);
495 
496  D_ee_inv(m_v2, m_w1, 0);
497 
498  for(int isa = 0; isa < m_Ns; ++isa){
499  for(int ida = 0; ida < m_Nd2; ++ida){
500  int ina = 2 * Nc * (ida * 2 + m_Nd * isa);
501  int idxa = m_index_eo.idxh(ina, m_NinF, site, 0);
502  m_mat_inv[mat_index(ida, isa, idb, isb)] = m_v2.cmp(idxa);
503  }
504  }
505 
506  }
507  }
508 
509 }
510 
511 
512 //====================================================================
513 template<typename AFIELD>
515 {
516  int nth = ThreadManager::get_num_threads();
517 
518  vout.detailed(m_vl, "%s: set_config is called: num_threads = %d\n",
519  class_name.c_str(), nth);
520 
521  if (nth > 1) {
522  set_config_impl(u);
523  } else {
524  set_config_omp(u);
525  }
526 
527  vout.detailed(m_vl, "%s: set_config finished\n", class_name.c_str());
528 }
529 
530 //====================================================================
531 template<typename AFIELD>
533 {
534  vout.detailed(m_vl, " set_config_omp is called.\n");
535 
536 #pragma omp parallel
537  {
538  set_config_impl(u);
539  }
540 }
541 
542 //====================================================================
543 template<typename AFIELD>
545 {
546 #pragma omp barrier
547 
548  convert_gauge(m_index_eo, m_Ueo, *u);
549 
550 #pragma omp barrier
551 }
552 
553 //====================================================================
554 template<typename AFIELD>
556 {
557 #pragma omp barrier
558 
559  int Nst = v.nvol();
560  if(w.nvol() != Nst){
561  vout.crucial(m_vl, "%s: convert: field size irrelevant\n");
562  exit(EXIT_FAILURE);
563  }
564 
565  int ith, nth, isite, nsite;
566  set_threadtask_afopr(ith, nth, isite, nsite, Nst);
567 
569 
570  for (int site = isite; site < nsite; ++site) {
571  for (int is = 0; is < m_Ns; ++is) {
572  for (int ivcd = 0; ivcd < NVCD; ++ivcd) {
573  int in_alt = ivcd + NVCD * is;
574  real_t vt = real_t(w.cmp(ivcd, site, is));
575  v.set_host(index.idx(in_alt, m_NinF, site, 0), vt);
576  }
577  }
578  } // site-llop
579 
580  v.update_device();
581 
582 }
583 
584 
585 //====================================================================
586 template<typename AFIELD>
588 {
589 #pragma omp barrier
590 
591  int Nst = v.nvol();
592  if(w.nvol() != Nst){
593  vout.crucial(m_vl, "%s: convert: field size irrelevant\n");
594  exit(EXIT_FAILURE);
595  }
596 
597  int ith, nth, isite, nsite;
598  set_threadtask_afopr(ith, nth, isite, nsite, Nst);
599 
600  w.update_host();
601 
603 
604  for (int site = isite; site < nsite; ++site) {
605  for (int is = 0; is < m_Ns; ++is) {
606  for (int ivcd = 0; ivcd < NVCD; ++ivcd) {
607  int in_alt = ivcd + NVCD * is;
608  double vt = double(w.cmp_host(index.idx(in_alt, m_NinF, site, 0)));
609  v.set(ivcd, site, is, vt);
610  }
611  }
612  } // site-loop
613 
614 #pragma omp barrier
615 }
616 
617 
618 //====================================================================
619 template<typename AFIELD>
621 {
622 #pragma omp barrier
623 
624  int ith = ThreadManager::get_thread_id();
625  if (ith == 0) m_mode = mode;
626  vout.paranoiac(m_vl, " mode is set to %s\n", mode.c_str());
627 
628 #pragma omp barrier
629 }
630 
631 
632 //====================================================================
633 template<typename AFIELD>
635 {
636  if (m_mode == "D") {
637  D(v, w);
638  //D_alt(v, w);
639  } else if (m_mode == "Ddag") {
640  Ddag(v, w);
641  //Ddag_alt(v, w);
642  } else if (m_mode == "DdagD") {
643  DdagD(v, w);
644  // DdagD_alt(v, w);
645  //} else if (m_mode == "DDdag") { // this mode is not available
646  // Ddag(m_w1, w); // this fails because m_w1 is used in D.
647  // D(v, m_w1);
648  //} else if (m_mode == "H") {
649  // H(v, w);
650  }else if (m_mode == "Deo") {
651  D_eo(v, w, 0);
652  } else if (m_mode == "Doe") {
653  D_eo(v, w, 1);
654  } else if (m_mode == "Dee") {
655  D_ee(v, w, 0);
656  } else if (m_mode == "Doo") {
657  D_ee(v, w, 1);
658  } else if (m_mode == "Dee_inv") {
659  D_ee_inv(v, w, 0);
660  //D_ee_inv_alt(v, w, 0);
661  } else if (m_mode == "Doo_inv") {
662  D_ee_inv(v, w, 1);
663  //D_ee_inv_alt(v, w, 1);
664  } else {
665  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
666  vout.crucial(m_vl, "in mult, mode=%s.\n", m_mode.c_str());
667  exit(EXIT_FAILURE);
668  }
669 }
670 
671 
672 //====================================================================
673 template<typename AFIELD>
675 {
676  if (m_mode == "D") {
677  Ddag(v, w);
678  //Ddag_alt(v, w);
679  } else if (m_mode == "Ddag") {
680  D(v, w);
681  //D_alt(v, w);
682  } else if (m_mode == "DdagD") {
683  DdagD(v, w);
684  //DdagD_alt(v, w);
685  //} else if (m_mode == "H") {
686  //Hdag(v, w);
687  } else if (m_mode == "Deo") {
688  Ddag_eo(v, w, 1);
689  } else if (m_mode == "Doe") {
690  Ddag_eo(v, w, 0);
691  } else if (m_mode == "Dee") {
692  Ddag_ee(v, w, 0);
693  } else if (m_mode == "Doo") {
694  Ddag_ee(v, w, 1);
695  } else if (m_mode == "Dee_inv") {
696  Ddag_ee_inv(v, w, 0);
697  //Ddag_ee_inv_alt(v, w, 0);
698  } else if (m_mode == "Doo_inv") {
699  Ddag_ee_inv(v, w, 1);
700  //Ddag_ee_inv_alt(v, w, 1);
701  } else {
702  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
703  vout.crucial(m_vl, "in mult_dag, mode=%s.\n", m_mode.c_str());
704  exit(EXIT_FAILURE);
705  }
706 }
707 
708 
709 //====================================================================
710 template<typename AFIELD>
712  std::string mode)
713 {
714  assert(w.check_size(m_NinF, m_Nst2, 1));
715  assert(v.check_size(m_NinF, m_Nst2, 1));
716 
717  if (mode == "Deo") {
718  D_eo(v, w, 0);
719  } else if (mode == "Doe") {
720  D_eo(v, w, 1);
721  } else if (mode == "Dee") {
722  D_ee(v, w, 0);
723  } else if (mode == "Doo") {
724  D_ee(v, w, 1);
725  } else if (mode == "Dee_inv") {
726  D_ee_inv(v, w, 0);
727  //D_ee_inv_alt(v, w, 0);
728  } else if (mode == "Doo_inv") {
729  D_ee_inv(v, w, 1);
730  //D_ee_inv_alt(v, w, 1);
731  } else {
732  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
733  vout.crucial(m_vl, "in mult, mode=%s.\n", m_mode.c_str());
734  exit(EXIT_FAILURE);
735  }
736 
737 }
738 
739 //====================================================================
740 template<typename AFIELD>
742  const AFIELD& w,
743  std::string mode)
744 {
745  assert(w.check_size(m_NinF, m_Nst2, 1));
746  assert(v.check_size(m_NinF, m_Nst2, 1));
747 
748  if (mode == "Deo") {
749  Ddag_eo(v, w, 1);
750  } else if (mode == "Doe") {
751  Ddag_eo(v, w, 0);
752  } else if (mode == "Dee") {
753  Ddag_ee(v, w, 0);
754  } else if (mode == "Doo") {
755  Ddag_ee(v, w, 1);
756  } else if (mode == "Dee_inv") {
757  Ddag_ee_inv(v, w, 0);
758  //Ddag_ee_inv_alt(v, w, 0);
759  } else if (mode == "Doo_inv") {
760  Ddag_ee_inv(v, w, 1);
761  //Ddag_ee_inv_alt(v, w, 1);
762  } else {
763  vout.crucial(m_vl, "mode undeifined in %s.\n", class_name.c_str());
764  vout.crucial(m_vl, "in mult, mode=%s.\n", m_mode.c_str());
765  exit(EXIT_FAILURE);
766  }
767 }
768 
769 //====================================================================
770 template<typename AFIELD>
772 {
773  D_eo(m_v1, w, 1);
774  LU_inv(m_v2, m_v1);
775  D_eo(m_v1, m_v2, 0);
776  LU_inv(m_v2, m_v1);
777 
778  copy(m_v1, w);
779  axpy(m_v1, real_t(-1.0), m_v2);
780 
781  LUdag_inv(v, m_v1);
782  Ddag_eo(m_v2, v, 1);
783  LUdag_inv(v, m_v2);
784  Ddag_eo(m_v2, v, 0);
785 
786  copy(v, m_v1);
787  axpy(v, real_t(-1.0), m_v2);
788 }
789 
790 
791 //====================================================================
792 template<typename AFIELD>
794 {
795  D_eo(m_v1, w, 1);
796  LU_inv(m_v2, m_v1);
797  D_eo(m_v1, m_v2, 0);
798  LU_inv(m_v2, m_v1);
799 
800  copy(v, w);
801  axpy(v, real_t(-1.0), m_v2);
802 }
803 
804 //====================================================================
805 template<typename AFIELD>
807 {
808  assert(w.check_size(m_NinF, m_Nst2, 1));
809  assert(v.check_size(m_NinF, m_Nst2, 1));
810 
811 #pragma omp barrier
812 
813  LUdag_inv(m_v1, w);
814  Ddag_eo(m_v2, m_v1, 1);
815  LUdag_inv(m_v1, m_v2);
816  Ddag_eo(m_v2, m_v1, 0);
817 
818  copy(v, w);
819  axpy(v, real_t(-1.0), m_v2);
820 
821 }
822 
823 
824 //====================================================================
825 template<typename AFIELD>
827 {
828  D_eo(m_v1, w, 1);
829  D_ee_inv_alt(m_v2, m_v1, 1);
830  D_eo(m_v1, m_v2, 0);
831  D_ee_inv_alt(m_v2, m_v1, 0);
832 
833  copy(m_v1, w);
834  axpy(m_v1, real_t(-1.0), m_v2);
835 
836  Ddag_ee_inv_alt(v, m_v1, 0);
837  Ddag_eo(m_v2, v, 1);
838  Ddag_ee_inv_alt(v, m_v2, 1);
839  Ddag_eo(m_v2, v, 0);
840 
841  copy(v, m_v1);
842  axpy(v, real_t(-1.0), m_v2);
843 }
844 
845 //====================================================================
846 template<typename AFIELD>
848 {
849  D_eo(m_v1, w, 1);
850  D_ee_inv_alt(m_v2, m_v1, 1);
851  D_eo(m_v1, m_v2, 0);
852  D_ee_inv_alt(m_v2, m_v1, 0);
853 
854  copy(v, w);
855  axpy(v, real_t(-1.0), m_v2);
856 }
857 
858 
859 //====================================================================
860 template<typename AFIELD>
862 {
863  assert(w.check_size(m_NinF, m_Nst2, 1));
864  assert(v.check_size(m_NinF, m_Nst2, 1));
865 
866 #pragma omp barrier
867 
868  Ddag_ee_inv_alt(m_v1, w, 0);
869  Ddag_eo(m_v2, m_v1, 1);
870  Ddag_ee_inv_alt(m_v1, m_v2, 1);
871  Ddag_eo(m_v2, m_v1, 0);
872 
873  copy(v, w);
874  axpy(v, real_t(-1.0), m_v2);
875 
876 }
877 
878 
879 //====================================================================
880 template<typename AFIELD>
882  const AFIELD& w)
883 {
884  real_t *vp = v.ptr(0);
885  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
886 
887 #pragma omp barrier
888 
889  int ith = ThreadManager::get_thread_id();
890 
891  if (ith == 0)
893  vp, wp, m_Ns, m_Nsize);
894 
895 #pragma omp barrier
896 }
897 
898 //====================================================================
899 template<typename AFIELD>
901  const AFIELD& w)
902 {
903  assert(w.check_size(m_NinF, m_Nst2, 1));
904  assert(v.check_size(m_NinF, m_Nst2, 1));
905 
906 #pragma omp barrier
907 
908  real_t *vp = v.ptr(0);
909  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
910 
911  int ith = ThreadManager::get_thread_id();
912 
913  if (ith == 0)
914  BridgeACC::mult_domainwall_5din_mult_R(vp, wp, m_Ns, m_Nsize);
915 
916 #pragma omp barrier
917 }
918 
919 
920 //====================================================================
921 template<typename AFIELD>
923  const AFIELD& w)
924 {
925 #pragma omp barrier
926 
927  real_t *vp = v.ptr(0);
928  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
929 
930  int ith = ThreadManager::get_thread_id();
931 
932  if (ith == 0)
934  vp, wp, m_Ns, m_Nsize);
935 
936 #pragma omp barrier
937 }
938 
939 
940 //====================================================================
941 template<typename AFIELD>
943  const AFIELD& w,
944  const int ieo)
945 { // ieo = 0: even < odd, ieo = 1: odd <- even
946 #pragma omp barrier
947 
948  int ith = ThreadManager::get_thread_id();
949  if(ith == 0) START_TIMER(timer_mult_Deo);
950 
951  int nth = ThreadManager::get_num_threads();
952  int ith_kernel = 0;
953  if(nth > 1) ith_kernel = 1;
954 
955  real_t *vp = v.ptr(0);
956  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
957  real_t *yp = m_w1.ptr(0); // working vector
958  real_t *up = m_Ueo.ptr(0);
959  int jeo = (m_Ieo_origin + ieo) % 2;
960 
961  if(do_comm_any > 0 && ith == 0){
962  START_TIMER(timer_comm_recv_start);
963  chset_recv.start();
964  STOP_TIMER(timer_comm_recv_start);
965  }
966 
967  if (ith == ith_kernel){
968 
970  yp, wp, m_mq, m_M0, m_Ns,
971  &m_b[0], &m_c[0], m_alpha, m_Nsize);
972 
973  if (do_comm_any > 0) {
974  START_TIMER(timer_pack);
975  real_t *buf1_xp = (real_t *)chsend_dn[0].ptr();
976  real_t *buf1_xm = (real_t *)chsend_up[0].ptr();
977  real_t *buf1_yp = (real_t *)chsend_dn[1].ptr();
978  real_t *buf1_ym = (real_t *)chsend_up[1].ptr();
979  real_t *buf1_zp = (real_t *)chsend_dn[2].ptr();
980  real_t *buf1_zm = (real_t *)chsend_up[2].ptr();
981  real_t *buf1_tp = (real_t *)chsend_dn[3].ptr();
982  real_t *buf1_tm = (real_t *)chsend_up[3].ptr();
983 
985  buf1_xp, buf1_xm, buf1_yp, buf1_ym,
986  buf1_zp, buf1_zm, buf1_tp, buf1_tm,
987  up, yp,
988  m_Ns, m_bc, m_Nsize, do_comm, ieo, jeo, 0);
989  STOP_TIMER(timer_pack);
990  }
991  }
992 
993 #pragma omp barrier
994 
995  if(do_comm_any > 0 && ith == 0){
996  START_TIMER(timer_comm);
997  START_TIMER(timer_comm_send_start);
998  chset_send.start();
999  STOP_TIMER(timer_comm_send_start);
1000  }
1001 
1002  if (ith == ith_kernel){
1003  START_TIMER(timer_bulk);
1004  if(m_impl == "5d"){
1006  vp, up, yp, m_Ns, m_bc2,
1007  m_Nsize, do_comm, ieo, jeo, 0);
1008  }else{
1009 #ifdef USE_DOMAINWALL_5DIN_4D_KERNEL
1011  vp, up, yp, m_Ns, m_bc2,
1012  m_Nsize, do_comm, ieo, jeo, 0);
1013 #else
1014  vout.crucial(m_vl, "%s: 4D kernel not compiled\n",
1015  class_name.c_str());
1016  exit(EXIT_FAILURE);
1017 #endif
1018  }
1019  STOP_TIMER(timer_bulk);
1020  }
1021 
1022  if(do_comm_any > 0 && ith == 0){
1023  START_TIMER(timer_comm_recv_wait);
1024  chset_recv.wait();
1025  STOP_TIMER(timer_comm_recv_wait);
1026  }
1027 
1028 #pragma omp barrier
1029 
1030  if(do_comm_any > 0 && ith == 0){
1031  STOP_TIMER(timer_comm);
1032  }
1033 
1034  if(do_comm_any > 0 && ith == ith_kernel){
1035  START_TIMER(timer_boundary);
1036 
1037  real_t *buf2_xp = (real_t *)chrecv_up[0].ptr();
1038  real_t *buf2_xm = (real_t *)chrecv_dn[0].ptr();
1039  real_t *buf2_yp = (real_t *)chrecv_up[1].ptr();
1040  real_t *buf2_ym = (real_t *)chrecv_dn[1].ptr();
1041  real_t *buf2_zp = (real_t *)chrecv_up[2].ptr();
1042  real_t *buf2_zm = (real_t *)chrecv_dn[2].ptr();
1043  real_t *buf2_tp = (real_t *)chrecv_up[3].ptr();
1044  real_t *buf2_tm = (real_t *)chrecv_dn[3].ptr();
1045 
1047  vp, up, yp,
1048  buf2_xp, buf2_xm, buf2_yp, buf2_ym,
1049  buf2_zp, buf2_zm, buf2_tp, buf2_tm,
1050  m_Ns, m_bc,
1051  m_Nsize, do_comm, ieo, jeo);
1052 
1053  STOP_TIMER(timer_boundary);
1054  }
1055  if(do_comm_any > 0 && ith == 0){
1056  START_TIMER(timer_comm_send_wait);
1057  chset_send.wait();
1058  STOP_TIMER(timer_comm_send_wait);
1059  }
1060 
1061 #pragma omp barrier
1062 
1063  if(ith == 0) STOP_TIMER(timer_mult_Deo);
1064 
1065 }
1066 
1067 //====================================================================
1068 template<typename AFIELD>
1070  const AFIELD& w,
1071  const int ieo)
1072  // ieo = 0: even < odd, ieo = 1: odd <- even
1073 {
1074 #pragma omp barrier
1075 
1076  int ith = ThreadManager::get_thread_id();
1077  if(ith == 0) START_TIMER(timer_mult_Deo);
1078 
1079  int nth = ThreadManager::get_num_threads();
1080  int ith_kernel = 0;
1081  if(nth > 1) ith_kernel = 1;
1082 
1083  real_t *vp = v.ptr(0);
1084  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1085  real_t *yp = m_w1.ptr(0);
1086  real_t *up = m_Ueo.ptr(0);
1087  int jeo = (m_Ieo_origin + ieo) % 2;
1088 
1089  if(do_comm_any > 0 && ith == 0){
1090  START_TIMER(timer_comm_recv_start);
1091  chset_recv.start();
1092  STOP_TIMER(timer_comm_recv_start);
1093  }
1094 
1095  if (ith == ith_kernel){
1096 
1097  if (do_comm_any > 0) {
1098  START_TIMER(timer_pack);
1099 
1100  real_t *buf1_xp = (real_t *)chsend_dn[0].ptr();
1101  real_t *buf1_xm = (real_t *)chsend_up[0].ptr();
1102  real_t *buf1_yp = (real_t *)chsend_dn[1].ptr();
1103  real_t *buf1_ym = (real_t *)chsend_up[1].ptr();
1104  real_t *buf1_zp = (real_t *)chsend_dn[2].ptr();
1105  real_t *buf1_zm = (real_t *)chsend_up[2].ptr();
1106  real_t *buf1_tp = (real_t *)chsend_dn[3].ptr();
1107  real_t *buf1_tm = (real_t *)chsend_up[3].ptr();
1108 
1110  buf1_xp, buf1_xm, buf1_yp, buf1_ym,
1111  buf1_zp, buf1_zm, buf1_tp, buf1_tm,
1112  up, wp,
1113  m_Ns, m_bc, m_Nsize, do_comm, ieo, jeo, 1);
1114  STOP_TIMER(timer_pack);
1115  }
1116  }
1117 #pragma omp barrier
1118 
1119  if(do_comm_any > 0 && ith == 0){
1120  START_TIMER(timer_comm);
1121  START_TIMER(timer_comm_send_start);
1122  chset_send.start();
1123  STOP_TIMER(timer_comm_send_start);
1124  }
1125 
1126  if(ith == ith_kernel){
1127  START_TIMER(timer_bulk);
1128  if(m_impl == "5d"){
1130  yp, up, wp, m_Ns, m_bc2,
1131  m_Nsize, do_comm, ieo, jeo, 1);
1132  }else{
1133 #ifdef USE_DOMAINWALL_5DIN_4D_KERNEL
1135  yp, up, wp, m_Ns, m_bc2,
1136  m_Nsize, do_comm, ieo, jeo, 1);
1137 #else
1138  vout.crucial(m_vl, "%s: 4D kernel not compiled\n",
1139  class_name.c_str());
1140  exit(EXIT_FAILURE);
1141 #endif
1142  }
1143  STOP_TIMER(timer_bulk);
1144  }
1145 
1146  if(do_comm_any > 0 && ith == 0){
1147  START_TIMER(timer_comm_recv_wait);
1148  chset_recv.wait();
1149  STOP_TIMER(timer_comm_recv_wait);
1150  }
1151 
1152 #pragma omp barrier
1153 
1154  if(do_comm_any > 0 && ith == 0){
1155  STOP_TIMER(timer_comm);
1156  }
1157 
1158  if (ith == ith_kernel) {
1159  if (do_comm_any > 0) {
1160  START_TIMER(timer_boundary);
1161 
1162  real_t *buf2_xp = (real_t *)chrecv_up[0].ptr();
1163  real_t *buf2_xm = (real_t *)chrecv_dn[0].ptr();
1164  real_t *buf2_yp = (real_t *)chrecv_up[1].ptr();
1165  real_t *buf2_ym = (real_t *)chrecv_dn[1].ptr();
1166  real_t *buf2_zp = (real_t *)chrecv_up[2].ptr();
1167  real_t *buf2_zm = (real_t *)chrecv_dn[2].ptr();
1168  real_t *buf2_tp = (real_t *)chrecv_up[3].ptr();
1169  real_t *buf2_tm = (real_t *)chrecv_dn[3].ptr();
1170 
1172  yp, up, vp,
1173  buf2_xp, buf2_xm, buf2_yp, buf2_ym,
1174  buf2_zp, buf2_zm, buf2_tp, buf2_tm,
1175  m_Ns, m_bc,
1176  m_Nsize, do_comm, ieo, jeo);
1177  STOP_TIMER(timer_boundary);
1178 
1179  }
1181  vp, yp, m_mq, m_M0, m_Ns,
1182  &m_b[0], &m_c[0], m_alpha, m_Nsize);
1183 
1184  }
1185  if(do_comm_any > 0 && ith == 0){
1186  START_TIMER(timer_comm_send_wait);
1187  chset_send.wait();
1188  STOP_TIMER(timer_comm_send_wait);
1189  }
1190 
1191 #pragma omp barrier
1192 
1193  if(ith == 0) STOP_TIMER(timer_mult_Deo);
1194 
1195 }
1196 
1197 //====================================================================
1198 template<typename AFIELD>
1200  const int ieo)
1201 {
1202 #pragma omp barrier
1203 
1204  int ith = ThreadManager::get_thread_id();
1205  if (ith == 0){
1206  real_t *vp = v.ptr(0);
1207  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1208 
1210  vp, wp, m_mq, m_M0, m_Ns,
1211  &m_b[0], &m_c[0], m_alpha, m_Nsize);
1212  }
1213 
1214 #pragma omp barrier
1215 
1216 }
1217 
1218 
1219 //====================================================================
1220 template<typename AFIELD>
1222  const int ieo)
1223 {
1224 #pragma omp barrier
1225 
1226  int ith = ThreadManager::get_thread_id();
1227  if (ith == 0){
1228  real_t *vp = v.ptr(0);
1229  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1230 
1232  vp, wp, m_mq, m_M0, m_Ns,
1233  &m_b[0], &m_c[0], m_alpha, m_Nsize);
1234  }
1235 
1236 #pragma omp barrier
1237 
1238 
1239 }
1240 
1241 //====================================================================
1242 template<typename AFIELD>
1244  const AFIELD& w,
1245  const int ieo)
1246 {
1247 #pragma omp barrier
1248 
1249  int ith = ThreadManager::get_thread_id();
1250 
1251  if(ith == 0) START_TIMER(timer_mult_Dee_inv);
1252 
1253 #ifdef USE_DOMAINWALL_5DIN_EE_MATINV_KERNEL
1254  if (ith == 0){
1255  real_t *vp = v.ptr(0);
1256  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1257 
1258  if(m_impl == "5d"){
1260  vp, wp, 1, m_Ns, &m_mat_inv[0], m_Nsize);
1261  }else{
1263  vp, wp, 1, m_Ns, &m_mat_inv[0], m_Nsize);
1264  }
1265  }
1266 #else
1267  vout.crucial(m_vl, "%s: Dee inverse matrix mult not compiled\n",
1268  class_name.c_str());
1269  exit(EXIT_FAILURE);
1270 #endif
1271 
1272 #pragma omp barrier
1273 
1274  if(ith == 0) STOP_TIMER(timer_mult_Dee_inv);
1275 
1276 }
1277 
1278 //====================================================================
1279 template<typename AFIELD>
1281  const AFIELD& w,
1282  const int ieo)
1283 {
1284 #pragma omp barrier
1285 
1286  int ith = ThreadManager::get_thread_id();
1287 
1288  if(ith == 0) START_TIMER(timer_mult_Dee_inv);
1289 
1290 #ifdef USE_DOMAINWALL_5DIN_EE_MATINV_KERNEL
1291  if (ith == 0){
1292  real_t *vp = v.ptr(0);
1293  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1294 
1295  if(m_impl == "5d"){
1297  vp, wp, -1, m_Ns, &m_mat_inv[0], m_Nsize);
1298  }else{
1300  vp, wp, -1, m_Ns, &m_mat_inv[0], m_Nsize);
1301  }
1302  }
1303 #else
1304  vout.crucial(m_vl, "%s: Dee inverse matrix mult not compiled\n",
1305  class_name.c_str());
1306  exit(EXIT_FAILURE);
1307 #endif
1308 
1309 #pragma omp barrier
1310 
1311  if(ith == 0) STOP_TIMER(timer_mult_Dee_inv);
1312 
1313 }
1314 
1315 //====================================================================
1316 template<typename AFIELD>
1318  const int ieo)
1319 {
1320 #pragma omp barrier
1321 
1322  LU_inv(v, w);
1323 
1324 #pragma omp barrier
1325 
1326 }
1327 
1328 //====================================================================
1329 template<typename AFIELD>
1331  const int ieo)
1332 {
1333 #pragma omp barrier
1334 
1335  LUdag_inv(v, w);
1336 
1337 #pragma omp barrier
1338 
1339 }
1340 
1341 //====================================================================
1342 template<typename AFIELD>
1344 {
1345 
1346  int ith = ThreadManager::get_thread_id();
1347 
1348  if (ith == 0){
1349 
1350  START_TIMER(timer_mult_Dee_inv);
1351 
1352  real_t *vp = v.ptr(0);
1353  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1354 
1356  vp, wp, m_Ns, m_Nsize,
1357  &m_e[0], &m_f[0], &m_dpinv[0], &m_dm[0], m_alpha);
1358 
1359  }
1360 #pragma omp barrier
1361 
1362  if(ith == 0) STOP_TIMER(timer_mult_Dee_inv);
1363 
1364 }
1365 
1366 
1367 //====================================================================
1368 template<typename AFIELD>
1370 {
1371  int ith = ThreadManager::get_thread_id();
1372 
1373  if (ith == 0){
1374 
1375  START_TIMER(timer_mult_Dee_inv);
1376 
1377  real_t *vp = v.ptr(0);
1378  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
1379 
1381  vp, wp, m_Ns, m_Nsize,
1382  &m_e[0], &m_f[0], &m_dpinv[0], &m_dm[0], m_alpha);
1383 
1384  }
1385 #pragma omp barrier
1386 
1387  if(ith == 0) STOP_TIMER(timer_mult_Dee_inv);
1388 
1389 }
1390 
1391 
1392 //====================================================================
1393 template<typename AFIELD>
1395 {
1396  // flop counting of this class is not confirmed.
1397 
1398  int Lvol = CommonParameters::Lvol();
1399  double vsite = static_cast<double>(Lvol);
1400  double vNs = static_cast<double>(m_Ns);
1401 
1402  double flop_Wilson_hop;
1403  double flop_LU_inv;
1404  double flop_pre;
1405  int Nc = CommonParameters::Nc();
1406  int Nd = CommonParameters::Nd();
1407  if (m_repr == "Dirac") {
1408  flop_Wilson_hop = static_cast<double>(
1409  vNs * Nc * Nd * (6 * (4 * Nc + 2) + 2 * (4 * Nc + 1))) * vsite;
1410  flop_pre = static_cast<double>(vNs * Nc * Nd * 14) * vsite;
1411  // chiral projectin:
1412  // FLOP as (1 +- gm5), because the factor (1/2) can be absorbed
1413  // into the other coefficients
1414  //
1415  // L_inv: chiral projection: Ns [each: 2*Nc*Nd]
1416  // axpy: 2(Ns-1) [each: 4*Nc*Nd]
1417  // U_inv: chiral projection: Ns
1418  // axpy: 2(Ns-1)
1419  // scal: Ns [each: 2*Nc*Nd]
1420  flop_LU_inv = static_cast<double>(Nc * Nd * (22 * vNs -16)) * vsite;
1421  } else if (m_repr == "Chiral") {
1422  flop_Wilson_hop = static_cast<double>(
1423  vNs * Nc * Nd * (8 * (4 * Nc + 2))) * vsite;
1424  flop_pre = static_cast<double>(vNs * Nc * Nd * 6) * vsite;
1425  // no FLOP is needed for chiral projection
1426  // L_inv: chiral projection: Ns
1427  // axpy: 2(Ns-1) [each: 2*Nc*Nd for chirality+/-]
1428  // U_inv: chiral projection: Ns
1429  // axpy: 2(Ns-1)
1430  // scal: Ns [each: 2*Nc*Nd]
1431  flop_LU_inv = static_cast<double>(Nc * Nd * (10 * vNs -8)) * vsite;
1432  }
1433  double flop_axpy = static_cast<double>(vNs * 2 * Nc * Nd) * vsite;
1434 
1435  // double axpy1 = static_cast<double>(2 * m_NinF);
1436  // double scal1 = static_cast<double>(1 * m_NinF);
1437  // double flop_DW = vNs * (flop_Wilson + vsite*(6*axpy1 + 2*scal1));
1438  // In Ddag case, flop_Wilson + 7 axpy which equals flop_DW.
1439  double flop_Deo = flop_Wilson_hop + flop_pre;
1440 
1441 
1442  // double flop_LU_inv = 2.0 * vsite *
1443  // ((3.0*axpy1 + scal1)*(vNs-1.0) + axpy1 + 2.0*scal1);
1444 
1445  double flop = 0.0;
1446  if ((mode == "D") || (mode == "Ddag")) {
1447  flop = flop_Deo + flop_LU_inv + 0.5 * flop_axpy; // axpy is for 1 - (Deo etc)
1448  } else if (mode == "DdagD") {
1449  flop = 2 * flop_Deo + 2 * flop_LU_inv + flop_axpy;
1450  } else if ((mode == "Dee_inv") || (mode == "Doo_inv")) {
1451  flop = 0.5 * flop_LU_inv; // 0.5 is for even-odd
1452  } else if ((mode == "Dee") || (mode == "Doo")) {
1453  flop = 0.5 * flop_pre;
1454  } else if ((mode == "Doe") || (mode == "Doe")) {
1455  flop = 0.5 * flop_Deo; // 0.5 is for even-odd
1456  } else {
1457  vout.crucial(m_vl, "Error at %s: input mode is undefined: %s.\n",
1458  class_name.c_str(), mode.c_str());
1459  exit(EXIT_FAILURE);
1460  }
1461 
1462  return flop;
1463 }
1464 
1465 //============================================================END=====
AFopr_Domainwall_5din_eo::flop_count
double flop_count()
this returns the number of floating point number operations.
Definition: afopr_Domainwall_5din_eo.h:189
BridgeACC::mult_domainwall_5din_mult_gm5_dirac
void mult_domainwall_5din_mult_gm5_dirac(double *RESTRICT vp, double *RESTRICT wp, int Ns, int *Nsize)
Definition: mult_Domainwall_5din_openacc-inc.h:293
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
AFopr_Domainwall_5din_eo::Ddag_alt
void Ddag_alt(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:861
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
convert_gauge
void convert_gauge(INDEX &index, AFIELD &v, const Field &w)
Definition: afield-inc.h:223
AFopr_Domainwall_5din_eo::init
void init(const Parameters &params)
initial setup.
Definition: afopr_Domainwall_5din_eo-tmpl.h:31
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
BridgeACC::mult_domainwall_5din_mult_R
void mult_domainwall_5din_mult_R(double *RESTRICT vp, double *RESTRICT wp, int Ns, int *Nsize)
Definition: mult_Domainwall_5din_openacc-inc.h:341
START_TIMER
#define START_TIMER(var_timer)
Definition: afopr_Domainwall_5din_eo-tmpl.h:22
AFopr_Domainwall_5din_eo::set_matrix5d_inverse
void set_matrix5d_inverse()
set parameters for preconditioning.
Definition: afopr_Domainwall_5din_eo-tmpl.h:481
ThreadManager::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager.cpp:246
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
BridgeACC::afield_tidyup
void afield_tidyup(double *data, const int size)
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
AFopr_Domainwall_5din_eo::set_precond_parameters
void set_precond_parameters()
set parameters for preconditioning.
Definition: afopr_Domainwall_5din_eo-tmpl.h:440
AIndex_lex
Definition: aindex_lex_base.h:17
AFopr_Domainwall_5din_eo::mult
void mult(AFIELD &v, const AFIELD &w)
multiplies fermion operator to a given field.
Definition: afopr_Domainwall_5din_eo-tmpl.h:634
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:518
BridgeACC::mult_domainwall_5din_ee_5dirdag_dirac
void mult_domainwall_5din_ee_5dirdag_dirac(double *RESTRICT vp, double *RESTRICT wp, double mq, double M0, int Ns, double *b, double *c, double alpha, int *Nsize)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:230
NVCD
#define NVCD
Definition: define_params_SU3.h:20
BridgeACC::mult_domainwall_5din_eo_hopb_dirac_5d
void mult_domainwall_5din_eo_hopb_dirac_5d(double *RESTRICT vp, double *RESTRICT up, double *RESTRICT wp, int Ns, int *bc, int *Nsize, int *do_comm, int ieo, int jeo, int jgm5)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:450
BridgeACC::mult_domainwall_5din_mult_gm5R_dirac
void mult_domainwall_5din_mult_gm5R_dirac(double *RESTRICT vp, double *RESTRICT wp, int Ns, int *Nsize)
Definition: mult_Domainwall_5din_openacc-inc.h:385
AFopr_Domainwall_5din_eo::mult_dag
void mult_dag(AFIELD &v, const AFIELD &w)
hermitian conjugate of mult.
Definition: afopr_Domainwall_5din_eo-tmpl.h:674
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:281
Field::check_size
bool check_size(const int nin, const int nvol, const int nex) const
checking size parameters. [23 May 2016 H.Matsufuru]
Definition: field.h:135
AFopr_Domainwall_5din_eo::convert
void convert(AFIELD &, const Field &)
convert Field to AField for this class.
Definition: afopr_Domainwall_5din_eo-tmpl.h:555
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
AFopr_Domainwall_5din_eo::set_mode
void set_mode(std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: afopr_Domainwall_5din_eo-tmpl.h:620
AFopr_Domainwall_5din_eo
Domain-wall fermion operator with even-odd site index.
Definition: afopr_Domainwall_5din_eo.h:41
BridgeACC::mult_domainwall_5din_eo_5dir_dirac
void mult_domainwall_5din_eo_5dir_dirac(double *RESTRICT yp, double *RESTRICT wp, double mq, double M0, int Ns, double *b, double *c, double alpha, int *Nsize)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:122
BridgeACC::mult_domainwall_5din_LUinv_dirac
void mult_domainwall_5din_LUinv_dirac(double *RESTRICT vp, double *RESTRICT wp, int Ns, int *Nsize, double *e, double *f, double *dpinv, double *dm, double alpha)
Definition: mult_Domainwall_5din_LUinv_openacc-inc.h:14
AFopr_Domainwall_5din_eo::reverse
void reverse(Field &, const AFIELD &)
reverse AField to Field.
Definition: afopr_Domainwall_5din_eo-tmpl.h:587
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
AFopr_Domainwall_5din_eo::Ddag_ee_inv
void Ddag_ee_inv(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1330
BridgeACC::mult_domainwall_5din_eo_5dirdag_dirac
void mult_domainwall_5din_eo_5dirdag_dirac(double *RESTRICT vp, double *RESTRICT yp, double mq, double M0, int Ns, double *b, double *c, double alpha, int *Nsize)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:346
AFopr_Domainwall_5din_eo::mult_gm5R
void mult_gm5R(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:922
AFopr_Domainwall_5din_eo::Ddag
void Ddag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:806
STOP_TIMER
#define STOP_TIMER(var_timer)
Definition: afopr_Domainwall_5din_eo-tmpl.h:23
Timer
Definition: timer.h:31
AFopr_Domainwall_5din_eo::setup_channels
void setup_channels()
Definition: afopr_Domainwall_5din_eo-tmpl.h:194
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
BridgeACC::mult_domainwall_5din_eo_hop2_dirac
void mult_domainwall_5din_eo_hop2_dirac(double *RESTRICT vp, double *RESTRICT up, double *RESTRICT wp, double *RESTRICT buf2_xp, double *RESTRICT buf2_xm, double *RESTRICT buf2_yp, double *RESTRICT buf2_ym, double *RESTRICT buf2_zp, double *RESTRICT buf2_zm, double *RESTRICT buf2_tp, double *RESTRICT buf2_tm, int Ns, int *bc, int *Nsize, int *do_comm, int ieo, int jeo)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:889
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:300
BridgeACC::mult_domainwall_5din_eo_hop1_dirac
void mult_domainwall_5din_eo_hop1_dirac(double *RESTRICT buf1_xp, double *RESTRICT buf1_xm, double *RESTRICT buf1_yp, double *RESTRICT buf1_ym, double *RESTRICT buf1_zp, double *RESTRICT buf1_zm, double *RESTRICT buf1_tp, double *RESTRICT buf1_tm, double *RESTRICT up, double *RESTRICT wp, int Ns, int *bc, int *Nsize, int *do_comm, int ieo, int jeo, int jgm5)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:553
AFopr_Domainwall_5din_eo::set_coefficients
void set_coefficients(const std::vector< real_t > b, const std::vector< real_t > c)
set coefficients if they depend in s.
Definition: afopr_Domainwall_5din_eo-tmpl.h:409
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
AFopr_Domainwall_5din_eo::mult_gm5
void mult_gm5(AFIELD &v, const AFIELD &w)
multiplies gamma_5 matrix.
Definition: afopr_Domainwall_5din_eo-tmpl.h:881
AFopr_Domainwall_5din_eo::Ddag_ee_inv_alt
void Ddag_ee_inv_alt(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1280
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
AFopr_Domainwall_5din_eo::D_ee_inv
void D_ee_inv(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1317
timer.h
AFopr_Domainwall_5din_eo::DdagD
void DdagD(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:771
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
BridgeACC::mult_domainwall_5din_ee_5dir_dirac
void mult_domainwall_5din_ee_5dir_dirac(double *RESTRICT vp, double *RESTRICT wp, double mq, double M0, int Ns, double *b, double *c, double alpha, int *Nsize)
Definition: mult_Domainwall_5din_eo_openacc-inc.h:14
Communicator::npe
static int npe(const int dir)
logical grid extent
Definition: communicator.cpp:112
afopr_Domainwall_5din_eo.h
AFopr_Domainwall_5din_eo::Ddag_ee
void Ddag_ee(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1221
AFopr_Domainwall_5din_eo::tidyup
void tidyup()
final tidyup.
Definition: afopr_Domainwall_5din_eo-tmpl.h:151
AFopr_Domainwall_5din_eo::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall_5din_eo-tmpl.h:310
BridgeACC::mult_domainwall_5din_ee_inv_dirac_5d
void mult_domainwall_5din_ee_inv_dirac_5d(double *RESTRICT vp, double *RESTRICT wp, int jd, int Ns, double *RESTRICT mat_inv, int *Nsize)
Definition: mult_Domainwall_5din_matinv_openacc-inc.h:14
AFopr_Domainwall_5din_eo::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall_5din_eo-tmpl.h:246
AFopr_Domainwall_5din_eo::D_eo
void D_eo(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:942
Field::nvol
int nvol() const
Definition: field.h:127
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
real_t
double real_t
Definition: bridgeACC_AField_double.cpp:14
AFopr_Domainwall_5din_eo::LU_inv
void LU_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1343
BridgeACC::mult_domainwall_5din_ee_inv_dirac_4d
void mult_domainwall_5din_ee_inv_dirac_4d(double *RESTRICT vp, double *RESTRICT wp, int jd, int Ns, double *RESTRICT mat_inv, int *Nsize)
Definition: mult_Domainwall_5din_matinv_openacc-inc.h:104
BridgeACC::mult_domainwall_5din_LUdaginv_dirac
void mult_domainwall_5din_LUdaginv_dirac(double *RESTRICT vp, double *RESTRICT wp, int Ns, int *Nsize, double *e, double *f, double *dpinv, double *dm, double alpha)
Definition: mult_Domainwall_5din_LUinv_openacc-inc.h:204
AFopr_Domainwall_5din_eo::mult_R
void mult_R(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:900
AFopr_Domainwall_5din_eo::set_config_impl
void set_config_impl(Field *U)
Definition: afopr_Domainwall_5din_eo-tmpl.h:544
Field::cmp
double cmp(const int jin, const int site, const int jex) const
Definition: field.h:143
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
AFopr_Domainwall_5din_eo::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: afopr_Domainwall_5din_eo-tmpl.h:514
AFopr_Domainwall_5din_eo::Ddag_eo
void Ddag_eo(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1069
AFopr_Domainwall_5din_eo::real_t
AFIELD::real_t real_t
Definition: afopr_Domainwall_5din_eo.h:44
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
AFopr_Domainwall_5din_eo::LUdag_inv
void LUdag_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1369
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
BridgeACC::mult_domainwall_5din_eo_hopb_dirac_4d
void mult_domainwall_5din_eo_hopb_dirac_4d(double *RESTRICT vp, double *RESTRICT up, double *RESTRICT wp, int Ns, int *bc, int *Nsize, int *do_comm, int ieo, int jeo, int jgm5)
Definition: mult_Domainwall_5din_eo_4d_openacc-inc.h:14
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
AFopr_Domainwall_5din_eo::D_alt
void D_alt(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:847
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
AFopr_Domainwall_5din_eo::D_ee
void D_ee(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1199
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
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
AFopr_Domainwall_5din_eo::D
void D(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:793
afopr_common_th-inc.h
AFopr_Domainwall_5din_eo::set_config_omp
void set_config_omp(Field *U)
Definition: afopr_Domainwall_5din_eo-tmpl.h:532
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
AFopr_Domainwall_5din_eo::DdagD_alt
void DdagD_alt(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din_eo-tmpl.h:826
Field
Container of Field-type object.
Definition: field.h:46
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
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
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:572
AFopr_Domainwall_5din_eo::D_ee_inv_alt
void D_ee_inv_alt(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Domainwall_5din_eo-tmpl.h:1243
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:216
BridgeACC::afield_init
void afield_init(double *data, const int size)