Bridge++  Ver. 2.0.2
afopr_Clover_eo-tmpl.h
Go to the documentation of this file.
1 
10 #include "lib/Field/index_eo.h"
11 #include "lib/ResourceManager/threadManager.h"
12 
13 template<typename AFIELD>
14 const std::string AFopr_Clover_eo<AFIELD>::class_name = "AFopr_Clover_eo";
15 
16 //====================================================================
17 template<typename AFIELD>
18 void AFopr_Clover_eo<AFIELD>::init(const Parameters& params)
19 {
20  ThreadManager::assert_single_thread(class_name);
21 
22  int req_comm = 1; // set 1 if communication forced any time
23  //int req_comm = 0;
24 
25  std::string vlevel;
26  if (!params.fetch_string("verbose_level", vlevel)) {
27  m_vl = vout.set_verbose_level(vlevel);
28  } else {
29  m_vl = CommonParameters::Vlevel();
30  }
31 
32  vout.general(m_vl, "%s: construction\n", class_name.c_str());
33 
34  m_repr = "Dirac"; // now only the Dirac repr is available.
35 
36  std::string repr;
37  if (!params.fetch_string("gamma_matrix_type", repr)) {
38  if (repr != "Dirac") {
39  vout.crucial(" Error at %s: unsupported gamma-matrix type: %s\n",
40  class_name.c_str(), repr.c_str());
41  exit(EXIT_FAILURE);
42  }
43  }
44 
45  m_Nc = CommonParameters::Nc();
46  if (m_Nc != 3) {
47  vout.crucial("%s: only applicable to Nc = 3\n",
48  class_name.c_str());
49  exit(EXIT_FAILURE);
50  }
51 
52  m_Nd = CommonParameters::Nd();
53  m_vl = CommonParameters::Vlevel();
54 
55  m_Ndf = 2 * m_Nc * m_Nc;
56  m_Nvc = m_Nc * 2;
57  m_Ndim = CommonParameters::Ndim();
58  m_Nx = CommonParameters::Nx();
59  m_Ny = CommonParameters::Ny();
60  m_Nz = CommonParameters::Nz();
61  m_Nt = CommonParameters::Nt();
62  m_Nst = CommonParameters::Nvol();
63 
64  m_Nx2 = m_Nx / 2;
65  m_Nst2 = m_Nst / 2;
66  m_Nx2v = m_Nx2 / VLENX;
67  m_Nyv = m_Ny / VLENY;
68  m_Nst2v = m_Nst2 / VLEN;
69 
70  // condition check
71  if (m_Nx % 2 != 0) {
72  vout.crucial(m_vl, "%s: Nx must be even.\n",
73  class_name.c_str());
74  exit(EXIT_FAILURE);
75  }
76  if (m_Nx % (2 * VLENX) != 0) {
77  vout.crucial(m_vl, "%s: Nx must be mulriple of 2*VLEN.\n",
78  class_name.c_str());
79  exit(EXIT_FAILURE);
80  }
81  if (m_Ny % (VLENY) != 0) {
82  vout.crucial(m_vl, "%s: Ny must be multiple of VLENY.\n",
83  class_name.c_str());
84  exit(EXIT_FAILURE);
85  }
86 
87  vout.general(m_vl, " VLENX = %2d Nx2v = %d\n", VLENX, m_Nx2v);
88  vout.general(m_vl, " VLENY = %2d Nyv = %d\n", VLENY, m_Nyv);
89  vout.general(m_vl, " VLEN = %2d Nst2v = %d\n", VLEN, m_Nst2v);
90 
91  m_Leo.resize(m_Ny * m_Nz * m_Nt);
92 
93  int ipe3 = Communicator::ipe(3);
94  int ipe2 = Communicator::ipe(2);
95  int ipe1 = Communicator::ipe(1);
96  for (int t = 0; t < m_Nt; ++t) {
97  for (int z = 0; z < m_Nz; ++z) {
98  for (int y = 0; y < m_Ny; ++y) {
99  int t2 = ipe3 * m_Nt + t;
100  int z2 = ipe2 * m_Nz + z;
101  int y2 = ipe1 * m_Ny + y;
102  m_Leo[y + m_Ny * (z + m_Nz * t)] = (y2 + z2 + t2) % 2;
103  }
104  }
105  }
106 
107  m_Nsize[0] = m_Nx2v;
108  m_Nsize[1] = m_Nyv;
109  m_Nsize[2] = m_Nz;
110  m_Nsize[3] = m_Nt;
111 
112  do_comm_any = 0;
113  for (int mu = 0; mu < m_Ndim; ++mu) {
114  do_comm[mu] = 1;
115  if ((req_comm == 0) && (Communicator::npe(mu) == 1)) do_comm[mu] = 0;
116  do_comm_any += do_comm[mu];
117  vout.general(" do_comm[%d] = %d\n", mu, do_comm[mu]);
118  }
119 
120  // setup of communication buffers
121  m_bdsize.resize(m_Ndim);
122  int Nd2 = m_Nd / 2;
123  m_bdsize[0] = m_Nvc * Nd2 * ((m_Ny * m_Nz * m_Nt + 1) / 2);
124  m_bdsize[1] = m_Nvc * Nd2 * m_Nx2 * m_Nz * m_Nt;
125  m_bdsize[2] = m_Nvc * Nd2 * m_Nx2 * m_Ny * m_Nt;
126  m_bdsize[3] = m_Nvc * Nd2 * m_Nx2 * m_Ny * m_Nz;
127 
128  setup_channels();
129 
130  // gauge configuration.
131  m_Ueo.reset(m_Ndf, m_Nst, m_Ndim);
132  m_Ulex.reset(m_Ndf, m_Nst, m_Ndim);
133 
134  int Ndm2 = m_Nd * m_Nd / 2;
135 
136  vout.increase_indent();
137 
138  Parameters params_ct = params;
139 
140 #ifdef CHIRAL_ROTATION
141  m_Te_inv.reset(NDF * Ndm2 / 2, m_Nst2, 1);
142  m_To_inv.reset(NDF * Ndm2 / 2, m_Nst2, 1);
143  params_ct.set_string("gamma_matrix_type", "Chiral");
144  m_fopr_ct = new Fopr_CloverTerm_eo(params_ct);
145 #else
146  m_Te_inv.reset(NDF * Ndm2, m_Nst2, 1);
147  m_To_inv.reset(NDF * Ndm2, m_Nst2, 1);
148  m_fopr_ct = new Fopr_CloverTerm_eo(params_ct);
149 #endif
150 
151  set_parameters(params);
152 
153  vout.decrease_indent();
154 
155  // working vectors.
156  int NinF = 2 * m_Nc * m_Nd;
157  m_v1.reset(NinF, m_Nst2, 1);
158  m_v2.reset(NinF, m_Nst2, 1);
159 
160  m_z1.reset(NinF, m_Nst, 1); // used in convert/revese.
161 
162  vout.general(m_vl, "%s: construction finished.\n",
163  class_name.c_str());
164 }
165 
166 
167 //====================================================================
168 template<typename AFIELD>
169 void AFopr_Clover_eo<AFIELD>::setup_channels()
170 {
171  chsend_up.resize(m_Ndim);
172  chrecv_up.resize(m_Ndim);
173  chsend_dn.resize(m_Ndim);
174  chrecv_dn.resize(m_Ndim);
175 
176  for (int mu = 0; mu < m_Ndim; ++mu) {
177  size_t Nvsize = m_bdsize[mu] * sizeof(real_t);
178 
179  chsend_dn[mu].send_init(Nvsize, mu, -1);
180  chsend_up[mu].send_init(Nvsize, mu, 1);
181 #ifdef USE_MPI
182  chrecv_up[mu].recv_init(Nvsize, mu, 1);
183  chrecv_dn[mu].recv_init(Nvsize, mu, -1);
184 #else
185  void *buf_up = (void *)chsend_dn[mu].ptr();
186  chrecv_up[mu].recv_init(Nvsize, mu, 1, buf_up);
187  void *buf_dn = (void *)chsend_up[mu].ptr();
188  chrecv_dn[mu].recv_init(Nvsize, mu, -1, buf_dn);
189 #endif
190 
191  if (do_comm[mu] == 1) {
192  chset_send.append(chsend_up[mu]);
193  chset_send.append(chsend_dn[mu]);
194  chset_recv.append(chrecv_up[mu]);
195  chset_recv.append(chrecv_dn[mu]);
196  }
197  }
198 }
199 
200 
201 //====================================================================
202 template<typename AFIELD>
203 void AFopr_Clover_eo<AFIELD>::tidyup()
204 {
205  ThreadManager::assert_single_thread(class_name);
206 
207  delete m_fopr_ct;
208 }
209 
210 
211 //====================================================================
212 template<typename AFIELD>
213 void AFopr_Clover_eo<AFIELD>::set_parameters(const Parameters& params)
214 {
215  const string str_vlevel = params.get_string("verbose_level");
216  m_vl = vout.set_verbose_level(str_vlevel);
217 
218  //- fetch and check input parameters
219  double kappa, csw;
220  std::vector<int> bc;
221 
222  int err = 0;
223  err += params.fetch_double("hopping_parameter", kappa);
224  err += params.fetch_double("clover_coefficient", csw);
225  err += params.fetch_int_vector("boundary_condition", bc);
226  if (err) {
227  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
228  class_name.c_str());
229  exit(EXIT_FAILURE);
230  }
231 
232  set_parameters(real_t(kappa), csw, bc);
233 
234  Parameters params_ct = params;
235 #ifdef CHIRAL_ROTATION
236  params_ct.set_string("gamma_matrix_type", "Chiral");
237 #endif
238  m_fopr_ct->set_parameters(params_ct);
239 }
240 
241 
242 //====================================================================
243 template<typename AFIELD>
244 void AFopr_Clover_eo<AFIELD>::set_parameters(const real_t CKs,
245  const real_t csw,
246  const std::vector<int> bc)
247 {
248  assert(bc.size() == m_Ndim);
249 
250 #pragma omp barrier
251 
252  int ith = ThreadManager::get_thread_id();
253 
254  if (ith == 0) {
255  m_CKs = CKs;
256  m_csw = csw;
257  m_boundary.resize(m_Ndim);
258  for (int mu = 0; mu < m_Ndim; ++mu) {
259  m_boundary[mu] = bc[mu];
260  }
261  }
262 
263  vout.general(m_vl, "%s: set parameters\n", class_name.c_str());
264  vout.general(m_vl, " gamma-matrix type = %s\n", m_repr.c_str());
265  vout.general(m_vl, " kappa = %8.4f\n", m_CKs);
266  vout.general(m_vl, " cSW = %12.8f\n", m_csw);
267  for (int mu = 0; mu < m_Ndim; ++mu) {
268  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
269  }
270 
271 #pragma omp barrier
272 }
273 
274 
275 //====================================================================
276 template<typename AFIELD>
277 void AFopr_Clover_eo<AFIELD>::get_parameters(Parameters& params) const
278 {
279  params.set_double("hopping_parameter", double(m_CKs));
280  params.set_double("clover_coefficient", double(m_csw));
281  params.set_int_vector("boundary_condition", m_boundary);
282  params.set_string("gamma_matrix_type", m_repr);
283 
284  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
285 }
286 
287 
288 //====================================================================
289 template<typename AFIELD>
290 void AFopr_Clover_eo<AFIELD>::set_config(Field *u)
291 {
292  int nth = ThreadManager::get_num_threads();
293 
294  vout.detailed(m_vl, "%s: set_config is called: num_threads = %d\n",
295  class_name.c_str(), nth);
296 
297  if (nth > 1) {
298  set_config_impl(u);
299  } else {
300  set_config_omp(u);
301  }
302 
303  vout.detailed(m_vl, "%s: set_config finished\n", class_name.c_str());
304 }
305 
306 
307 //====================================================================
308 template<typename AFIELD>
309 void AFopr_Clover_eo<AFIELD>::set_config_omp(Field *u)
310 {
311  vout.detailed(m_vl, " set_config_omp is called.\n");
312 
313 #pragma omp parallel
314  {
315  set_config_impl(u);
316  }
317 }
318 
319 
320 //====================================================================
321 template<typename AFIELD>
322 void AFopr_Clover_eo<AFIELD>::set_config_impl(Field *u)
323 {
324 #pragma omp barrier
325 
326  int ith, nth, is, ns;
327  ith = ThreadManager::get_thread_id();
328 
329  if (ith == 0) m_conf = u;
330 
331  AIndex_lex<real_t, QXS> index_lex;
332 
333  convert_gauge(index_lex, m_Ulex, *u);
334 
335  QXS_Gauge::set_boundary(m_Ulex, m_boundary);
336 
337  AIndex_lex<real_t, QXS> index_lex2;
338  AIndex_eo<real_t, QXS> index_eo;
339 
340  set_threadtask_mult(ith, nth, is, ns, m_Nst);
341 
342  for (int ex = 0; ex < m_Ndim; ++ex) {
343  for (int site = is; site < ns; ++site) {
344  for (int in = 0; in < m_Ndf; ++in) {
345  int iv1 = index_lex2.idx(in, m_Ndf, site, ex);
346  int iv2 = index_eo.idx(in, m_Ndf, site, ex);
347  m_Ueo.e(iv2) = m_Ulex.cmp(iv1);
348  }
349  }
350  }
351 
352  // setup of clover term inverse
353  m_fopr_ct->set_config(u);
354 
355 #ifdef CHIRAL_ROTATION
356  set_csw_chrot(m_Te_inv, 0);
357  set_csw_chrot(m_To_inv, 1);
358 #else
359  set_csw();
360 #endif
361 
362 #pragma omp barrier
363 }
364 
365 
366 //====================================================================
367 template<typename AFIELD>
368 void AFopr_Clover_eo<AFIELD>::set_csw()
369 {
370  vout.crucial(m_vl, "%s: set_csw() has not been implemented.\n",
371  class_name.c_str());
372  exit(EXIT_FAILURE);
373 }
374 
375 
376 //====================================================================
377 template<typename AFIELD>
378 void AFopr_Clover_eo<AFIELD>::set_csw_chrot(AFIELD& m_T_inv, int ieo)
379 {
380  // This method set the clover term matrix assuming the chiral
381  // representation for gamma-matrix in Bridge++ convention.
382  // The conversion from the Bridge++ to QWS convention and
383  // from Dirac to chiral representations are assumed to be cared
384  // in BridgeQXS functions.
385  // [04 Sep 2021 H.Matsufuru]
386 
387  // Field_F v(m_Nst2, 1), w(m_Nst2, 1);
388 
389  AIndex_eo<real_t, AFIELD::IMPL> index;
390 
391  // const int Nin = NDF * ND * 2;
392  const int Nin = NDF * ND;
393 
394  m_fopr_ct->set_mode("D");
395 
396  int ith, nth, is, ns;
397  set_threadtask_mult(ith, nth, is, ns, m_Nst2);
398 
399 //11 22 33 44 55 66 12 34 56 23 45 13 24 35 46 15 26 14 36 25 16
400  constexpr int idx[72] = {
401  0, -1, 6, 7, 16, 17, 28, 29, 24, 25, 34, 35,
402  -1, -1, 1, -1, 12, 13, 18, 19, 32, 33, 26, 27,
403  -1, -1, -1, -5, 2, -1, 8, 9, 20, 21, 30, 31,
404  -1, -1, -1, -1, -1, -1, 3, -1, 14, 15, 22, 23,
405  -1, -1, -1, -1, -1, -1, -1, -1, 4, -1, 10, 11,
406  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 5, -1,
407  };
408 //#pragma omp parallel
409 // {
410  for (int id = 0; id < m_Nd / 2; ++id) {
411  for (int ic = 0; ic < m_Nc; ++ic) {
412  m_w1.set(0.0);
413 #pragma omp barrier
414 
415  for (int site = is; site < ns; ++site) {
416  m_w1.set_r(ic, id, site, 0, 1.0);
417  m_w1.set_r(ic, id + 2, site, 0, 1.0);
418  }
419 #pragma omp barrier
420 
421  m_fopr_ct->mult_csw_inv(m_w2, m_w1, ieo);
422 #pragma omp barrier
423 
424  for (int site = is; site < ns; ++site) {
425  for (int id2 = 0; id2 < m_Nd; ++id2) {
426  for (int ic2 = 0; ic2 < m_Nc; ++ic2) {
427  real_t vt_r = 0.5 * m_w2.cmp_r(ic2, id2, site, 0);
428  real_t vt_i = 0.5 * m_w2.cmp_i(ic2, id2, site, 0);
429  int i = ic2 + m_Nc * (id2 % 2);
430  int j = ic + m_Nc * id;
431  int ij = m_Nc * 2 * i + j;
432  int in_r = idx[2 * ij];
433  int in_i = idx[2 * ij + 1];
434  if (in_r >= 0) {
435  in_r += 36 * (id2 / 2);
436  int idx_r = index.idxh(in_r, Nin, site, 0);
437  m_T_inv.set(idx_r, vt_r);
438  }
439  if (in_i >= 0) {
440  in_i += 36 * (id2 / 2);
441  int idx_i = index.idxh(in_i, Nin, site, 0);
442  m_T_inv.set(idx_i, vt_i);
443  }
444  }
445  }
446  } // site loop
447 #pragma omp barrier
448  }
449  }
450 #pragma omp barrier
451 
452  // } // omp parallel
453 }
454 
455 
456 //====================================================================
457 template<typename AFIELD>
458 void AFopr_Clover_eo<AFIELD>::convert(AFIELD& v, const Field& w)
459 {
460  AIndex_lex<real_t, AFIELD::IMPL> index_lex;
461  convert_spinor(index_lex, m_z1, w);
462 
463  mult_gm4(v, m_z1);
464 }
465 
466 
467 //====================================================================
468 template<typename AFIELD>
469 void AFopr_Clover_eo<AFIELD>::reverse(Field& v, const AFIELD& w)
470 {
471  mult_gm4(m_z1, w);
472 
473  AIndex_lex<real_t, AFIELD::IMPL> index_lex;
474  reverse_spinor(index_lex, v, m_z1);
475 }
476 
477 
478 //====================================================================
479 template<typename AFIELD>
480 void AFopr_Clover_eo<AFIELD>::set_mode(std::string mode)
481 {
482 #pragma omp barrier
483 
484  int ith = ThreadManager::get_thread_id();
485  if (ith == 0) m_mode = mode;
486 
487 #pragma omp barrier
488 }
489 
490 
491 //====================================================================
492 template<typename AFIELD>
493 void AFopr_Clover_eo<AFIELD>::mult(AFIELD& v, const AFIELD& w)
494 {
495  if (m_mode == "D") {
496  return D(v, w);
497  } else if (m_mode == "DdagD") {
498  return DdagD(v, w);
499  } else if (m_mode == "Ddag") {
500  return Ddag(v, w);
501  } else {
502  vout.crucial(m_vl, "%s: mode undefined.\n", class_name.c_str());
503  exit(EXIT_FAILURE);
504  }
505 }
506 
507 
508 //====================================================================
509 template<typename AFIELD>
510 void AFopr_Clover_eo<AFIELD>::mult_dag(AFIELD& v, const AFIELD& w)
511 {
512  if (m_mode == "D") {
513  return Ddag(v, w);
514  } else if (m_mode == "DdagD") {
515  return DdagD(v, w);
516  } else if (m_mode == "Ddag") {
517  return D(v, w);
518  } else {
519  vout.crucial(m_vl, "%s: mode undefined.\n", class_name.c_str());
520  exit(EXIT_FAILURE);
521  }
522 }
523 
524 
525 //====================================================================
526 template<typename AFIELD>
527 void AFopr_Clover_eo<AFIELD>::mult_gm5(AFIELD& v, const AFIELD& w)
528 {
529  real_t *vp = v.ptr(0);
530  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
531 
532 #pragma omp barrier
533 
534  BridgeQXS::mult_wilson_gm5_dirac(vp, wp, m_Nsize);
535 
536 #pragma omp barrier
537 
538  //#pragma omp barrier
539  // mult_gm5(vp, wp);
540 }
541 
542 
543 //====================================================================
544 template<typename AFIELD>
545 void AFopr_Clover_eo<AFIELD>::mult_gm4(AFIELD& v, const AFIELD& w)
546 {
547  int Nst = w.nvol();
548  if (Nst != v.nvol()) {
549  vout.crucial(m_vl, "%s: sizes of fields inconsisutent in mult_gm4.\n",
550  class_name.c_str());
551  exit(EXIT_FAILURE);
552  }
553 
554  real_t *vp = v.ptr(0);
555  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
556 
557  int ith, nth, is, ns;
558  set_threadtask_mult(ith, nth, is, ns, m_Nst);
559 
560 #pragma omp barrier
561 
562  // Dirac representation.
563  for (int site = is; site < ns; ++site) {
564  for (int ic = 0; ic < NC; ++ic) {
565  vp[m_index.idxh_SPr(ic, 0, site, 0)] = wp[m_index.idxh_SPr(ic, 0, site, 0)];
566  vp[m_index.idxh_SPi(ic, 0, site, 0)] = wp[m_index.idxh_SPi(ic, 0, site, 0)];
567  vp[m_index.idxh_SPr(ic, 1, site, 0)] = wp[m_index.idxh_SPr(ic, 1, site, 0)];
568  vp[m_index.idxh_SPi(ic, 1, site, 0)] = wp[m_index.idxh_SPi(ic, 1, site, 0)];
569  vp[m_index.idxh_SPr(ic, 2, site, 0)] = -wp[m_index.idxh_SPr(ic, 2, site, 0)];
570  vp[m_index.idxh_SPi(ic, 2, site, 0)] = -wp[m_index.idxh_SPi(ic, 2, site, 0)];
571  vp[m_index.idxh_SPr(ic, 3, site, 0)] = -wp[m_index.idxh_SPr(ic, 3, site, 0)];
572  vp[m_index.idxh_SPi(ic, 3, site, 0)] = -wp[m_index.idxh_SPi(ic, 3, site, 0)];
573  }
574  }
575 
576 #pragma omp barrier
577 }
578 
579 
580 //====================================================================
581 template<typename AFIELD>
582 void AFopr_Clover_eo<AFIELD>::mult(AFIELD& v, const AFIELD& w,
583  const std::string mode)
584 {
585  if (mode == "Dee_inv") {
586  mult_cswinv(v, w, 0);
587  } else if (mode == "Doo_inv") {
588  mult_cswinv(v, w, 1);
589  } else if (mode == "Deo") {
590  Meo(v, w, 0);
591  } else if (mode == "Doe") {
592  Meo(v, w, 1);
593  } else {
594  vout.crucial(m_vl, "%s: illegal mode is given to mult with mode\n",
595  class_name.c_str());
596  exit(EXIT_FAILURE);
597  }
598 }
599 
600 
601 //====================================================================
602 template<typename AFIELD>
603 void AFopr_Clover_eo<AFIELD>::DdagD(AFIELD& v, const AFIELD& w)
604 {
605  D(m_v2, w);
606  Ddag(v, m_v2);
607 }
608 
609 
610 //====================================================================
611 template<typename AFIELD>
612 void AFopr_Clover_eo<AFIELD>::Ddag(AFIELD& v, const AFIELD& w)
613 {
614 #pragma omp barrier
615 
616  mult_gm5(m_v1, w);
617  mult_cswinv(v, m_v1, 0);
618  Meo(m_v1, v, 1);
619  mult_cswinv(v, m_v1, 1);
620  Meo(m_v1, v, 0);
621  mult_gm5(v, m_v1);
622  aypx(-1.0, v, w);
623 }
624 
625 
626 //====================================================================
627 template<typename AFIELD>
628 void AFopr_Clover_eo<AFIELD>::D(AFIELD& v, const AFIELD& w)
629 {
630 #pragma omp barrier
631  Meo(m_v1, w, 1);
632  mult_cswinv(v, m_v1, 1);
633  Meo(m_v1, v, 0);
634  mult_cswinv(v, m_v1, 0);
635  aypx(-1.0, v, w);
636 }
637 
638 
639 //====================================================================
640 template<typename AFIELD>
641 void AFopr_Clover_eo<AFIELD>::aypx(real_t a, AFIELD& v, const AFIELD& w)
642 {
643  real_t *vp = v.ptr(0);
644  real_t *wp = const_cast<AFIELD *>(&w)->ptr(0);
645 
646 #pragma omp barrier
647 
648  aypx(a, vp, wp);
649 }
650 
651 
652 //====================================================================
653 template<typename AFIELD>
654 void AFopr_Clover_eo<AFIELD>::Meo(AFIELD& v, const AFIELD& w,
655  const int ieo)
656 {
657  mult_Meo_qxs(v, w, w, ieo, 0);
658 }
659 
660 
661 //====================================================================
662 template<typename AFIELD>
663 void AFopr_Clover_eo<AFIELD>::Meo(AFIELD& v,
664  const AFIELD& w, const AFIELD& x,
665  const int ieo, const int iflag)
666 {
667  mult_Meo_qxs(v, w, x, ieo, iflag);
668 }
669 
670 
671 //====================================================================
672 template<typename AFIELD>
673 void AFopr_Clover_eo<AFIELD>::mult_Meo_qxs(AFIELD& v,
674  const AFIELD& w, const AFIELD& x,
675  const int ieo, const int iflag)
676 {
677  real_t *v2 = v.ptr(0);
678  real_t *v1 = const_cast<AFIELD *>(&w)->ptr(0);
679  real_t *xp = const_cast<AFIELD *>(&x)->ptr(0);
680 
681  real_t *up = m_Ueo.ptr(0);
682 
683  int ith = ThreadManager::get_thread_id();
684 
685 #pragma omp barrier
686 
687  if (do_comm_any > 0) {
688  if (ith == 0) chset_recv.start();
689 
690  real_t *buf1_xp = (real_t *)chsend_dn[0].ptr();
691  real_t *buf1_xm = (real_t *)chsend_up[0].ptr();
692  real_t *buf1_yp = (real_t *)chsend_dn[1].ptr();
693  real_t *buf1_ym = (real_t *)chsend_up[1].ptr();
694  real_t *buf1_zp = (real_t *)chsend_dn[2].ptr();
695  real_t *buf1_zm = (real_t *)chsend_up[2].ptr();
696  real_t *buf1_tp = (real_t *)chsend_dn[3].ptr();
697  real_t *buf1_tm = (real_t *)chsend_up[3].ptr();
698  BridgeQXS::mult_wilson_eo_1_dirac(buf1_xp, buf1_xm, buf1_yp, buf1_ym,
699  buf1_zp, buf1_zm, buf1_tp, buf1_tm,
700  up, v1, &m_boundary[0],
701  m_Nsize, do_comm, &m_Leo[0], ieo, iflag);
702 
703 #pragma omp barrier
704 
705  if (ith == 0) chset_send.start();
706  }
707 
708  BridgeQXS::mult_wilson_eo_bulk_dirac(v2, up, v1, xp, m_CKs, &m_boundary[0],
709  m_Nsize, do_comm, &m_Leo[0], ieo, iflag);
710 
711  if (do_comm_any > 0) {
712  if (ith == 0) chset_recv.wait();
713 
714 #pragma omp barrier
715 
716  real_t *buf2_xp = (real_t *)chrecv_up[0].ptr();
717  real_t *buf2_xm = (real_t *)chrecv_dn[0].ptr();
718  real_t *buf2_yp = (real_t *)chrecv_up[1].ptr();
719  real_t *buf2_ym = (real_t *)chrecv_dn[1].ptr();
720  real_t *buf2_zp = (real_t *)chrecv_up[2].ptr();
721  real_t *buf2_zm = (real_t *)chrecv_dn[2].ptr();
722  real_t *buf2_tp = (real_t *)chrecv_up[3].ptr();
723  real_t *buf2_tm = (real_t *)chrecv_dn[3].ptr();
724 
725  BridgeQXS::mult_wilson_eo_2_dirac(v2, up, v1, xp,
726  buf2_xp, buf2_xm, buf2_yp, buf2_ym,
727  buf2_zp, buf2_zm, buf2_tp, buf2_tm,
728  m_CKs, &m_boundary[0],
729  m_Nsize, do_comm, &m_Leo[0], ieo, iflag);
730 
731  if (ith == 0) chset_send.wait();
732  }
733 
734 #pragma omp barrier
735 }
736 
737 
738 //====================================================================
739 template<typename AFIELD>
740 void AFopr_Clover_eo<AFIELD>::mult_cswinv(AFIELD& v, const AFIELD& w,
741  const int ieo)
742 {
743  real_t *v2 = v.ptr(0);
744  real_t *v1 = const_cast<AFIELD *>(&w)->ptr(0);
745  real_t *ct;
746  if (ieo == 0) {
747  ct = m_Te_inv.ptr(0);
748  } else {
749  ct = m_To_inv.ptr(0);
750  }
751 
752 #pragma omp barrier
753 
754 #ifdef CHIRAL_ROTATION
755  BridgeQXS::mult_clover_cswinv_dirac_chrot(v2, ct, v1, m_Nsize);
756 #endif
757 
758 #pragma omp barrier
759 }
760 
761 
762 //====================================================================
763 template<typename AFIELD>
764 void AFopr_Clover_eo<AFIELD>::aypx(real_t a, real_t *__restrict v,
765  real_t *__restrict w)
766 {
767  int Ncd = m_Nvc * m_Nd;
768 #ifdef USE_QXS_ACLE
769  svbool_t pg = set_predicate();
770 
771  int ith, nth, is, ns;
772  set_threadtask(ith, nth, is, ns, m_Nst2v);
773 
774  for (int site = is; site < ns; ++site) {
775  real_t *__restrict v2 = v + VLEN * Ncd * site;
776  real_t *__restrict w2 = w + VLEN * Ncd * site;
777  for (int icd = 0; icd < Ncd / 2; ++icd) {
778  svreal_t vtr, vti;
779  svreal_t wtr, wti;
780  load_vec(pg, wtr, &w2[2 * VLEN * icd]);
781  load_vec(pg, wti, &w2[2 * VLEN * icd + VLEN]);
782  load_vec(pg, vtr, &v2[2 * VLEN * icd]);
783  load_vec(pg, vti, &v2[2 * VLEN * icd + VLEN]);
784  aypx_vec(pg, a, vtr, wtr);
785  aypx_vec(pg, a, vti, wti);
786  save_vec(pg, &v2[2 * VLEN * icd], vtr);
787  save_vec(pg, &v2[2 * VLEN * icd + VLEN], vti);
788  }
789  }
790 #else
791  int ith, nth, is, ns;
792  set_threadtask_mult(ith, nth, is, ns, m_Nst2);
793  for (int i = is; i < ns; ++i) {
794  for (int icd = 0; icd < Ncd; ++icd) {
795  int i2 = icd + Ncd * i;
796  v[i2] = a * v[i2] + w[i2];
797  }
798  }
799 #endif
800 
801 #pragma omp barrier
802 }
803 
804 
805 //====================================================================
806 template<typename AFIELD>
807 void AFopr_Clover_eo<AFIELD>::clear(real_t *v)
808 {
809  int ith, nth, is, ns;
810  set_threadtask_mult(ith, nth, is, ns, m_Nst2);
811 
812  real_t zero = 0.0;
813  int Ncd = m_Nvc * m_Nd;
814 
815  for (int i = is; i < ns; ++i) {
816  for (int icd = 0; icd < Ncd; ++icd) {
817  v[icd + Ncd * i] = zero;
818  }
819  }
820 
821 #pragma omp barrier
822 }
823 
824 
825 //====================================================================
826 template<typename AFIELD>
827 void AFopr_Clover_eo<AFIELD>::scal(real_t *v, const real_t a)
828 {
829  int ith, nth, is, ns;
830  set_threadtask_mult(ith, nth, is, ns, m_Nst2);
831 
832  int Ncd = m_Nvc * m_Nd;
833 
834  for (int i = is; i < ns; ++i) {
835  for (int icd = 0; icd < Ncd; ++icd) {
836  v[icd + Ncd * i] *= a;
837  }
838  }
839 
840 #pragma omp barrier
841 }
842 
843 
844 //====================================================================
845 template<typename AFIELD>
846 double AFopr_Clover_eo<AFIELD>::flop_count(const std::string mode)
847 {
848  // The following counting explicitly depends on the implementation.
849  // It will be recalculated when the code is modified.
850  // The present counting is copied from AFopr_Clover in QXS impl.
851 
852  size_t Lvol = CommonParameters::Lvol();
853  double flop_wilson, flop_clover, flop_site, flop;
854 
855  if (m_repr == "Dirac") {
856  flop_wilson = static_cast<double>(
857  m_Nc * m_Nd * (4 // aypx
858  + 6 * (4 * m_Nc + 2) // spatial hopping
859  + 2 * (4 * m_Nc + 1))); // temporal hopping
860 
861  // clover term mult assumes rotation to chiral repr.
862  flop_clover = static_cast<double>(
863  m_Nc * m_Nd * (2 // Dirac -> chiral
864  + 2 * (2 * (m_Nc * m_Nd - 1) + 1) // clover term mult
865  + 2 // chiral -> Dirac
866  + 2)); // addition to vector
867  } else if (m_repr == "Chiral") {
868  flop_wilson = static_cast<double>(
869  m_Nc * m_Nd * (4 + 8 * (4 * m_Nc + 2)));
870 
871  flop_clover = static_cast<double>(
872  m_Nc * m_Nd * (2 * (2 * (m_Nc * m_Nd - 1) + 1) // clover term mult
873  + 2)); // addition to vector
874  } else {
875  vout.crucial(m_vl, "%s: input repr is undefined.\n");
876  exit(EXIT_FAILURE);
877  }
878 
879  flop_site = flop_wilson + flop_clover;
880 
881  flop = flop_site * static_cast<double>(Lvol);
882  if ((mode == "DdagD") || (mode == "DDdag")) flop *= 2.0;
883 
884  return flop;
885 }
886 
887 
888 //============================================================END=====
AFopr_Clover_eo::mult
void mult(AFIELD &, const AFIELD &)
multiplies fermion operator to a given field.
Definition: afopr_Clover_eo-tmpl.h:493
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
BridgeQXS::mult_wilson_gm5_dirac
void mult_wilson_gm5_dirac(double *v2, double *v1, int *Nsize)
Definition: mult_Wilson_qxs-inc.h:411
AFopr_Clover_eo::set_csw_chrot
void set_csw_chrot(AFIELD &m_T_inv, int ieo)
setting clover term inverse for chiral rotation case.
Definition: afopr_Clover_eo-tmpl.h:378
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
AFopr_Clover_eo::D
void D(AFIELD &, const AFIELD &)
Definition: afopr_Clover_eo-tmpl.h:628
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
Parameters
Class for parameters.
Definition: parameters.h:46
AIndex_lex
Definition: aindex_lex_base.h:17
AFopr_Clover_eo::convert
void convert(AFIELD &, const Field &)
convert Field to AField for this class.
Definition: afopr_Clover_eo-tmpl.h:458
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.h:86
AFopr_Clover_eo::set_config_impl
void set_config_impl(Field *u)
setting gauge configuration (implementation).
Definition: afopr_Clover_eo-tmpl.h:322
AFopr_Clover_eo::Meo
void Meo(AFIELD &, const AFIELD &, const int ieo)
Fermion matrix with ieo = 0: even <– odd, 1: odd <– even.
Definition: afopr_Clover_eo-tmpl.h:654
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
VLEN
#define VLEN
Definition: bridgeQXS_Clover_coarse_double.cpp:12
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
AFopr_Clover_eo::get_parameters
void get_parameters(Parameters &params) const
get parameters via a Parameter object
Definition: afopr_Clover_eo-tmpl.h:277
NDF
#define NDF
Definition: field_F_imp_SU2-inc.h:4
AFopr_Clover_eo::mult_cswinv
void mult_cswinv(AFIELD &, const AFIELD &, const int ieo)
Definition: afopr_Clover_eo-tmpl.h:740
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
AFopr_Clover_eo::mult_gm4
void mult_gm4(AFIELD &, const AFIELD &)
Definition: afopr_Clover_eo-tmpl.h:545
convert_gauge
void convert_gauge(INDEX &index, FIELD &v, const Field &w)
Definition: afield-inc.h:224
aypx
void aypx(const double a, Field &y, const Field &x)
aypx(y, a, x): y := a * y + x
Definition: field.cpp:509
Vsimd_t
Definition: vsimd_double-inc.h:13
AFopr_Clover_eo::real_t
AFIELD::real_t real_t
Definition: afopr_Clover_eo.h:43
AFopr_Clover_eo::clear
void clear(real_t *)
Definition: afopr_Clover_eo-tmpl.h:807
AFopr_Clover_eo::flop_count
double flop_count()
returns floating operation counts.
Definition: afopr_Clover_eo.h:165
AFopr_Clover_eo::mult_gm5
void mult_gm5(AFIELD &, const AFIELD &)
multiplies gamma_5 matrix.
Definition: afopr_Clover_eo-tmpl.h:527
AFopr_Clover_eo::setup_channels
void setup_channels()
setup channels for communication.
Definition: afopr_Clover_eo-tmpl.h:169
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
AFopr_Clover_eo::set_csw
void set_csw()
setting clover term inverse.
Definition: afopr_Clover_eo-tmpl.h:368
AFopr_Clover_eo::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
Definition: afopr_Clover_eo-tmpl.h:213
AFopr_Clover_eo::set_mode
void set_mode(std::string mode)
setting mult mode.
Definition: afopr_Clover_eo-tmpl.h:480
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
AFopr_Clover_eo::aypx
void aypx(real_t, AFIELD &, const AFIELD &)
Definition: afopr_Clover_eo-tmpl.h:641
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
NC
#define NC
Definition: field_F_imp_SU2-inc.h:2
AFopr_Clover_eo::scal
void scal(real_t *, const real_t)
Definition: afopr_Clover_eo-tmpl.h:827
QXS_Gauge::set_boundary
void set_boundary(AField< REALTYPE, QXS > &ulex, const std::vector< int > &boundary)
Definition: afield_Gauge-inc.h:24
reverse_spinor
void reverse_spinor(INDEX &index, Field &v, FIELD &w)
Definition: afield-inc.h:380
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_Clover_eo::set_config
void set_config(Field *u)
setting gauge configuration (common interface).
Definition: afopr_Clover_eo-tmpl.h:290
BridgeQXS::mult_wilson_eo_bulk_dirac
void mult_wilson_eo_bulk_dirac(double *v2, double *up, double *v1, double *xp, double kappa, int *bc, int *Nsize, int *do_comm, int *Leo, const int ieo, const int iflag)
Definition: mult_Wilson_eo_qxs-inc.h:20
AIndex_eo_qxs::idx
int idx(const int in, const int Nin, const int ist, const int Nx2, const int Ny, const int leo, const int Nvol2, const int ex)
Definition: aindex_eo.h:27
threadManager.h
Fopr_CloverTerm_eo
Org::Fopr_CloverTerm_eo Fopr_CloverTerm_eo
Definition: fopr_CloverTerm_eo.h:32
ND
#define ND
Definition: field_F_imp_SU2-inc.h:5
Field::nvol
int nvol() const
Definition: field.h:127
AFopr_Clover_eo::mult_dag
void mult_dag(AFIELD &, const AFIELD &)
hermitian conjugate of mult.
Definition: afopr_Clover_eo-tmpl.h:510
AFopr_Clover_eo::Ddag
void Ddag(AFIELD &, const AFIELD &)
Definition: afopr_Clover_eo-tmpl.h:612
AFopr_Clover_eo::mult_Meo_qxs
void mult_Meo_qxs(AFIELD &, const AFIELD &, const AFIELD &, const int ieo, const int iflag)
Meo implementation: using qxs library.
Definition: afopr_Clover_eo-tmpl.h:673
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
BridgeQXS::mult_wilson_eo_2_dirac
void mult_wilson_eo_2_dirac(double *v2, double *up, double *v1, double *xp, double *buf_xp, double *buf_xm, double *buf_yp, double *buf_ym, double *buf_zp, double *buf_zm, double *buf_tp, double *buf_tm, double kappa, int *bc, int *Nsize, int *do_comm, int *Leo, const int ieo, const int iflag)
Definition: mult_Wilson_eo_qxs-inc.h:429
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_Clover_eo::set_config_omp
void set_config_omp(Field *u)
setting gauge configuration (setting omp parallel).
Definition: afopr_Clover_eo-tmpl.h:309
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:133
svbool_t
Definition: vsimd_double-inc.h:30
AFopr_Clover_eo::init
void init(const Parameters &params)
initial setup.
Definition: afopr_Clover_eo-tmpl.h:18
VLENX
#define VLENX
Definition: bridgeQXS_Clover_coarse_double.cpp:13
AFopr_Clover_eo::reverse
void reverse(Field &, const AFIELD &)
reverse AField to Field.
Definition: afopr_Clover_eo-tmpl.h:469
AFopr_Clover_eo::DdagD
void DdagD(AFIELD &, const AFIELD &)
Definition: afopr_Clover_eo-tmpl.h:603
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_Clover_eo
Definition: afopr_Clover_eo.h:40
index_eo.h
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
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
BridgeQXS::mult_clover_cswinv_dirac_chrot
void mult_clover_cswinv_dirac_chrot(double *v2, double *ct, double *v1, int *Nsize)
Definition: mult_Clover_qxs-inc.h:288
AFopr_Clover_eo::tidyup
void tidyup()
final tidy-up.
Definition: afopr_Clover_eo-tmpl.h:203
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
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154
AIndex_eo
Definition: aindex_eo_base.h:17
convert_spinor
void convert_spinor(INDEX &index, FIELD &v, const Field &w)
Definition: afield-inc.h:187
BridgeQXS::mult_wilson_eo_1_dirac
void mult_wilson_eo_1_dirac(double *buf_xp, double *buf_xm, double *buf_yp, double *buf_ym, double *buf_zp, double *buf_zm, double *buf_tp, double *buf_tm, double *up, double *v1, int *bc, int *Nsize, int *do_comm, int *Leo, const int ieo, const int iflag)
Definition: mult_Wilson_eo_qxs-inc.h:238