Bridge++  Ver. 2.0.4
afopr_Domainwall_5din-tmpl.h
Go to the documentation of this file.
1 
10 template<typename AFIELD>
11 const std::string AFopr_Domainwall_5din<AFIELD>::class_name
12  = "AFopr_Domainwall_5din";
13 //====================================================================
14 template<typename AFIELD>
15 void AFopr_Domainwall_5din<AFIELD>::init(const Parameters& params)
16 {
17  ThreadManager::assert_single_thread(class_name);
18 
19  int req_comm = 1; // set 1 if communication forced any time
20  //int req_comm = 0; // set 1 if communication forced any time
21 
22  std::string vlevel;
23  if (!params.fetch_string("verbose_level", vlevel)) {
24  m_vl = vout.set_verbose_level(vlevel);
25  } else {
26  m_vl = CommonParameters::Vlevel();
27  }
28 
29  vout.general(m_vl, "%s: construction\n", class_name.c_str());
30 
31  m_repr = "Dirac"; // now only the Dirac repr is available.
32 
33  std::string repr;
34  if (!params.fetch_string("gamma_matrix_type", repr)) {
35  if (repr != "Dirac") {
36  vout.crucial(" Error at %s: unsupported gamma-matrix type: %s\n",
37  class_name.c_str(), repr.c_str());
38  exit(EXIT_FAILURE);
39  }
40  }
41 
42  int Nc = CommonParameters::Nc();
43  if (Nc != 3) {
44  vout.crucial("%s: only applicable to Nc = 3\n",
45  class_name.c_str());
46  exit(EXIT_FAILURE);
47  }
48 
49  int Nd = CommonParameters::Nd();
50  m_Nvcd = 2 * Nc * Nd;
51  m_Ndf = 2 * Nc * Nc;
52 
53  m_Nvol = CommonParameters::Nvol();
54  m_Ndim = CommonParameters::Ndim();
55  m_Nx = CommonParameters::Nx();
56  m_Ny = CommonParameters::Ny();
57  m_Nz = CommonParameters::Nz();
58  m_Nt = CommonParameters::Nt();
59 
60  m_Nxv = m_Nx / VLENX;
61  m_Nyv = m_Ny / VLENY;
62  m_Nstv = m_Nvol / VLEN;
63 
64  if (VLENX * m_Nxv != m_Nx) {
65  vout.crucial(m_vl, "%s: Nx must be multiple of VLENX.\n",
66  class_name.c_str());
67  exit(EXIT_FAILURE);
68  }
69  if (VLENY * m_Nyv != m_Ny) {
70  vout.crucial(m_vl, "%s: Ny must be multiple of VLENY.\n",
71  class_name.c_str());
72  exit(EXIT_FAILURE);
73  }
74 
75  vout.general(m_vl, " VLENX = %2d Nxv = %d\n", VLENX, m_Nxv);
76  vout.general(m_vl, " VLENY = %2d Nyv = %d\n", VLENY, m_Nyv);
77  vout.general(m_vl, " VLEN = %2d Nstv = %d\n", VLEN, m_Nstv);
78 
79  m_Nsize[0] = m_Nxv;
80  m_Nsize[1] = m_Nyv;
81  m_Nsize[2] = m_Nz;
82  m_Nsize[3] = m_Nt;
83 
84  do_comm_any = 0;
85  for (int mu = 0; mu < m_Ndim; ++mu) {
86  do_comm[mu] = 1;
87  if ((req_comm == 0) && (Communicator::npe(mu) == 1)) do_comm[mu] = 0;
88  do_comm_any += do_comm[mu];
89  vout.general(" do_comm[%d] = %d\n", mu, do_comm[mu]);
90  }
91 
92  vout.increase_indent();
93 
94  m_Ns = 0; // temporary set
95  set_parameters(params);
96 
97  vout.decrease_indent();
98 
99  // gauge configuration.
100  m_U.reset(m_Ndf, m_Nvol, m_Ndim);
101 
102  m_Nbdsize.resize(m_Ndim);
103  int Nvst = (m_Nvcd / 2) * m_Ns * m_Nvol;
104  m_Nbdsize[0] = Nvst / m_Nx;
105  m_Nbdsize[1] = Nvst / m_Ny;
106  m_Nbdsize[2] = Nvst / m_Nz;
107  m_Nbdsize[3] = Nvst / m_Nt;
108 
109  setup_channels();
110 
111  vout.general(m_vl, "%s: construction finished.\n",
112  class_name.c_str());
113 }
114 
115 
116 //====================================================================
117 template<typename AFIELD>
118 void AFopr_Domainwall_5din<AFIELD>::tidyup()
119 {
120  // need to do nothing.
121 }
122 
123 
124 //====================================================================
125 template<typename AFIELD>
126 void AFopr_Domainwall_5din<AFIELD>::setup_channels()
127 {
128  chsend_up.resize(m_Ndim);
129  chrecv_up.resize(m_Ndim);
130  chsend_dn.resize(m_Ndim);
131  chrecv_dn.resize(m_Ndim);
132 
133  for (int mu = 0; mu < m_Ndim; ++mu) {
134  size_t Nvsize = m_Nbdsize[mu] * sizeof(real_t);
135 
136  chsend_dn[mu].send_init(Nvsize, mu, -1);
137  chsend_up[mu].send_init(Nvsize, mu, 1);
138 #ifdef USE_MPI
139  chrecv_up[mu].recv_init(Nvsize, mu, 1);
140  chrecv_dn[mu].recv_init(Nvsize, mu, -1);
141 #else
142  void *buf_up = (void *)chsend_dn[mu].ptr();
143  chrecv_up[mu].recv_init(Nvsize, mu, 1, buf_up);
144  void *buf_dn = (void *)chsend_up[mu].ptr();
145  chrecv_dn[mu].recv_init(Nvsize, mu, -1, buf_dn);
146 #endif
147 
148  if (do_comm[mu] == 1) {
149  chset_send.append(chsend_up[mu]);
150  chset_send.append(chsend_dn[mu]);
151  chset_recv.append(chrecv_up[mu]);
152  chset_recv.append(chrecv_dn[mu]);
153  }
154  }
155 }
156 
157 
158 //====================================================================
159 template<typename AFIELD>
160 void AFopr_Domainwall_5din<AFIELD>::set_parameters(
161  const Parameters& params)
162 {
163  string gmset_type;
164  double mq, M0;
165  int Ns;
166  std::vector<int> bc;
167  double b, c;
168 
169  int err_optional = 0;
170  err_optional += params.fetch_string("gamma_matrix_type", gmset_type);
171  if (!err_optional) {
172  if (gmset_type != m_repr) {
173  vout.crucial(m_vl, "%s: unsupported gamma_matrix_type: %s\n",
174  class_name.c_str(), gmset_type.c_str());
175  exit(EXIT_FAILURE);
176  }
177  }
178 
179  int err = 0;
180  err += params.fetch_double("quark_mass", mq);
181  err += params.fetch_double("domain_wall_height", M0);
182  err += params.fetch_int("extent_of_5th_dimension", Ns);
183  err += params.fetch_int_vector("boundary_condition", bc);
184  err += params.fetch_double("coefficient_b", b);
185  err += params.fetch_double("coefficient_c", c);
186 
187  if (err) {
188  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
189  class_name.c_str());
190  exit(EXIT_FAILURE);
191  }
192 
193  set_parameters(mq, M0, Ns, bc, b, c);
194 }
195 
196 
197 //====================================================================
198 template<typename AFIELD>
199 void AFopr_Domainwall_5din<AFIELD>::set_parameters(
200  const double mq,
201  const double M0,
202  const int Ns,
203  const vector<int> bc,
204  const double b,
205  const double c)
206 {
207  assert(bc.size() == m_Ndim);
208 
209 #pragma omp barrier
210 
211  int ith = ThreadManager::get_thread_id();
212 
213  if (ith == 0) {
214  m_M0 = real_t(M0);
215  m_mq = real_t(mq);
216  m_Ns = Ns;
217  m_NinF = m_Nvcd * m_Ns;
218 
219  if (m_boundary.size() != m_Ndim) m_boundary.resize(m_Ndim);
220  for (int mu = 0; mu < m_Ndim; ++mu) {
221  m_boundary[mu] = bc[mu];
222  }
223 
224  if (m_b.size() != m_Ns) {
225  m_b.resize(m_Ns);
226  m_c.resize(m_Ns);
227  }
228  for (int is = 0; is < m_Ns; ++is) {
229  m_b[is] = real_t(b);
230  m_c[is] = real_t(c);
231  }
232  }
233 
234  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
235  vout.general(m_vl, " mq = %8.4f\n", m_mq);
236  vout.general(m_vl, " M0 = %8.4f\n", m_M0);
237  vout.general(m_vl, " Ns = %4d\n", m_Ns);
238  for (int mu = 0; mu < m_Ndim; ++mu) {
239  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
240  }
241  vout.general(m_vl, " coefficients b = %16.10f c = %16.10f\n",
242  m_b[0], m_c[0]);
243 
244  set_precond_parameters();
245 
246  // 5D working vectors.
247  if (m_w1.nin() != m_NinF) {
248  if (ith == 0) {
249  m_w1.reset(m_NinF, m_Nvol, 1);
250  m_v1.reset(m_NinF, m_Nvol, 1);
251  m_v2.reset(m_NinF, m_Nvol, 1);
252  }
253  }
254 
255 #pragma omp barrier
256 }
257 
258 
259 //====================================================================
260 template<typename AFIELD>
261 void AFopr_Domainwall_5din<AFIELD>::set_precond_parameters()
262 {
263 #pragma omp barrier
264 
265  int ith = ThreadManager::get_thread_id();
266  if (ith == 0) {
267  if (m_dp.size() != m_Ns) {
268  m_dp.resize(m_Ns);
269  m_dm.resize(m_Ns);
270  m_dpinv.resize(m_Ns);
271  m_e.resize(m_Ns - 1);
272  m_f.resize(m_Ns - 1);
273  }
274 
275  for (int is = 0; is < m_Ns; ++is) {
276  m_dp[is] = 1.0 + m_b[is] * (4.0 - m_M0);
277  m_dm[is] = 1.0 - m_c[is] * (4.0 - m_M0);
278  }
279 
280  m_e[0] = m_mq * m_dm[m_Ns - 1] / m_dp[0];
281  m_f[0] = m_mq * m_dm[0];
282  for (int is = 1; is < m_Ns - 1; ++is) {
283  m_e[is] = m_e[is - 1] * m_dm[is - 1] / m_dp[is];
284  m_f[is] = m_f[is - 1] * m_dm[is] / m_dp[is - 1];
285  }
286 
287  m_g = m_e[m_Ns - 2] * m_dm[m_Ns - 2];
288 
289  for (int is = 0; is < m_Ns - 1; ++is) {
290  m_dpinv[is] = 1.0 / m_dp[is];
291  }
292  m_dpinv[m_Ns - 1] = 1.0 / (m_dp[m_Ns - 1] + m_g);
293  }
294 
295 #pragma omp barrier
296 }
297 
298 
299 //====================================================================
300 template<typename AFIELD>
301 void AFopr_Domainwall_5din<AFIELD>::set_coefficients(
302  const std::vector<double> vec_b,
303  const std::vector<double> vec_c)
304 {
305 #pragma omp barrier
306 
307  if ((vec_b.size() != m_Ns) || (vec_c.size() != m_Ns)) {
308  vout.crucial(m_vl, "%s: size of coefficient vectors incorrect.\n",
309  class_name.c_str());
310  }
311 
312  vout.general(m_vl, "%s: coefficient vectors are set:\n",
313  class_name.c_str());
314 
315  int ith = ThreadManager::get_thread_id();
316 
317  if (ith == 0) {
318  for (int is = 0; is < m_Ns; ++is) {
319  m_b[is] = real_t(vec_b[is]);
320  m_c[is] = real_t(vec_c[is]);
321  }
322  }
323 
324 #pragma omp barrier
325 
326  for (int is = 0; is < m_Ns; ++is) {
327  vout.general(m_vl, "b[%2d] = %16.10f c[%2d] = %16.10f\n",
328  is, m_b[is], is, m_c[is]);
329  }
330 
331  set_precond_parameters();
332 }
333 
334 
335 //====================================================================
336 template<typename AFIELD>
337 void AFopr_Domainwall_5din<AFIELD>::set_config(Field *u)
338 {
339  int nth = ThreadManager::get_num_threads();
340 
341  vout.detailed(m_vl, "%s: set_config is called: num_threads = %d\n",
342  class_name.c_str(), nth);
343 
344  if (nth > 1) {
345  set_config_impl(u);
346  } else {
347  set_config_omp(u);
348  }
349 
350  vout.detailed(m_vl, "%s: set_config finished\n", class_name.c_str());
351 }
352 
353 
354 //====================================================================
355 template<typename AFIELD>
356 void AFopr_Domainwall_5din<AFIELD>::set_config_omp(Field *u)
357 {
358  vout.detailed(m_vl, " set_config_omp is called.\n");
359 
360 #pragma omp parallel
361  {
362  set_config_impl(u);
363  }
364 }
365 
366 
367 //====================================================================
368 template<typename AFIELD>
369 void AFopr_Domainwall_5din<AFIELD>::set_config_impl(Field *u)
370 {
371 #pragma omp barrier
372 
373  AIndex_lex<real_t, AFIELD::IMPL> index_lex;
374  convert_gauge(index_lex, m_U, *u);
375 
376  QXS_Gauge::set_boundary(m_U, m_boundary);
377 
378 #pragma omp barrier
379 }
380 
381 
382 //====================================================================
383 template<typename AFIELD>
384 void AFopr_Domainwall_5din<AFIELD>::convert(AFIELD& v,
385  const Field& w)
386 { // the following implementation only applies Dirac repr.
387 #pragma omp barrier
388 
389  int ith, nth, isite, nsite;
390  set_threadtask_mult(ith, nth, isite, nsite, m_Nvol);
391 
392  AIndex_lex<real_t, AFIELD::IMPL> index;
393 
394  for (int site = isite; site < nsite; ++site) {
395  for (int is = 0; is < m_Ns; ++is) {
396  for (int ic = 0; ic < NC; ++ic) {
397  for (int id = 0; id < ND2; ++id) {
398  for (int iri = 0; iri < 2; ++iri) {
399  int in_org = iri + 2 * (ic + NC * id);
400  int in_alt = iri + 2 * (id + ND * ic) + NVCD * is;
401  v.set(index.idx(in_alt, m_NinF, site, 0),
402  real_t(w.cmp(in_org, site, is)));
403  }
404  }
405 
406  for (int id = ND2; id < ND; ++id) {
407  for (int iri = 0; iri < 2; ++iri) {
408  int in_org = iri + 2 * (ic + NC * id);
409  int in_alt = iri + 2 * (id + ND * ic) + NVCD * is;
410  v.set(index.idx(in_alt, m_NinF, site, 0),
411  real_t(-w.cmp(in_org, site, is)));
412  }
413  }
414  }
415  }
416  } // site-llop
417 
418 #pragma omp barrier
419 }
420 
421 
422 //====================================================================
423 template<typename AFIELD>
424 void AFopr_Domainwall_5din<AFIELD>::reverse(Field& v,
425  const AFIELD& w)
426 { // the following implementation only applies Dirac repr.
427 #pragma omp barrier
428 
429  int ith, nth, isite, nsite;
430  set_threadtask_mult(ith, nth, isite, nsite, m_Nvol);
431 
432  AIndex_lex<real_t, AFIELD::IMPL> index;
433 
434  for (int site = isite; site < nsite; ++site) {
435  for (int is = 0; is < m_Ns; ++is) {
436  for (int ic = 0; ic < NC; ++ic) {
437  for (int id = 0; id < ND2; ++id) {
438  for (int iri = 0; iri < 2; ++iri) {
439  int in_org = iri + 2 * (ic + NC * id);
440  int in_alt = iri + 2 * (id + ND * ic) + NVCD * is;
441  v.set(in_org, site, is,
442  double( w.cmp(index.idx(in_alt, m_NinF, site, 0))));
443  }
444  }
445 
446  for (int id = ND2; id < ND; ++id) {
447  for (int iri = 0; iri < 2; ++iri) {
448  int in_org = iri + 2 * (ic + NC * id);
449  int in_alt = iri + 2 * (id + ND * ic) + NVCD * is;
450  v.set(in_org, site, is,
451  double( -w.cmp(index.idx(in_alt, m_NinF, site, 0))));
452  }
453  }
454  }
455  }
456  } // site-loopa
457 
458 #pragma omp barrier
459 }
460 
461 
462 //====================================================================
463 template<typename AFIELD>
464 void AFopr_Domainwall_5din<AFIELD>::set_mode(std::string mode)
465 {
466 #pragma omp barrier
467 
468  int ith = ThreadManager::get_thread_id();
469  if (ith == 0) m_mode = mode;
470 
471 #pragma omp barrier
472 }
473 
474 
475 //====================================================================
476 template<typename AFIELD>
477 void AFopr_Domainwall_5din<AFIELD>::mult(AFIELD& v, const AFIELD& w,
478  std::string mode)
479 {
480  assert(w.check_size(m_NinF, m_Nvol, 1));
481  assert(v.check_size(m_NinF, m_Nvol, 1));
482 
483  if (mode == "D") {
484  D(v, w);
485  } else if (mode == "Ddag") {
486  Ddag(v, w);
487  } else if (mode == "DdagD") {
488  DdagD(v, w);
489  } else if (mode == "D_prec") {
490  D_prec(v, w);
491  } else if (mode == "Ddag_prec") {
492  Ddag_prec(v, w);
493  } else if (mode == "DdagD_prec") {
494  DdagD_prec(v, w);
495  } else if (mode == "Prec") {
496  Prec(v, w);
497  } else if (mode == "Precdag") {
498  Precdag(v, w);
499  } else {
500  vout.crucial(m_vl, "%s: undefined mode = %s\n",
501  class_name.c_str(), mode.c_str());
502  exit(EXIT_FAILURE);
503  }
504 }
505 
506 
507 //====================================================================
508 template<typename AFIELD>
509 void AFopr_Domainwall_5din<AFIELD>::mult_dag(AFIELD& v,
510  const AFIELD& w,
511  std::string mode)
512 {
513  assert(w.check_size(m_NinF, m_Nvol, 1));
514  assert(v.check_size(m_NinF, m_Nvol, 1));
515 
516  if (mode == "D") {
517  Ddag(v, w);
518  } else if (mode == "Ddag") {
519  D(v, w);
520  } else if (mode == "DdagD") {
521  DdagD(v, w);
522  } else if (mode == "D_prec") {
523  Ddag_prec(v, w);
524  } else if (mode == "Ddag_prec") {
525  D_prec(v, w);
526  } else if (mode == "DdagD_prec") {
527  DdagD_prec(v, w);
528  } else if (mode == "Prec") {
529  Precdag(v, w);
530  } else if (mode == "Precdag") {
531  Prec(v, w);
532  } else {
533  std::cout << "mode undeifined in AFopr_Domainwall_5din.\n";
534  abort();
535  }
536 }
537 
538 
539 //====================================================================
540 template<typename AFIELD>
541 void AFopr_Domainwall_5din<AFIELD>::mult_gm5(AFIELD& v, const AFIELD& w)
542 {
543  real_t *vp = v.ptr(0);
544  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
545 
546 #pragma omp barrier
547 
548  BridgeQXS::mult_domainwall_5din_mult_gm5_dirac(vp, wp, m_Ns, m_Nsize);
549 
550 #pragma omp barrier
551 }
552 
553 
554 //====================================================================
555 template<typename AFIELD>
556 void AFopr_Domainwall_5din<AFIELD>::DdagD(AFIELD& v, const AFIELD& w)
557 {
558  D(m_v1, w);
559  Ddag(v, m_v1);
560 }
561 
562 
563 //====================================================================
564 template<typename AFIELD>
565 void AFopr_Domainwall_5din<AFIELD>::DdagD_prec(AFIELD& v,
566  const AFIELD& w)
567 {
568  D_prec(m_v1, w);
569  Ddag_prec(v, m_v1);
570 }
571 
572 
573 //====================================================================
574 template<typename AFIELD>
575 void AFopr_Domainwall_5din<AFIELD>::D_prec(AFIELD& v,
576  const AFIELD& w)
577 {
578 #pragma omp barrier
579  L_inv(v, w);
580  U_inv(m_v2, v);
581  D(v, m_v2);
582 }
583 
584 
585 //====================================================================
586 template<typename AFIELD>
587 void AFopr_Domainwall_5din<AFIELD>::Ddag_prec(AFIELD& v,
588  const AFIELD& w)
589 {
590  Ddag(v, w);
591  Udag_inv(m_v2, v);
592  Ldag_inv(v, m_v2);
593 }
594 
595 
596 //====================================================================
597 template<typename AFIELD>
598 void AFopr_Domainwall_5din<AFIELD>::Prec(AFIELD& v,
599  const AFIELD& w)
600 {
601 #pragma omp barrier
602  L_inv(m_v2, w);
603  U_inv(v, m_v2);
604 }
605 
606 
607 //====================================================================
608 template<typename AFIELD>
609 void AFopr_Domainwall_5din<AFIELD>::Precdag(AFIELD& v,
610  const AFIELD& w)
611 {
612 #pragma omp barrier
613  Udag_inv(m_v2, w);
614  Ldag_inv(v, m_v2);
615 }
616 
617 
618 //====================================================================
619 template<typename AFIELD>
620 void AFopr_Domainwall_5din<AFIELD>::D(AFIELD& v, const AFIELD& w)
621 {
622  mult_D(v, w);
623 }
624 
625 
626 //====================================================================
627 template<typename AFIELD>
628 void AFopr_Domainwall_5din<AFIELD>::Ddag(AFIELD& v, const AFIELD& w)
629 {
630  mult_Ddag(v, w);
631 }
632 
633 
634 //====================================================================
635 template<typename AFIELD>
636 void AFopr_Domainwall_5din<AFIELD>::mult_D(AFIELD& v, const AFIELD& w)
637 {
638  int Nin4 = VLEN * NVCD;
639  int Nin5 = Nin4 * m_Ns;
640 
641  real_t *vp = v.ptr(0);
642  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
643  real_t *yp = m_w1.ptr(0); // working vector
644  real_t *up = m_U.ptr(0);
645 
646  int ith = ThreadManager::get_thread_id();
647 
648 #pragma omp barrier
649 
650  BridgeQXS::mult_domainwall_5din_5dir_dirac(vp, yp, wp,
651  m_mq, m_M0, m_Ns, &m_boundary[0],
652  &m_b[0], &m_c[0],
653  m_Nsize, do_comm);
654 
655 #pragma omp barrier
656 
657  if (do_comm_any > 0) {
658  if (ith == 0) chset_recv.start();
659 
660  real_t *buf1_xp = (real_t *)chsend_dn[0].ptr();
661  real_t *buf1_xm = (real_t *)chsend_up[0].ptr();
662  real_t *buf1_yp = (real_t *)chsend_dn[1].ptr();
663  real_t *buf1_ym = (real_t *)chsend_up[1].ptr();
664  real_t *buf1_zp = (real_t *)chsend_dn[2].ptr();
665  real_t *buf1_zm = (real_t *)chsend_up[2].ptr();
666  real_t *buf1_tp = (real_t *)chsend_dn[3].ptr();
667  real_t *buf1_tm = (real_t *)chsend_up[3].ptr();
668 
669  BridgeQXS::mult_domainwall_5din_hop1_dirac(
670  buf1_xp, buf1_xm, buf1_yp, buf1_ym,
671  buf1_zp, buf1_zm, buf1_tp, buf1_tm,
672  up, yp,
673  m_mq, m_M0, m_Ns, &m_boundary[0],
674  m_Nsize, do_comm);
675 
676 #pragma omp barrier
677 
678  if (ith == 0) chset_send.start();
679  }
680 
681  BridgeQXS::mult_domainwall_5din_hopb_dirac(vp, up, yp,
682  m_mq, m_M0, m_Ns, &m_boundary[0],
683  &m_b[0], &m_c[0],
684  m_Nsize, do_comm);
685 
686  if (do_comm_any > 0) {
687  if (ith == 0) chset_recv.wait();
688 
689 #pragma omp barrier
690 
691  real_t *buf2_xp = (real_t *)chrecv_up[0].ptr();
692  real_t *buf2_xm = (real_t *)chrecv_dn[0].ptr();
693  real_t *buf2_yp = (real_t *)chrecv_up[1].ptr();
694  real_t *buf2_ym = (real_t *)chrecv_dn[1].ptr();
695  real_t *buf2_zp = (real_t *)chrecv_up[2].ptr();
696  real_t *buf2_zm = (real_t *)chrecv_dn[2].ptr();
697  real_t *buf2_tp = (real_t *)chrecv_up[3].ptr();
698  real_t *buf2_tm = (real_t *)chrecv_dn[3].ptr();
699 
700  BridgeQXS::mult_domainwall_5din_hop2_dirac(vp, up, yp,
701  buf2_xp, buf2_xm, buf2_yp, buf2_ym,
702  buf2_zp, buf2_zm, buf2_tp, buf2_tm,
703  m_mq, m_M0, m_Ns, &m_boundary[0],
704  m_Nsize, do_comm);
705 
706  if (ith == 0) chset_send.wait();
707  }
708 
709 #pragma omp barrier
710 }
711 
712 
713 //====================================================================
714 template<typename AFIELD>
715 void AFopr_Domainwall_5din<AFIELD>::mult_Ddag(AFIELD& v, const AFIELD& w)
716 {
717  real_t *vp = v.ptr(0);
718  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
719  real_t *yp = m_w1.ptr(0);
720  real_t *up = m_U.ptr(0);
721 
722  int ith = ThreadManager::get_thread_id();
723 
724  BridgeQXS::mult_domainwall_5din_mult_gm5_dirac(vp, wp, m_Ns, m_Nsize);
725 
726 #pragma omp barrier
727 
728  if (do_comm_any > 0) {
729  if (ith == 0) chset_recv.start();
730 
731  real_t *buf1_xp = (real_t *)chsend_dn[0].ptr();
732  real_t *buf1_xm = (real_t *)chsend_up[0].ptr();
733  real_t *buf1_yp = (real_t *)chsend_dn[1].ptr();
734  real_t *buf1_ym = (real_t *)chsend_up[1].ptr();
735  real_t *buf1_zp = (real_t *)chsend_dn[2].ptr();
736  real_t *buf1_zm = (real_t *)chsend_up[2].ptr();
737  real_t *buf1_tp = (real_t *)chsend_dn[3].ptr();
738  real_t *buf1_tm = (real_t *)chsend_up[3].ptr();
739 
740  BridgeQXS::mult_domainwall_5din_hop1_dirac(
741  buf1_xp, buf1_xm, buf1_yp, buf1_ym,
742  buf1_zp, buf1_zm, buf1_tp, buf1_tm,
743  up, vp,
744  m_mq, m_M0, m_Ns, &m_boundary[0],
745  m_Nsize, do_comm);
746 
747 #pragma omp barrier
748 
749  if (ith == 0) chset_send.start();
750  }
751 
752  BridgeQXS::mult_domainwall_5din_clear(yp, m_Ns, m_Nsize);
753 
754  BridgeQXS::mult_domainwall_5din_hopb_dirac(yp, up, vp,
755  m_mq, m_M0, m_Ns, &m_boundary[0],
756  &m_b[0], &m_c[0],
757  m_Nsize, do_comm);
758 
759  if (do_comm_any > 0) {
760  if (ith == 0) chset_recv.wait();
761 
762 #pragma omp barrier
763 
764  real_t *buf2_xp = (real_t *)chrecv_up[0].ptr();
765  real_t *buf2_xm = (real_t *)chrecv_dn[0].ptr();
766  real_t *buf2_yp = (real_t *)chrecv_up[1].ptr();
767  real_t *buf2_ym = (real_t *)chrecv_dn[1].ptr();
768  real_t *buf2_zp = (real_t *)chrecv_up[2].ptr();
769  real_t *buf2_zm = (real_t *)chrecv_dn[2].ptr();
770  real_t *buf2_tp = (real_t *)chrecv_up[3].ptr();
771  real_t *buf2_tm = (real_t *)chrecv_dn[3].ptr();
772 
773  BridgeQXS::mult_domainwall_5din_hop2_dirac(yp, up, vp,
774  buf2_xp, buf2_xm, buf2_yp, buf2_ym,
775  buf2_zp, buf2_zm, buf2_tp, buf2_tm,
776  m_mq, m_M0, m_Ns, &m_boundary[0],
777  m_Nsize, do_comm);
778 
779  if (ith == 0) chset_send.wait();
780  }
781 
782 #pragma omp barrier
783 
784  BridgeQXS::mult_domainwall_5din_5dirdag_dirac(vp, yp, wp,
785  m_mq, m_M0, m_Ns, &m_boundary[0],
786  &m_b[0], &m_c[0],
787  m_Nsize, do_comm);
788 
789 #pragma omp barrier
790 }
791 
792 
793 //====================================================================
794 template<typename AFIELD>
795 void AFopr_Domainwall_5din<AFIELD>::L_inv(AFIELD& v,
796  const AFIELD& w)
797 {
798  real_t *vp = v.ptr(0);
799  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
800  BridgeQXS::mult_domainwall_5din_L_inv_dirac(vp, wp,
801  m_Ns, m_Nsize,
802  &m_e[0], &m_dpinv[0], &m_dm[0]);
803 }
804 
805 
806 //====================================================================
807 template<typename AFIELD>
808 void AFopr_Domainwall_5din<AFIELD>::U_inv(AFIELD& v,
809  const AFIELD& w)
810 {
811  real_t *vp = v.ptr(0);
812  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
813  BridgeQXS::mult_domainwall_5din_U_inv_dirac(
814  vp, wp, m_Ns, m_Nsize,
815  &m_f[0], &m_dpinv[0], &m_dm[0]);
816 }
817 
818 
819 //====================================================================
820 template<typename AFIELD>
821 void AFopr_Domainwall_5din<AFIELD>::Ldag_inv(AFIELD& v,
822  const AFIELD& w)
823 {
824  real_t *vp = v.ptr(0);
825  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
826 
827  BridgeQXS::mult_domainwall_5din_Ldag_inv_dirac(
828  vp, wp, m_Ns, m_Nsize,
829  &m_e[0], &m_dpinv[0], &m_dm[0]);
830 }
831 
832 
833 //====================================================================
834 template<typename AFIELD>
835 void AFopr_Domainwall_5din<AFIELD>::Udag_inv(AFIELD& v,
836  const AFIELD& w)
837 {
838  real_t *vp = v.ptr(0);
839  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
840 
841  BridgeQXS::mult_domainwall_5din_Udag_inv_dirac(
842  vp, wp, m_Ns, m_Nsize,
843  &m_f[0], &m_dpinv[0], &m_dm[0]);
844 }
845 
846 
847 //====================================================================
848 template<typename AFIELD>
849 double AFopr_Domainwall_5din<AFIELD>::flop_count(std::string mode)
850 {
851  int Lvol = CommonParameters::Lvol();
852  double vsite = static_cast<double>(Lvol);
853  double vNs = static_cast<double>(m_Ns);
854  int Nc = CommonParameters::Nc();
855  int Nd = CommonParameters::Nd();
856 
857  double flop_Wilson;
858  double flop_LU_inv;
859  double axpy1 = static_cast<double>(2 * m_NinF);
860  double scal1 = static_cast<double>(1 * m_NinF);
861  if (m_repr == "Dirac") {
862  flop_Wilson = static_cast<double>(
863  Nc * Nd * (4 + 6 * (4 * Nc + 2) + 2 * (4 * Nc + 1))) * vsite;
864  flop_LU_inv = static_cast<double>(Nc * Nd * (2 + (vNs - 1) * 26)) * vsite;
865  } else if (m_repr == "Chiral") {
866  flop_Wilson = static_cast<double>(
867  Nc * Nd * (4 + 8 * (4 * Nc + 2))) * vsite;
868  flop_LU_inv = static_cast<double>(Nc * Nd * (2 + (vNs - 1) * 10)) * vsite;
869  }
870 
871  // Note that m_NinF := m_Nvcd * m_Ns;
872  double flop_DW = vNs * flop_Wilson + vsite * (6 * axpy1 + 2 * scal1);
873  // In Ddag case, flop_Wilson + 7 axpy which equals flop_DW.
874 
875  // double flop_LU_inv = 2.0 * vsite *
876  // ((3.0*axpy1 + scal1)*(vNs-1.0) + axpy1 + 2.0*scal1);
877 
878  double flop = 0.0;
879  if (mode == "Prec") {
880  flop = flop_LU_inv;
881  } else if ((mode == "D") || (mode == "Ddag")) {
882  flop = flop_DW;
883  } else if (mode == "DdagD") {
884  flop = 2.0 * flop_DW;
885  } else if ((mode == "D_prec") || (mode == "Ddag_prec")) {
886  flop = flop_LU_inv + flop_DW;
887  } else if (mode == "DdagD_prec") {
888  flop = 2.0 * (flop_LU_inv + flop_DW);
889  } else {
890  vout.crucial(m_vl, "Error at %s: input mode is undefined.\n",
891  class_name.c_str());
892  exit(EXIT_FAILURE);
893  }
894 
895  return flop;
896 }
897 
898 
899 //============================================================END=====
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
AFopr_Domainwall_5din::set_coefficients
void set_coefficients(const std::vector< double > b, const std::vector< double > c)
set coefficients if they depend in s.
Definition: afopr_Domainwall_5din-tmpl.h:301
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
AFopr_Domainwall_5din::Ldag_inv
void Ldag_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:821
BridgeQXS::mult_domainwall_5din_hopb_dirac
void mult_domainwall_5din_hopb_dirac(double *vp, double *up, double *wp, double mq, double M0, int Ns, int *bc, double *b, double *c, int *Nsize, int *do_comm)
Definition: mult_Domainwall_5din_qxs-inc.h:246
AFopr_Domainwall_5din::mult
void mult(AFIELD &v, const AFIELD &w)
multiplies fermion operator to a given field.
Definition: afopr_Domainwall_5din.h:126
AFopr_Domainwall_5din::Precdag
void Precdag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:609
BridgeQXS::mult_domainwall_5din_hop2_dirac
void mult_domainwall_5din_hop2_dirac(double *vp, double *up, double *wp, double *buf2_xp, double *buf2_xm, double *buf2_yp, double *buf2_ym, double *buf2_zp, double *buf2_zm, double *buf2_tp, double *buf2_tm, double mq, double M0, int Ns, int *bc, int *Nsize, int *do_comm)
Definition: mult_Domainwall_5din_qxs-inc.h:575
AFopr_Domainwall_5din::mult_Ddag
void mult_Ddag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:715
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
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
AFopr_Domainwall_5din::mult_dag
void mult_dag(AFIELD &v, const AFIELD &w)
hermitian conjugate of mult.
Definition: afopr_Domainwall_5din.h:129
Parameters
Class for parameters.
Definition: parameters.h:46
AIndex_lex
Definition: aindex_lex_base.h:17
AFopr_Domainwall_5din::Ddag
void Ddag(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:628
BridgeQXS::mult_domainwall_5din_Udag_inv_dirac
void mult_domainwall_5din_Udag_inv_dirac(real_t *vp, real_t *wp, int Ns, int *Nsize, real_t *f, real_t *dpinv, real_t *dm)
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:517
BridgeQXS::mult_domainwall_5din_5dir_dirac
void mult_domainwall_5din_5dir_dirac(double *vp, double *yp, double *wp, double mq, double M0, int Ns, int *bc, double *b, double *c, int *Nsize, int *do_comm)
Definition: mult_Domainwall_5din_qxs-inc.h:16
NVCD
#define NVCD
Definition: define_params_SU3.h:20
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
VLEN
#define VLEN
Definition: bridgeQXS_Clover_coarse_double.cpp:12
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:281
AFopr_Domainwall_5din::set_config_omp
void set_config_omp(Field *u)
setting gauge configuration (setting omp parallel).
Definition: afopr_Domainwall_5din-tmpl.h:356
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
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
convert_gauge
void convert_gauge(INDEX &index, FIELD &v, const Field &w)
Definition: afield-inc.h:224
BridgeQXS::mult_domainwall_5din_L_inv_dirac
void mult_domainwall_5din_L_inv_dirac(real_t *vp, real_t *wp, int Ns, int *Nsize, real_t *e, real_t *dpinv, real_t *dm)
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
AFopr_Domainwall_5din::D
void D(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:620
AFopr_Domainwall_5din::setup_channels
void setup_channels()
setup channels for communication.
Definition: afopr_Domainwall_5din-tmpl.h:126
AFopr_Domainwall_5din::D_prec
void D_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:575
AFopr_Domainwall_5din::convert
void convert(AFIELD &, const Field &)
convert Field to AField for this class.
Definition: afopr_Domainwall_5din-tmpl.h:384
AFopr_Domainwall_5din::mult_D
void mult_D(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:636
BridgeQXS::mult_domainwall_5din_clear
void mult_domainwall_5din_clear(double *vp, int Ns, int *Nsize)
Definition: mult_Domainwall_5din_qxs-inc.h:218
AFopr_Domainwall_5din::Ddag_prec
void Ddag_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:587
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
AFopr_Domainwall_5din::set_config_impl
void set_config_impl(Field *u)
setting gauge configuration (implementation).
Definition: afopr_Domainwall_5din-tmpl.h:369
AFopr_Domainwall_5din
Optimal Domain-wall fermion operator.
Definition: afopr_Domainwall_5din.h:35
AFopr_Domainwall_5din::set_precond_parameters
void set_precond_parameters()
set parameters for preconditioning.
Definition: afopr_Domainwall_5din-tmpl.h:261
AFopr_Domainwall_5din::L_inv
void L_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:795
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
NC
#define NC
Definition: field_F_imp_SU2-inc.h:15
BridgeQXS::mult_domainwall_5din_U_inv_dirac
void mult_domainwall_5din_U_inv_dirac(real_t *vp, real_t *wp, int Ns, int *Nsize, real_t *f, real_t *dpinv, real_t *dm)
QXS_Gauge::set_boundary
void set_boundary(AField< REALTYPE, QXS > &ulex, const std::vector< int > &boundary)
Definition: afield_Gauge-inc.h:24
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
Communicator::npe
static int npe(const int dir)
logical grid extent
Definition: communicator.cpp:112
AFopr_Domainwall_5din::Prec
void Prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:598
BridgeQXS::mult_domainwall_5din_hop1_dirac
void mult_domainwall_5din_hop1_dirac(double *buf1_xp, double *buf1_xm, double *buf1_yp, double *buf1_ym, double *buf1_zp, double *buf1_zm, double *buf1_tp, double *buf1_tm, double *up, double *wp, double mq, double M0, int Ns, int *bc, int *Nsize, int *do_comm)
Definition: mult_Domainwall_5din_qxs-inc.h:400
AFopr_Domainwall_5din::real_t
AFIELD::real_t real_t
Definition: afopr_Domainwall_5din.h:38
ND
#define ND
Definition: field_F_imp_SU2-inc.h:18
AFopr_Domainwall_5din::DdagD
void DdagD(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:556
AFopr_Domainwall_5din::DdagD_prec
void DdagD_prec(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:565
Field::cmp
double cmp(const int jin, const int site, const int jex) const
Definition: field.h:143
AFopr_Domainwall_5din::reverse
void reverse(Field &, const AFIELD &)
reverse AField to Field.
Definition: afopr_Domainwall_5din-tmpl.h:424
VLENY
#define VLENY
Definition: bridgeQXS_Clover_coarse_double.cpp:14
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
AFopr_Domainwall_5din::tidyup
void tidyup()
final tidyup.
Definition: afopr_Domainwall_5din-tmpl.h:118
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
AFopr_Domainwall_5din::U_inv
void U_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:808
BridgeQXS::mult_domainwall_5din_mult_gm5_dirac
void mult_domainwall_5din_mult_gm5_dirac(double *vp, double *wp, int Ns, int *Nsize)
Definition: mult_Domainwall_5din_qxs-inc.h:162
AFopr_Domainwall_5din::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-tmpl.h:464
AFopr_Domainwall_5din::mult_gm5
void mult_gm5(AFIELD &v, const AFIELD &w)
mult_dag with specified mode.
Definition: afopr_Domainwall_5din-tmpl.h:541
VLENX
#define VLENX
Definition: bridgeQXS_Clover_coarse_double.cpp:13
AFopr_Domainwall_5din::Udag_inv
void Udag_inv(AFIELD &, const AFIELD &)
Definition: afopr_Domainwall_5din-tmpl.h:835
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::flop_count
double flop_count()
this returns the number of floating point number operations.
Definition: afopr_Domainwall_5din.h:166
AFopr_Domainwall_5din::init
void init(const Parameters &params)
initial setup.
Definition: afopr_Domainwall_5din-tmpl.h:15
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
AFopr_Domainwall_5din::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Domainwall_5din-tmpl.h:160
BridgeQXS::mult_domainwall_5din_Ldag_inv_dirac
void mult_domainwall_5din_Ldag_inv_dirac(real_t *vp, real_t *wp, int Ns, int *Nsize, real_t *e, real_t *dpinv, real_t *dm)
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
ND2
#define ND2
Definition: define_params_SU3.h:18
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
BridgeQXS::mult_domainwall_5din_5dirdag_dirac
void mult_domainwall_5din_5dirdag_dirac(double *vp, double *yp, double *wp, double mq, double M0, int Ns, int *bc, double *b, double *c, int *Nsize, int *do_comm)
Definition: mult_Domainwall_5din_qxs-inc.h:79
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:569
AFopr_Domainwall_5din::set_config
void set_config(Field *U)
setting gauge configuration (common interface).
Definition: afopr_Domainwall_5din-tmpl.h:337