Bridge++  Ver.2.1.3
fopr_Wilson_impl.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Wilson_impl.h"
15 
16 #include "bridge_defs.h"
17 #include "Fopr/fopr_thread-inc.h"
18 
19 #if defined USE_GROUP_SU3
20 #include "fopr_Wilson_impl_SU3-inc.h"
21 #elif defined USE_GROUP_SU2
22 #include "fopr_Wilson_impl_SU2-inc.h"
23 #elif defined USE_GROUP_SU_N
24 #include "fopr_Wilson_impl_SU_N-inc.h"
25 #endif
26 
27 #include "Fopr/Imp/fopr_Wilson_impl_common-inc.h"
28 
29 namespace Imp {
30 #ifdef USE_FACTORY_AUTOREGISTER
31  namespace {
32  bool init = Fopr_Wilson::register_factory();
33  }
34 #endif
35 
36  const std::string Fopr_Wilson::class_name = "Imp::Fopr_Wilson";
37 
38 //====================================================================
39  void Fopr_Wilson::init(const Parameters& params)
40  {
41  ThreadManager::assert_single_thread(class_name);
42 
43  m_vl = CommonParameters::Vlevel();
44 
45  vout.general(m_vl, "%s: construction\n", class_name.c_str());
46  vout.increase_indent();
47 
48  setup();
49 
50  std::string repr;
51  if (!params.fetch_string("gamma_matrix_type", repr)) {
52  m_repr = repr;
53  } else {
54  m_repr = "Dirac"; // default gamma-matrix type
55  vout.general(m_vl, "gamma_matrix_type is not given: defalt = %s\n",
56  m_repr.c_str());
57  }
58  if ((m_repr != "Dirac") && (m_repr != "Chiral")) {
59  vout.crucial("Error at %s: unsupported gamma-matrix type: %s\n",
60  class_name.c_str(), m_repr.c_str());
61  exit(EXIT_FAILURE);
62  }
63 
64  set_parameters(params);
65 
66  vout.decrease_indent();
67  vout.general(m_vl, "%s: construction finished.\n",
68  class_name.c_str());
69  }
70 
71 
72 //====================================================================
73  void Fopr_Wilson::init(const std::string repr)
74  {
75  ThreadManager::assert_single_thread(class_name);
76 
77  m_vl = CommonParameters::Vlevel();
78 
79  vout.general(m_vl, "%s: construction (obsolete)\n",
80  class_name.c_str());
81  vout.increase_indent();
82 
83  setup();
84 
85  m_repr = repr;
86  vout.general(m_vl, "gamma-matrix type is set to %s\n",
87  m_repr.c_str());
88 
89  vout.decrease_indent();
90  vout.general(m_vl, "%s: construction finished.\n",
91  class_name.c_str());
92  }
93 
94 
95 //====================================================================
96  void Fopr_Wilson::setup()
97  {
98  // check_Nc();
99  std::string imple_gauge = imple_Nc();
100  vout.general(m_vl, "Gauge group implementation: %s\n",
101  imple_gauge.c_str());
102 
103  m_Nc = CommonParameters::Nc();
104  m_Nd = CommonParameters::Nd();
105  m_Nvc = 2 * m_Nc;
106  m_Ndf = 2 * m_Nc * m_Nc;
107 
108  m_Nx = CommonParameters::Nx();
109  m_Ny = CommonParameters::Ny();
110  m_Nz = CommonParameters::Nz();
111  m_Nt = CommonParameters::Nt();
112 
113  m_Nvol = CommonParameters::Nvol();
114  m_Ndim = CommonParameters::Ndim();
115  m_boundary.resize(m_Ndim);
116  m_boundary_each_node.resize(m_Ndim);
117 
118  m_U = 0;
119 
120  const int Nvx = m_Nvc * 2 * m_Ny * m_Nz * m_Nt;
121  vcp1_xp = new double[Nvx];
122  vcp2_xp = new double[Nvx];
123  vcp1_xm = new double[Nvx];
124  vcp2_xm = new double[Nvx];
125 
126  const int Nvy = m_Nvc * 2 * m_Nx * m_Nz * m_Nt;
127  vcp1_yp = new double[Nvy];
128  vcp2_yp = new double[Nvy];
129  vcp1_ym = new double[Nvy];
130  vcp2_ym = new double[Nvy];
131 
132  const int Nvz = m_Nvc * 2 * m_Nx * m_Ny * m_Nt;
133  vcp1_zp = new double[Nvz];
134  vcp2_zp = new double[Nvz];
135  vcp1_zm = new double[Nvz];
136  vcp2_zm = new double[Nvz];
137 
138  const int Nvt = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
139  vcp1_tp = new double[Nvt];
140  vcp2_tp = new double[Nvt];
141  vcp1_tm = new double[Nvt];
142  vcp2_tm = new double[Nvt];
143 
144  m_w1.reset(m_Nvc * m_Nd, m_Nvol, 1);
145  m_w2.reset(m_Nvc * m_Nd, m_Nvol, 1);
146  }
147 
148 
149 //====================================================================
150  void Fopr_Wilson::tidyup()
151  {
152  delete[] vcp1_xp;
153  delete[] vcp2_xp;
154  delete[] vcp1_xm;
155  delete[] vcp2_xm;
156 
157  delete[] vcp1_yp;
158  delete[] vcp2_yp;
159  delete[] vcp1_ym;
160  delete[] vcp2_ym;
161 
162  delete[] vcp1_zp;
163  delete[] vcp2_zp;
164  delete[] vcp1_zm;
165  delete[] vcp2_zm;
166 
167  delete[] vcp1_tp;
168  delete[] vcp2_tp;
169  delete[] vcp1_tm;
170  delete[] vcp2_tm;
171  }
172 
173 
174 //====================================================================
175  void Fopr_Wilson::set_parameters(const Parameters& params)
176  {
177 #pragma omp barrier
178  int ith = ThreadManager::get_thread_id();
179  std::string vlevel;
180  if (!params.fetch_string("verbose_level", vlevel)) {
181  if (ith == 0) m_vl = vout.set_verbose_level(vlevel);
182  }
183 #pragma omp barrier
184 
185  //- fetch and check input parameters
186  double kappa;
187  std::vector<int> bc;
188 
189  int err = 0;
190  err += params.fetch_double("hopping_parameter", kappa);
191  err += params.fetch_int_vector("boundary_condition", bc);
192 
193  if (err) {
194  vout.crucial("Error at %s: input parameter not found.\n",
195  class_name.c_str());
196  exit(EXIT_FAILURE);
197  }
198 
199  set_parameters(kappa, bc);
200  }
201 
202 
203 //====================================================================
204  void Fopr_Wilson::set_parameters(const double kappa,
205  const std::vector<int> bc)
206  {
207  assert(bc.size() == m_Ndim);
208 
209 #pragma omp barrier
210 
211  int ith = ThreadManager::get_thread_id();
212  if (ith == 0) {
213  m_kappa = kappa;
214  m_boundary = bc;
215 
216  for (int idir = 0; idir < m_Ndim; ++idir) {
217  m_boundary_each_node[idir] = 1.0;
218  if (Communicator::ipe(idir) == 0) {
219  m_boundary_each_node[idir] = m_boundary[idir];
220  }
221  }
222  }
223 #pragma omp barrier
224 
225  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
226  vout.general(m_vl, " gamma-matrix type = %s\n", m_repr.c_str());
227  vout.general(m_vl, " kappa = %12.8f\n", m_kappa);
228  for (int mu = 0; mu < m_Ndim; ++mu) {
229  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
230  }
231  }
232 
233 
234 //====================================================================
235  void Fopr_Wilson::get_parameters(Parameters& params) const
236  {
237  params.set_double("hopping_parameter", m_kappa);
238  params.set_int_vector("boundary_condition", m_boundary);
239  params.set_string("gamma_matrix_type", m_repr);
240 
241  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
242  }
243 
244 
245 //====================================================================
246  void Fopr_Wilson::set_config(Field *U)
247  {
248 #pragma omp barrier
249  m_U = (Field_G *)U;
250 #pragma omp barrier
251  }
252 
253 
254 //====================================================================
255  void Fopr_Wilson::set_mode(const std::string mode)
256  {
257 #pragma omp barrier
258  int ith = ThreadManager::get_thread_id();
259  if (ith == 0) m_mode = mode;
260 #pragma omp barrier
261  }
262 
263 
264 //====================================================================
265  void Fopr_Wilson::mult(Field& v, const Field& w)
266  {
267  if (m_mode == "D") {
268  D(v, w);
269  } else if (m_mode == "Ddag") {
270  Ddag(v, w);
271  } else if (m_mode == "DdagD") {
272  DdagD(v, w);
273  } else if (m_mode == "DDdag") {
274  DDdag(v, w);
275  } else if (m_mode == "H") {
276  H(v, w);
277  } else {
278  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
279  class_name.c_str(), m_mode.c_str());
280  exit(EXIT_FAILURE);
281  }
282  }
283 
284 
285 //====================================================================
286  void Fopr_Wilson::mult_dag(Field& v, const Field& w)
287  {
288  if (m_mode == "D") {
289  Ddag(v, w);
290  } else if (m_mode == "Ddag") {
291  D(v, w);
292  } else if (m_mode == "DdagD") {
293  DdagD(v, w);
294  } else if (m_mode == "DDdag") {
295  DDdag(v, w);
296  } else if (m_mode == "H") {
297  H(v, w);
298  } else {
299  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
300  class_name.c_str(), m_mode.c_str());
301  exit(EXIT_FAILURE);
302  }
303  }
304 
305 
306 //====================================================================
307  void Fopr_Wilson::mult(Field& v, const Field& w,
308  const std::string mode)
309  {
310  if (mode == "D") {
311  D(v, w);
312  } else if (mode == "Ddag") {
313  Ddag(v, w);
314  } else if (mode == "DdagD") {
315  DdagD(v, w);
316  } else if (mode == "DDdag") {
317  DDdag(v, w);
318  } else if (mode == "H") {
319  H(v, w);
320  } else {
321  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
322  class_name.c_str(), mode.c_str());
323  exit(EXIT_FAILURE);
324  }
325  }
326 
327 
328 //====================================================================
329  void Fopr_Wilson::mult_dag(Field& v, const Field& w,
330  const std::string mode)
331  {
332  if (mode == "D") {
333  Ddag(v, w);
334  } else if (mode == "Ddag") {
335  D(v, w);
336  } else if (mode == "DdagD") {
337  DdagD(v, w);
338  } else if (mode == "DDdag") {
339  DDdag(v, w);
340  } else if (mode == "H") {
341  H(v, w);
342  } else {
343  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
344  class_name.c_str(), mode.c_str());
345  exit(EXIT_FAILURE);
346  }
347  }
348 
349 
350 //====================================================================
351  void Fopr_Wilson::mult_up(const int mu, Field& v, const Field& w)
352  {
353  if (mu == 0) {
354  mult_xp(v, w);
355  } else if (mu == 1) {
356  mult_yp(v, w);
357  } else if (mu == 2) {
358  mult_zp(v, w);
359  } else if (mu == 3) {
360  if (m_repr == "Dirac") {
361  mult_tp_dirac(v, w);
362  } else {
363  mult_tp_chiral(v, w);
364  }
365  } else {
366  vout.crucial(m_vl, "Error at %s::mult_up: illegal mu=%d\n",
367  class_name.c_str(), mu);
368  exit(EXIT_FAILURE);
369  }
370  }
371 
372 
373 //====================================================================
374  void Fopr_Wilson::mult_dn(const int mu, Field& v, const Field& w)
375  {
376  if (mu == 0) {
377  mult_xm(v, w);
378  } else if (mu == 1) {
379  mult_ym(v, w);
380  } else if (mu == 2) {
381  mult_zm(v, w);
382  } else if (mu == 3) {
383  if (m_repr == "Dirac") {
384  mult_tm_dirac(v, w);
385  } else {
386  mult_tm_chiral(v, w);
387  }
388  } else {
389  vout.crucial(m_vl, "Error at %s::mult_dn: illegal mu=%d\n",
390  class_name.c_str(), mu);
391  exit(EXIT_FAILURE);
392  }
393  }
394 
395 
396 //====================================================================
397  void Fopr_Wilson::mult_gm5(Field& v, const Field& w)
398  {
399  if (m_repr == "Dirac") {
400  mult_gm5_dirac(v, w);
401  } else if (m_repr == "Chiral") {
402  mult_gm5_chiral(v, w);
403  }
404  }
405 
406 
407 //====================================================================
408  void Fopr_Wilson::D(Field& v, const Field& w)
409  {
410  if (m_repr == "Dirac") {
411  D_ex_dirac(v, 0, w, 0);
412  //D_ex_dirac_alt(v, 0, w, 0);
413  } else if (m_repr == "Chiral") {
414  D_ex_chiral(v, 0, w, 0);
415  //D_ex_chiral_alt(v, 0, w, 0);
416  }
417  }
418 
419 
420 //====================================================================
421  void Fopr_Wilson::D_ex(Field& v, const int ex1,
422  const Field& w, const int ex2)
423  {
424  if (m_repr == "Dirac") {
425  D_ex_dirac(v, ex1, w, ex2);
426  } else if (m_repr == "Chiral") {
427  D_ex_chiral(v, ex1, w, ex2);
428  }
429  }
430 
431 
432 //====================================================================
433  void Fopr_Wilson::Ddag(Field& v, const Field& w)
434  {
435  mult_gm5(v, w);
436  D(m_w1, v);
437  mult_gm5(v, m_w1);
438  }
439 
440 
441 //====================================================================
442  void Fopr_Wilson::DdagD(Field& v, const Field& w)
443  {
444  D(m_w1, w);
445  mult_gm5(v, m_w1);
446  D(m_w1, v);
447  mult_gm5(v, m_w1);
448  }
449 
450 
451 //====================================================================
452  void Fopr_Wilson::DDdag(Field& v, const Field& w)
453  {
454  mult_gm5(m_w1, w);
455  D(v, m_w1);
456  mult_gm5(m_w1, v);
457  D(v, m_w1);
458  }
459 
460 
461 //====================================================================
462  void Fopr_Wilson::H(Field& v, const Field& w)
463  {
464  D(m_w1, w);
465  mult_gm5(v, m_w1);
466  }
467 
468 
469 //====================================================================
470  void Fopr_Wilson::D_ex_dirac(Field& v, const int ex1,
471  const Field& w, const int ex2)
472  {
473  const int Ninvol = m_Nvc * m_Nd * m_Nvol;
474  double *RESTRICT vp = v.ptr(Ninvol * ex1);
475  const double *RESTRICT wp = w.ptr(Ninvol * ex2);
476  const double *RESTRICT up = m_U->ptr(0);
477 
478  int ith, nth, is, ns;
479  set_threadtask(ith, nth, is, ns, m_Nvol);
480 
481  int Nvc2 = m_Nvc * 2;
482  int Nvcd = m_Nvc * m_Nd;
483 
484  int Nxy = m_Nx * m_Ny;
485  int Nxyz = m_Nx * m_Ny * m_Nz;
486 
487 #pragma omp barrier
488 
489  for (int site = is; site < ns; ++site) {
490  int ix = site % m_Nx;
491  int iyzt = site / m_Nx;
492  int iy = iyzt % m_Ny;
493  int izt = iyzt / m_Ny;
494  int iz = izt % m_Nz;
495  int it = izt / m_Nz;
496 
497  int ixy = ix + m_Nx * iy;
498  int ixyz = ixy + Nxy * iz;
499  int ixyt = ixy + Nxy * it;
500  int ixzt = ix + m_Nx * izt;
501 
502  int in = Nvcd * site;
503 
504  if (ix == 0) {
505  int ix1 = Nvc2 * iyzt;
506  int ix2 = ix1 + m_Nvc;
507  double bc2 = m_boundary_each_node[0];
508  double vt1[NVC], vt2[NVC];
509  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
510  for (int ivc = 0; ivc < NVC; ++ivc) {
511  vcp1_xp[ivc + ix1] = bc2 * vt1[ivc];
512  vcp1_xp[ivc + ix2] = bc2 * vt2[ivc];
513  }
514  }
515 
516  if (ix == m_Nx - 1) {
517  int ig = m_Ndf * (site + 0 * m_Nvol);
518  int ix1 = Nvc2 * iyzt;
519  int ix2 = ix1 + NVC;
520  double vt1[NVC], vt2[NVC];
521  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
522  for (int ic = 0; ic < m_Nc; ++ic) {
523  int ic2 = 2 * ic;
524  int icr = 2 * ic;
525  int ici = 2 * ic + 1;
526  vcp1_xm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
527  vcp1_xm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
528  vcp1_xm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
529  vcp1_xm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
530  }
531  }
532 
533  if (iy == 0) {
534  int ix1 = Nvc2 * ixzt;
535  int ix2 = ix1 + NVC;
536  double bc2 = m_boundary_each_node[1];
537  double vt1[NVC], vt2[NVC];
538  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
539  for (int ivc = 0; ivc < NVC; ++ivc) {
540  vcp1_yp[ivc + ix1] = bc2 * vt1[ivc];
541  vcp1_yp[ivc + ix2] = bc2 * vt2[ivc];
542  }
543  }
544 
545  if (iy == m_Ny - 1) {
546  int ig = m_Ndf * (site + 1 * m_Nvol);
547  int ix1 = Nvc2 * ixzt;
548  int ix2 = ix1 + NVC;
549  double vt1[NVC], vt2[NVC];
550  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
551  for (int ic = 0; ic < m_Nc; ++ic) {
552  int ic2 = 2 * ic;
553  int icr = 2 * ic;
554  int ici = 2 * ic + 1;
555  vcp1_ym[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
556  vcp1_ym[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
557  vcp1_ym[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
558  vcp1_ym[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
559  }
560  }
561 
562  if (iz == 0) {
563  int ix1 = Nvc2 * ixyt;
564  int ix2 = ix1 + NVC;
565  double bc2 = m_boundary_each_node[2];
566  double vt1[NVC], vt2[NVC];
567  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
568  for (int ivc = 0; ivc < NVC; ++ivc) {
569  vcp1_zp[ivc + ix1] = bc2 * vt1[ivc];
570  vcp1_zp[ivc + ix2] = bc2 * vt2[ivc];
571  }
572  }
573 
574  if (iz == m_Nz - 1) {
575  int ig = m_Ndf * (site + 2 * m_Nvol);
576  int ix1 = Nvc2 * ixyt;
577  int ix2 = ix1 + NVC;
578  double vt1[NVC], vt2[NVC];
579  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
580  for (int ic = 0; ic < m_Nc; ++ic) {
581  int ic2 = 2 * ic;
582  int icr = 2 * ic;
583  int ici = 2 * ic + 1;
584  vcp1_zm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
585  vcp1_zm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
586  vcp1_zm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
587  vcp1_zm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
588  }
589  }
590 
591  if (it == 0) {
592  int ix1 = Nvc2 * ixyz;
593  int ix2 = ix1 + NVC;
594  double bc2 = m_boundary_each_node[3];
595  double vt1[NVC], vt2[NVC];
596  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
597  for (int ivc = 0; ivc < NVC; ++ivc) {
598  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
599  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
600  }
601  }
602 
603  if (it == m_Nt - 1) {
604  int ig = m_Ndf * (site + 3 * m_Nvol);
605  int ix1 = Nvc2 * ixyz;
606  int ix2 = ix1 + NVC;
607  double vt1[NVC], vt2[NVC];
608  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
609  for (int ic = 0; ic < m_Nc; ++ic) {
610  int ic2 = 2 * ic;
611  int icr = 2 * ic;
612  int ici = 2 * ic + 1;
613  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
614  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
615  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
616  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
617  }
618  }
619  }
620 
621 #pragma omp barrier
622 
623 #pragma omp master
624  {
625  int Nvx = m_Nvc * 2 * m_Ny * m_Nz * m_Nt;
626  Communicator::exchange(Nvx, vcp2_xp, vcp1_xp, 0, 1, 1);
627  Communicator::exchange(Nvx, vcp2_xm, vcp1_xm, 0, -1, 2);
628 
629  int Nvy = m_Nvc * 2 * m_Nx * m_Nz * m_Nt;
630  Communicator::exchange(Nvy, vcp2_yp, vcp1_yp, 1, 1, 3);
631  Communicator::exchange(Nvy, vcp2_ym, vcp1_ym, 1, -1, 4);
632 
633  int Nvz = m_Nvc * 2 * m_Nx * m_Ny * m_Nt;
634  Communicator::exchange(Nvz, vcp2_zp, vcp1_zp, 2, 1, 5);
635  Communicator::exchange(Nvz, vcp2_zm, vcp1_zm, 2, -1, 6);
636 
637  int Nvt = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
638  Communicator::exchange(Nvt, vcp2_tp, vcp1_tp, 3, 1, 7);
639  Communicator::exchange(Nvt, vcp2_tm, vcp1_tm, 3, -1, 8);
640  }
641 #pragma omp barrier
642 
643  for (int site = is; site < ns; ++site) {
644  int ix = site % m_Nx;
645  int iyzt = site / m_Nx;
646  int iy = iyzt % m_Ny;
647  int izt = iyzt / m_Ny;
648  int iz = izt % m_Nz;
649  int it = izt / m_Nz;
650 
651  int ixy = site % Nxy;
652  int ixyz = ixy + Nxy * iz;
653  int ixyt = ixy + Nxy * it;
654  int ixzt = ix + m_Nx * izt;
655 
656  int iv = Nvcd * site;
657 
658  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
659  vp[ivcd + iv] = 0.0;
660  }
661 
662  if (ix < m_Nx - 1) {
663  int nei = ix + 1 + m_Nx * iyzt;
664  int in = Nvcd * nei;
665  int ig = m_Ndf * (site + 0 * m_Nvol);
666  double vt1[NVC], vt2[NVC];
667  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
668  for (int ic = 0; ic < m_Nc; ++ic) {
669  int ic2 = ic * NVC;
670  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
671  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
672  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
673  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
674  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
675  }
676  } else {
677  int ix1 = Nvc2 * iyzt;
678  int ix2 = ix1 + NVC;
679  int ig = m_Ndf * (site + 0 * m_Nvol);
680  for (int ic = 0; ic < m_Nc; ++ic) {
681  int ic2 = ic * NVC;
682  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
683  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
684  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
685  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
686  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
687  }
688  }
689 
690  if (ix > 0) {
691  int nei = ix - 1 + m_Nx * iyzt;
692  int ig = m_Ndf * (nei + 0 * m_Nvol);
693  int in = Nvcd * nei;
694  double vt1[NVC], vt2[NVC];
695  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
696  for (int ic = 0; ic < m_Nc; ++ic) {
697  int ic2 = 2 * ic;
698  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
699  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
700  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
701  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
702  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
703  }
704  } else {
705  int ix1 = Nvc2 * iyzt;
706  int ix2 = ix1 + NVC;
707  double bc2 = m_boundary_each_node[0];
708  for (int ic = 0; ic < m_Nc; ++ic) {
709  double wt1r = bc2 * vcp2_xm[2 * ic + ix1];
710  double wt1i = bc2 * vcp2_xm[2 * ic + 1 + ix1];
711  double wt2r = bc2 * vcp2_xm[2 * ic + ix2];
712  double wt2i = bc2 * vcp2_xm[2 * ic + 1 + ix2];
713  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
714  }
715  }
716 
717  if (iy < m_Ny - 1) {
718  int nei = ix + m_Nx * (iy + 1 + m_Ny * izt);
719  int ig = m_Ndf * (site + 1 * m_Nvol);
720  int in = Nvcd * nei;
721  double vt1[NVC], vt2[NVC];
722  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
723  for (int ic = 0; ic < m_Nc; ++ic) {
724  int ic2 = ic * NVC;
725  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
726  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
727  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
728  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
729  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
730  }
731  } else {
732  int ig = m_Ndf * (site + 1 * m_Nvol);
733  int ix1 = Nvc2 * ixzt;
734  int ix2 = ix1 + NVC;
735  for (int ic = 0; ic < m_Nc; ++ic) {
736  int ic2 = ic * NVC;
737  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
738  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
739  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
740  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
741  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
742  }
743  }
744 
745  if (iy > 0) {
746  int nei = ix + m_Nx * (iy - 1 + m_Ny * izt);
747  int ig = m_Ndf * (nei + 1 * m_Nvol);
748  int in = Nvcd * nei;
749  double vt1[NVC], vt2[NVC];
750  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
751  for (int ic = 0; ic < m_Nc; ++ic) {
752  int ic2 = 2 * ic;
753  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
754  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
755  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
756  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
757  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
758  }
759  } else {
760  int ix1 = Nvc2 * ixzt;
761  int ix2 = ix1 + NVC;
762  double bc2 = m_boundary_each_node[1];
763  for (int ic = 0; ic < m_Nc; ++ic) {
764  double wt1r = bc2 * vcp2_ym[2 * ic + ix1];
765  double wt1i = bc2 * vcp2_ym[2 * ic + 1 + ix1];
766  double wt2r = bc2 * vcp2_ym[2 * ic + ix2];
767  double wt2i = bc2 * vcp2_ym[2 * ic + 1 + ix2];
768  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
769  }
770  }
771 
772  if (iz < m_Nz - 1) {
773  int nei = ixy + Nxy * (iz + 1 + m_Nz * it);
774  int ig = m_Ndf * (site + 2 * m_Nvol);
775  int in = Nvcd * nei;
776  double vt1[NVC], vt2[NVC];
777  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
778  for (int ic = 0; ic < m_Nc; ++ic) {
779  int ic2 = ic * NVC;
780  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
781  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
782  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
783  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
784  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
785  }
786  } else {
787  int ig = m_Ndf * (site + 2 * m_Nvol);
788  int ix1 = Nvc2 * ixyt;
789  int ix2 = ix1 + NVC;
790  for (int ic = 0; ic < m_Nc; ++ic) {
791  int ic2 = ic * NVC;
792  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
793  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
794  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
795  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
796  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
797  }
798  }
799 
800  if (iz > 0) {
801  int nei = ixy + Nxy * (iz - 1 + m_Nz * it);
802  int ig = m_Ndf * (nei + 2 * m_Nvol);
803  int in = Nvcd * nei;
804  double vt1[NVC], vt2[NVC];
805  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
806  for (int ic = 0; ic < m_Nc; ++ic) {
807  int ic2 = 2 * ic;
808  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
809  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
810  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
811  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
812  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
813  }
814  } else {
815  int ix1 = Nvc2 * ixyt;
816  int ix2 = ix1 + NVC;
817  double bc2 = m_boundary_each_node[2];
818  for (int ic = 0; ic < m_Nc; ++ic) {
819  double wt1r = bc2 * vcp2_zm[2 * ic + ix1];
820  double wt1i = bc2 * vcp2_zm[2 * ic + 1 + ix1];
821  double wt2r = bc2 * vcp2_zm[2 * ic + ix2];
822  double wt2i = bc2 * vcp2_zm[2 * ic + 1 + ix2];
823  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
824  }
825  }
826 
827  if (it < m_Nt - 1) {
828  int nei = ixyz + Nxyz * (it + 1);
829  int ig = m_Ndf * (site + 3 * m_Nvol);
830  int in = Nvcd * nei;
831  double vt1[NVC], vt2[NVC];
832  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
833  for (int ic = 0; ic < m_Nc; ++ic) {
834  int ic2 = ic * NVC;
835  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
836  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
837  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
838  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
839  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
840  }
841  } else {
842  int ig = m_Ndf * (site + 3 * m_Nvol);
843  int ix1 = Nvc2 * ixyz;
844  int ix2 = ix1 + NVC;
845  for (int ic = 0; ic < m_Nc; ++ic) {
846  int ic2 = ic * NVC;
847  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
848  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
849  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
850  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
851  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
852  }
853  }
854 
855  if (it > 0) {
856  int nei = ixyz + Nxyz * (it - 1);
857  int ig = m_Ndf * (nei + 3 * m_Nvol);
858  int in = Nvcd * nei;
859  double vt1[NVC], vt2[NVC];
860  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
861  for (int ic = 0; ic < m_Nc; ++ic) {
862  int ic2 = 2 * ic;
863  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
864  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
865  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
866  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
867  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
868  }
869  } else {
870  int ix1 = Nvc2 * ixyz;
871  int ix2 = ix1 + NVC;
872  double bc2 = m_boundary_each_node[3];
873  for (int ic = 0; ic < m_Nc; ++ic) {
874  int icr = 2 * ic;
875  int ici = 2 * ic + 1;
876  double wt1r = bc2 * vcp2_tm[icr + ix1];
877  double wt1i = bc2 * vcp2_tm[ici + ix1];
878  double wt2r = bc2 * vcp2_tm[icr + ix2];
879  double wt2i = bc2 * vcp2_tm[ici + ix2];
880  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
881  }
882  }
883 
884  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
885  vp[ivcd + iv] = -m_kappa * vp[ivcd + iv] + wp[ivcd + iv];
886  }
887  }
888 
889 #pragma omp barrier
890  }
891 
892 
893 //====================================================================
894  void Fopr_Wilson::D_ex_chiral(Field& v, const int ex1,
895  const Field& w, const int ex2)
896  {
897  const int Ninvol = m_Nvc * m_Nd * m_Nvol;
898  double *RESTRICT vp = v.ptr(Ninvol * ex1);
899  const double *RESTRICT wp = w.ptr(Ninvol * ex2);
900  const double *RESTRICT up = m_U->ptr(0);
901 
902  int ith, nth, is, ns;
903  set_threadtask(ith, nth, is, ns, m_Nvol);
904 
905  int Nvc2 = m_Nvc * 2;
906  int Nvcd = m_Nvc * m_Nd;
907 
908  int Nxy = m_Nx * m_Ny;
909  int Nxyz = m_Nx * m_Ny * m_Nz;
910 
911 #pragma omp barrier
912 
913  for (int site = is; site < ns; ++site) {
914  int ix = site % m_Nx;
915  int iyzt = site / m_Nx;
916  int iy = iyzt % m_Ny;
917  int izt = iyzt / m_Ny;
918  int iz = izt % m_Nz;
919  int it = izt / m_Nz;
920 
921  int ixy = ix + m_Nx * iy;
922  int ixyz = ixy + Nxy * iz;
923  int ixyt = ixy + Nxy * it;
924  int ixzt = ix + m_Nx * izt;
925 
926  int in = Nvcd * site;
927 
928  if (ix == 0) {
929  int ix1 = Nvc2 * iyzt;
930  int ix2 = ix1 + NVC;
931  double bc2 = m_boundary_each_node[0];
932  double vt1[NVC], vt2[NVC];
933  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
934  for (int ivc = 0; ivc < NVC; ++ivc) {
935  vcp1_xp[ivc + ix1] = bc2 * vt1[ivc];
936  vcp1_xp[ivc + ix2] = bc2 * vt2[ivc];
937  }
938  }
939 
940  if (ix == m_Nx - 1) {
941  int ig = m_Ndf * (site + 0 * m_Nvol);
942  int ix1 = Nvc2 * iyzt;
943  int ix2 = ix1 + NVC;
944  double vt1[NVC], vt2[NVC];
945  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
946  for (int ic = 0; ic < m_Nc; ++ic) {
947  int ic2 = 2 * ic;
948  int icr = 2 * ic;
949  int ici = 2 * ic + 1;
950  vcp1_xm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
951  vcp1_xm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
952  vcp1_xm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
953  vcp1_xm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
954  }
955  }
956 
957  if (iy == 0) {
958  int ix1 = Nvc2 * ixzt;
959  int ix2 = ix1 + NVC;
960  double bc2 = m_boundary_each_node[1];
961  double vt1[NVC], vt2[NVC];
962  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
963  for (int ivc = 0; ivc < NVC; ++ivc) {
964  vcp1_yp[ivc + ix1] = bc2 * vt1[ivc];
965  vcp1_yp[ivc + ix2] = bc2 * vt2[ivc];
966  }
967  }
968 
969  if (iy == m_Ny - 1) {
970  int ig = m_Ndf * (site + 1 * m_Nvol);
971  int ix1 = Nvc2 * ixzt;
972  int ix2 = ix1 + NVC;
973  double vt1[NVC], vt2[NVC];
974  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
975  for (int ic = 0; ic < m_Nc; ++ic) {
976  int ic2 = 2 * ic;
977  int icr = 2 * ic;
978  int ici = 2 * ic + 1;
979  vcp1_ym[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
980  vcp1_ym[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
981  vcp1_ym[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
982  vcp1_ym[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
983  }
984  }
985 
986  if (iz == 0) {
987  int ix1 = Nvc2 * ixyt;
988  int ix2 = ix1 + NVC;
989  double bc2 = m_boundary_each_node[2];
990  double vt1[NVC], vt2[NVC];
991  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
992  for (int ivc = 0; ivc < NVC; ++ivc) {
993  vcp1_zp[ivc + ix1] = bc2 * vt1[ivc];
994  vcp1_zp[ivc + ix2] = bc2 * vt2[ivc];
995  }
996  }
997 
998  if (iz == m_Nz - 1) {
999  int ig = m_Ndf * (site + 2 * m_Nvol);
1000  int ix1 = Nvc2 * ixyt;
1001  int ix2 = ix1 + NVC;
1002  double vt1[NVC], vt2[NVC];
1003  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
1004  for (int ic = 0; ic < m_Nc; ++ic) {
1005  int ic2 = 2 * ic;
1006  int icr = 2 * ic;
1007  int ici = 2 * ic + 1;
1008  vcp1_zm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1009  vcp1_zm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1010  vcp1_zm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1011  vcp1_zm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1012  }
1013  }
1014 
1015  if (it == 0) {
1016  int ix1 = Nvc2 * ixyz;
1017  int ix2 = ix1 + NVC;
1018  double bc2 = m_boundary_each_node[3];
1019  double vt1[NVC], vt2[NVC];
1020  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
1021  for (int ivc = 0; ivc < NVC; ++ivc) {
1022  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
1023  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
1024  }
1025  }
1026 
1027  if (it == m_Nt - 1) {
1028  int ig = m_Ndf * (site + 3 * m_Nvol);
1029  int ix1 = Nvc2 * ixyz;
1030  int ix2 = ix1 + NVC;
1031  double vt1[NVC], vt2[NVC];
1032  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
1033  for (int ic = 0; ic < m_Nc; ++ic) {
1034  int ic2 = 2 * ic;
1035  int icr = 2 * ic;
1036  int ici = 2 * ic + 1;
1037  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1038  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1039  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1040  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1041  }
1042  }
1043  }
1044 
1045 #pragma omp barrier
1046 
1047 #pragma omp master
1048  {
1049  int Nvx = m_Nvc * 2 * m_Ny * m_Nz * m_Nt;
1050  Communicator::exchange(Nvx, vcp2_xp, vcp1_xp, 0, 1, 1);
1051  Communicator::exchange(Nvx, vcp2_xm, vcp1_xm, 0, -1, 2);
1052 
1053  int Nvy = m_Nvc * 2 * m_Nx * m_Nz * m_Nt;
1054  Communicator::exchange(Nvy, vcp2_yp, vcp1_yp, 1, 1, 3);
1055  Communicator::exchange(Nvy, vcp2_ym, vcp1_ym, 1, -1, 4);
1056 
1057  int Nvz = m_Nvc * 2 * m_Nx * m_Ny * m_Nt;
1058  Communicator::exchange(Nvz, vcp2_zp, vcp1_zp, 2, 1, 5);
1059  Communicator::exchange(Nvz, vcp2_zm, vcp1_zm, 2, -1, 6);
1060 
1061  int Nvt = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
1062  Communicator::exchange(Nvt, vcp2_tp, vcp1_tp, 3, 1, 7);
1063  Communicator::exchange(Nvt, vcp2_tm, vcp1_tm, 3, -1, 8);
1064  }
1065 #pragma omp barrier
1066 
1067  for (int site = is; site < ns; ++site) {
1068  int ix = site % m_Nx;
1069  int iyzt = site / m_Nx;
1070  int iy = iyzt % m_Ny;
1071  int izt = iyzt / m_Ny;
1072  int iz = izt % m_Nz;
1073  int it = izt / m_Nz;
1074 
1075  int ixy = site % Nxy;
1076  int ixyz = ixy + Nxy * iz;
1077  int ixyt = ixy + Nxy * it;
1078  int ixzt = ix + m_Nx * izt;
1079 
1080  int iv = Nvcd * site;
1081 
1082  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
1083  vp[ivcd + iv] = 0.0;
1084  }
1085 
1086  if (ix < m_Nx - 1) {
1087  int nei = ix + 1 + m_Nx * iyzt;
1088  int in = Nvcd * nei;
1089  int ig = m_Ndf * (site + 0 * m_Nvol);
1090  double vt1[NVC], vt2[NVC];
1091  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
1092  for (int ic = 0; ic < m_Nc; ++ic) {
1093  int ic2 = ic * NVC;
1094  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1095  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1096  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1097  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1098  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1099  }
1100  } else {
1101  int ix1 = Nvc2 * iyzt;
1102  int ix2 = ix1 + NVC;
1103  int ig = m_Ndf * (site + 0 * m_Nvol);
1104  for (int ic = 0; ic < m_Nc; ++ic) {
1105  int ic2 = ic * NVC;
1106  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
1107  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
1108  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
1109  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
1110  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1111  }
1112  }
1113 
1114  if (ix > 0) {
1115  int nei = ix - 1 + m_Nx * iyzt;
1116  int ig = m_Ndf * (nei + 0 * m_Nvol);
1117  int in = Nvcd * nei;
1118  double vt1[NVC], vt2[NVC];
1119  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
1120  for (int ic = 0; ic < m_Nc; ++ic) {
1121  int ic2 = 2 * ic;
1122  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1123  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1124  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1125  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1126  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1127  }
1128  } else {
1129  int ix1 = Nvc2 * iyzt;
1130  int ix2 = ix1 + NVC;
1131  double bc2 = m_boundary_each_node[0];
1132  for (int ic = 0; ic < m_Nc; ++ic) {
1133  double wt1r = bc2 * vcp2_xm[2 * ic + ix1];
1134  double wt1i = bc2 * vcp2_xm[2 * ic + 1 + ix1];
1135  double wt2r = bc2 * vcp2_xm[2 * ic + ix2];
1136  double wt2i = bc2 * vcp2_xm[2 * ic + 1 + ix2];
1137  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1138  }
1139  }
1140 
1141  if (iy < m_Ny - 1) {
1142  int nei = ix + m_Nx * (iy + 1 + m_Ny * izt);
1143  int ig = m_Ndf * (site + 1 * m_Nvol);
1144  int in = Nvcd * nei;
1145  double vt1[NVC], vt2[NVC];
1146  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
1147  for (int ic = 0; ic < m_Nc; ++ic) {
1148  int ic2 = ic * NVC;
1149  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1150  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1151  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1152  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1153  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1154  }
1155  } else {
1156  int ig = m_Ndf * (site + 1 * m_Nvol);
1157  int ix1 = Nvc2 * ixzt;
1158  int ix2 = ix1 + NVC;
1159  for (int ic = 0; ic < m_Nc; ++ic) {
1160  int ic2 = ic * NVC;
1161  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1162  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1163  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1164  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1165  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1166  }
1167  }
1168 
1169  if (iy > 0) {
1170  int nei = ix + m_Nx * (iy - 1 + m_Ny * izt);
1171  int ig = m_Ndf * (nei + 1 * m_Nvol);
1172  int in = Nvcd * nei;
1173  double vt1[NVC], vt2[NVC];
1174  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
1175  for (int ic = 0; ic < m_Nc; ++ic) {
1176  int ic2 = 2 * ic;
1177  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1178  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1179  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1180  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1181  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1182  }
1183  } else {
1184  int ix1 = Nvc2 * ixzt;
1185  int ix2 = ix1 + NVC;
1186  double bc2 = m_boundary_each_node[1];
1187  for (int ic = 0; ic < m_Nc; ++ic) {
1188  double wt1r = bc2 * vcp2_ym[2 * ic + ix1];
1189  double wt1i = bc2 * vcp2_ym[2 * ic + 1 + ix1];
1190  double wt2r = bc2 * vcp2_ym[2 * ic + ix2];
1191  double wt2i = bc2 * vcp2_ym[2 * ic + 1 + ix2];
1192  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1193  }
1194  }
1195 
1196  if (iz < m_Nz - 1) {
1197  int nei = ixy + Nxy * (iz + 1 + m_Nz * it);
1198  int ig = m_Ndf * (site + 2 * m_Nvol);
1199  int in = Nvcd * nei;
1200  double vt1[NVC], vt2[NVC];
1201  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
1202  for (int ic = 0; ic < m_Nc; ++ic) {
1203  int ic2 = ic * NVC;
1204  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1205  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1206  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1207  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1208  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1209  }
1210  } else {
1211  int ig = m_Ndf * (site + 2 * m_Nvol);
1212  int ix1 = Nvc2 * ixyt;
1213  int ix2 = ix1 + NVC;
1214  for (int ic = 0; ic < m_Nc; ++ic) {
1215  int ic2 = ic * NVC;
1216  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1217  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1218  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1219  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1220  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1221  }
1222  }
1223 
1224  if (iz > 0) {
1225  int nei = ixy + Nxy * (iz - 1 + m_Nz * it);
1226  int ig = m_Ndf * (nei + 2 * m_Nvol);
1227  int in = Nvcd * nei;
1228  double vt1[NVC], vt2[NVC];
1229  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
1230  for (int ic = 0; ic < m_Nc; ++ic) {
1231  int ic2 = 2 * ic;
1232  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1233  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1234  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1235  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1236  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1237  }
1238  } else {
1239  int ix1 = Nvc2 * ixyt;
1240  int ix2 = ix1 + NVC;
1241  double bc2 = m_boundary_each_node[2];
1242  for (int ic = 0; ic < m_Nc; ++ic) {
1243  double wt1r = bc2 * vcp2_zm[2 * ic + ix1];
1244  double wt1i = bc2 * vcp2_zm[2 * ic + 1 + ix1];
1245  double wt2r = bc2 * vcp2_zm[2 * ic + ix2];
1246  double wt2i = bc2 * vcp2_zm[2 * ic + 1 + ix2];
1247  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1248  }
1249  }
1250 
1251  if (it < m_Nt - 1) {
1252  int nei = ixyz + Nxyz * (it + 1);
1253  int ig = m_Ndf * (site + 3 * m_Nvol);
1254  int in = Nvcd * nei;
1255  double vt1[NVC], vt2[NVC];
1256  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
1257  for (int ic = 0; ic < m_Nc; ++ic) {
1258  int ic2 = ic * NVC;
1259  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1260  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1261  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1262  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1263  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1264  }
1265  } else {
1266  int ig = m_Ndf * (site + 3 * m_Nvol);
1267  int ix1 = Nvc2 * ixyz;
1268  int ix2 = ix1 + NVC;
1269  for (int ic = 0; ic < m_Nc; ++ic) {
1270  int ic2 = ic * NVC;
1271  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1272  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1273  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1274  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1275  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1276  }
1277  }
1278 
1279  if (it > 0) {
1280  int nei = ixyz + Nxyz * (it - 1);
1281  int ig = m_Ndf * (nei + 3 * m_Nvol);
1282  int in = Nvcd * nei;
1283  double vt1[NVC], vt2[NVC];
1284  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
1285  for (int ic = 0; ic < m_Nc; ++ic) {
1286  int ic2 = 2 * ic;
1287  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1288  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1289  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1290  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1291  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1292  }
1293  } else {
1294  int ix1 = Nvc2 * ixyz;
1295  int ix2 = ix1 + NVC;
1296  double bc2 = m_boundary_each_node[3];
1297  for (int ic = 0; ic < m_Nc; ++ic) {
1298  int icr = 2 * ic;
1299  int ici = 2 * ic + 1;
1300  double wt1r = bc2 * vcp2_tm[icr + ix1];
1301  double wt1i = bc2 * vcp2_tm[ici + ix1];
1302  double wt2r = bc2 * vcp2_tm[icr + ix2];
1303  double wt2i = bc2 * vcp2_tm[ici + ix2];
1304  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1305  }
1306  }
1307 
1308  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
1309  vp[ivcd + iv] = -m_kappa * vp[ivcd + iv] + wp[ivcd + iv];
1310  }
1311  }
1312 
1313 #pragma omp barrier
1314  }
1315 
1316 
1317 //====================================================================
1318  void Fopr_Wilson::D_ex_dirac_alt(Field& w, const int ex1,
1319  const Field& f, const int ex2)
1320  {
1321  clear(w);
1322  mult_xp(w, f);
1323  mult_xm(w, f);
1324  mult_yp(w, f);
1325  mult_ym(w, f);
1326  mult_zp(w, f);
1327  mult_zm(w, f);
1328  mult_tp_dirac(w, f);
1329  mult_tm_dirac(w, f);
1330  daypx(w, -m_kappa, f); // w = -m_kappa * w + f.
1331  }
1332 
1333 
1334 //====================================================================
1335  void Fopr_Wilson::D_ex_chiral_alt(Field& w, const int ex1,
1336  const Field& f, const int ex2)
1337  {
1338  clear(w);
1339  mult_xp(w, f);
1340  mult_xm(w, f);
1341  mult_yp(w, f);
1342  mult_ym(w, f);
1343  mult_zp(w, f);
1344  mult_zm(w, f);
1345  mult_tp_chiral(w, f);
1346  mult_tm_chiral(w, f);
1347  daypx(w, -m_kappa, f); // w = -m_kappa * w + f.
1348  }
1349 
1350 
1351 //====================================================================
1352  void Fopr_Wilson::mult_gm5p(const int mu, Field& v, const Field& w)
1353  {
1354  clear(m_w2);
1355  mult_up(mu, m_w2, w);
1356  mult_gm5(v, m_w2);
1357  }
1358 
1359 
1360 //====================================================================
1361  void Fopr_Wilson::proj_chiral(Field& w, const int ex1,
1362  const Field& v, const int ex2, const int ipm)
1363  {
1364  double fpm = 0.0;
1365 
1366  if (ipm == 1) {
1367  fpm = 1.0;
1368  } else if (ipm == -1) {
1369  fpm = -1.0;
1370  } else {
1371  vout.crucial(m_vl, "Error at %s: illegal chirality = %d\n", class_name.c_str(), ipm);
1372  exit(EXIT_FAILURE);
1373  }
1374 
1375  copy(m_w1, 0, v, ex2);
1376  mult_gm5(m_w2, m_w1);
1377  axpy(m_w1, 0, fpm, m_w2, 0);
1378  copy(w, ex1, m_w1, 0);
1379 
1380 #pragma omp barrier
1381  }
1382 
1383 
1384 //====================================================================
1385  void Fopr_Wilson::clear(Field& w)
1386  {
1387  double *wp = w.ptr(0);
1388 
1389  int ith, nth, is, ns;
1390  set_threadtask(ith, nth, is, ns, m_Nvol);
1391 
1392  int Nvcd = m_Nvc * m_Nd;
1393 
1394 #pragma omp barrier
1395 
1396  for (int site = is; site < ns; ++site) {
1397  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
1398  wp[ivcd + Nvcd * site] = 0.0;
1399  }
1400  }
1401 
1402 #pragma omp barrier
1403  }
1404 
1405 
1406 //====================================================================
1407  void Fopr_Wilson::daypx(Field& v,
1408  const double fac, const Field& w)
1409  {
1410  double *vp = v.ptr(0);
1411  const double *wp = w.ptr(0);
1412 
1413  int ith, nth, is, ns;
1414  set_threadtask(ith, nth, is, ns, m_Nvol);
1415 
1416  int Nvcd = m_Nvc * m_Nd;
1417 
1418 #pragma omp barrier
1419 
1420  for (int site = is; site < ns; ++site) {
1421  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
1422  vp[ivcd + Nvcd * site]
1423  = fac * vp[ivcd + Nvcd * site] + wp[ivcd + Nvcd * site];
1424  }
1425  }
1426 
1427 #pragma omp barrier
1428  }
1429 
1430 
1431 //====================================================================
1432  void Fopr_Wilson::mult_gm5_dirac(Field& v, const Field& w)
1433  {
1434  double *vp = v.ptr(0);
1435  const double *wp = w.ptr(0);
1436 
1437  int ith, nth, is, ns;
1438  set_threadtask(ith, nth, is, ns, m_Nvol);
1439 
1440  int Nvcd = m_Nvc * m_Nd;
1441 
1442 #pragma omp barrier
1443 
1444  for (int site = is; site < ns; ++site) {
1445  mult_gamma5_dirac(&vp[Nvcd * site], &wp[Nvcd * site], m_Nc);
1446  }
1447 
1448 #pragma omp barrier
1449  }
1450 
1451 
1452 //====================================================================
1453  void Fopr_Wilson::mult_gm5_chiral(Field& v, const Field& w)
1454  {
1455  double *vp = v.ptr(0);
1456  const double *wp = w.ptr(0);
1457 
1458  int ith, nth, is, ns;
1459  set_threadtask(ith, nth, is, ns, m_Nvol);
1460 
1461  int Nvcd = m_Nvc * m_Nd;
1462 
1463 #pragma omp barrier
1464 
1465  for (int site = is; site < ns; ++site) {
1466  mult_gamma5_chiral(&vp[Nvcd * site], &wp[Nvcd * site], m_Nc);
1467  }
1468 
1469 #pragma omp barrier
1470  }
1471 
1472 
1473 //====================================================================
1474  void Fopr_Wilson::mult_xp(Field& v, const Field& w)
1475  {
1476  int idir = 0;
1477 
1478  int Nvc2 = m_Nvc * 2;
1479  int Nvcd = m_Nvc * m_Nd;
1480 
1481  double bc2 = m_boundary_each_node[idir];
1482 
1483  double *vp = v.ptr(0);
1484  const double *wp = w.ptr(0);
1485  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1486 
1487  int ith, nth, is, ns;
1488  set_threadtask(ith, nth, is, ns, m_Nvol);
1489 
1490 #pragma omp barrier
1491 
1492  for (int site = is; site < ns; ++site) {
1493  int ix = site % m_Nx;
1494  int iyzt = site / m_Nx;
1495  if (ix == 0) {
1496  int in = Nvcd * site;
1497  int ix1 = Nvc2 * iyzt;
1498  int ix2 = ix1 + NVC;
1499  double vt1[NVC], vt2[NVC];
1500  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
1501  for (int ivc = 0; ivc < NVC; ++ivc) {
1502  vcp1_xp[ivc + ix1] = bc2 * vt1[ivc];
1503  vcp1_xp[ivc + ix2] = bc2 * vt2[ivc];
1504  }
1505  }
1506  }
1507 
1508 #pragma omp barrier
1509 
1510 #pragma omp master
1511  {
1512  const int Nv = m_Nvc * 2 * m_Ny * m_Nz * m_Nt;
1513  Communicator::exchange(Nv, vcp2_xp, vcp1_xp, 0, 1, 1);
1514  }
1515 #pragma omp barrier
1516 
1517  for (int site = is; site < ns; ++site) {
1518  int ix = site % m_Nx;
1519  int iyzt = site / m_Nx;
1520  int nei = ix + 1 + m_Nx * iyzt;
1521  int iv = Nvcd * site;
1522  int ig = m_Ndf * site;
1523 
1524  if (ix < m_Nx - 1) {
1525  int in = Nvcd * nei;
1526  double vt1[NVC], vt2[NVC];
1527  set_sp2_xp(vt1, vt2, &wp[in], m_Nc);
1528  for (int ic = 0; ic < m_Nc; ++ic) {
1529  int ic2 = ic * NVC;
1530  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1531  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1532  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1533  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1534  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1535  }
1536  } else {
1537  int ix1 = Nvc2 * iyzt;
1538  int ix2 = ix1 + NVC;
1539  for (int ic = 0; ic < m_Nc; ++ic) {
1540  int ic2 = ic * NVC;
1541  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
1542  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
1543  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
1544  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
1545  set_sp4_xp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1546  }
1547  }
1548  }
1549 
1550 #pragma omp barrier
1551  }
1552 
1553 
1554 //====================================================================
1555  void Fopr_Wilson::mult_xm(Field& v, const Field& w)
1556  {
1557  int idir = 0;
1558 
1559  int Nvc2 = m_Nvc * 2;
1560  int Nvcd = m_Nvc * m_Nd;
1561 
1562  double bc2 = m_boundary_each_node[idir];
1563 
1564  double *vp = v.ptr(0);
1565  const double *wp = w.ptr(0);
1566  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1567 
1568  int ith, nth, is, ns;
1569  set_threadtask(ith, nth, is, ns, m_Nvol);
1570 
1571 #pragma omp barrier
1572 
1573  for (int site = is; site < ns; ++site) {
1574  int ix = site % m_Nx;
1575  int iyzt = site / m_Nx;
1576  if (ix == m_Nx - 1) {
1577  int in = Nvcd * site;
1578  int ig = m_Ndf * site;
1579  int ix1 = Nvc2 * iyzt;
1580  int ix2 = ix1 + NVC;
1581 
1582  double vt1[NVC], vt2[NVC];
1583  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
1584  for (int ic = 0; ic < m_Nc; ++ic) {
1585  int ic2 = 2 * ic;
1586  int icr = 2 * ic;
1587  int ici = 2 * ic + 1;
1588  vcp1_xm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1589  vcp1_xm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1590  vcp1_xm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1591  vcp1_xm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1592  }
1593  }
1594  }
1595 
1596 #pragma omp barrier
1597 
1598 #pragma omp master
1599  {
1600  const int Nv = m_Nvc * 2 * m_Ny * m_Nz * m_Nt;
1601  Communicator::exchange(Nv, vcp2_xm, vcp1_xm, 0, -1, 2);
1602  }
1603 #pragma omp barrier
1604 
1605  for (int site = is; site < ns; ++site) {
1606  int ix = site % m_Nx;
1607  int iyzt = site / m_Nx;
1608  int nei = ix - 1 + m_Nx * iyzt;
1609  int iv = Nvcd * site;
1610 
1611  if (ix > 0) {
1612  int ig = m_Ndf * nei;
1613  int in = Nvcd * nei;
1614  double vt1[NVC], vt2[NVC];
1615  set_sp2_xm(vt1, vt2, &wp[in], m_Nc);
1616  for (int ic = 0; ic < m_Nc; ++ic) {
1617  int ic2 = 2 * ic;
1618  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1619  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1620  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1621  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1622  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1623  }
1624  } else {
1625  int ix1 = Nvc2 * iyzt;
1626  int ix2 = ix1 + NVC;
1627  for (int ic = 0; ic < m_Nc; ++ic) {
1628  double wt1r = bc2 * vcp2_xm[2 * ic + ix1];
1629  double wt1i = bc2 * vcp2_xm[2 * ic + 1 + ix1];
1630  double wt2r = bc2 * vcp2_xm[2 * ic + ix2];
1631  double wt2i = bc2 * vcp2_xm[2 * ic + 1 + ix2];
1632  set_sp4_xm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1633  }
1634  }
1635  }
1636 
1637 #pragma omp barrier
1638  }
1639 
1640 
1641 //====================================================================
1642  void Fopr_Wilson::mult_yp(Field& v, const Field& w)
1643  {
1644  int idir = 1;
1645 
1646  int Nvc2 = m_Nvc * 2;
1647  int Nvcd = m_Nvc * m_Nd;
1648 
1649  double bc2 = m_boundary_each_node[idir];
1650 
1651  double *vp = v.ptr(0);
1652  const double *wp = w.ptr(0);
1653  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1654 
1655  int ith, nth, is, ns;
1656  set_threadtask(ith, nth, is, ns, m_Nvol);
1657 
1658 #pragma omp barrier
1659 
1660  for (int site = is; site < ns; ++site) {
1661  int ix = site % m_Nx;
1662  int iyzt = site / m_Nx;
1663  int iy = iyzt % m_Ny;
1664  int izt = iyzt / m_Ny;
1665  int ixzt = ix + m_Nx * izt;
1666  if (iy == 0) {
1667  int in = Nvcd * site;
1668  int ix1 = Nvc2 * ixzt;
1669  int ix2 = ix1 + NVC;
1670  double vt1[NVC], vt2[NVC];
1671  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
1672  for (int ivc = 0; ivc < NVC; ++ivc) {
1673  vcp1_yp[ivc + ix1] = bc2 * vt1[ivc];
1674  vcp1_yp[ivc + ix2] = bc2 * vt2[ivc];
1675  }
1676  }
1677  }
1678 
1679 #pragma omp barrier
1680 
1681 #pragma omp master
1682  {
1683  const int Nv = m_Nvc * 2 * m_Nx * m_Nz * m_Nt;
1684  Communicator::exchange(Nv, vcp2_yp, vcp1_yp, 1, 1, 3);
1685  }
1686 #pragma omp barrier
1687 
1688  for (int site = is; site < ns; ++site) {
1689  int ix = site % m_Nx;
1690  int iyzt = site / m_Nx;
1691  int iy = iyzt % m_Ny;
1692  int izt = iyzt / m_Ny;
1693  int ixzt = ix + m_Nx * izt;
1694  int nei = ix + m_Nx * (iy + 1 + m_Ny * izt);
1695  int iv = Nvcd * site;
1696  int ig = m_Ndf * site;
1697 
1698  if (iy < m_Ny - 1) {
1699  int in = Nvcd * nei;
1700  double vt1[NVC], vt2[NVC];
1701  set_sp2_yp(vt1, vt2, &wp[in], m_Nc);
1702  for (int ic = 0; ic < m_Nc; ++ic) {
1703  int ic2 = ic * NVC;
1704  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1705  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1706  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1707  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1708  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1709  }
1710  } else {
1711  int ix1 = Nvc2 * ixzt;
1712  int ix2 = ix1 + NVC;
1713  for (int ic = 0; ic < m_Nc; ++ic) {
1714  int ic2 = ic * NVC;
1715  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1716  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1717  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1718  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1719  set_sp4_yp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1720  }
1721  }
1722  }
1723 
1724 #pragma omp barrier
1725  }
1726 
1727 
1728 //====================================================================
1729  void Fopr_Wilson::mult_ym(Field& v, const Field& w)
1730  {
1731  int idir = 1;
1732 
1733  int Nvc2 = m_Nvc * 2;
1734  int Nvcd = m_Nvc * m_Nd;
1735 
1736  double bc2 = m_boundary_each_node[idir];
1737 
1738  double *vp = v.ptr(0);
1739  const double *wp = w.ptr(0);
1740  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1741 
1742  int ith, nth, is, ns;
1743  set_threadtask(ith, nth, is, ns, m_Nvol);
1744 
1745 #pragma omp barrier
1746 
1747  for (int site = is; site < ns; ++site) {
1748  int ix = site % m_Nx;
1749  int iyzt = site / m_Nx;
1750  int iy = iyzt % m_Ny;
1751  int izt = iyzt / m_Ny;
1752  int ixzt = ix + m_Nx * izt;
1753  if (iy == m_Ny - 1) {
1754  int in = Nvcd * site;
1755  int ig = m_Ndf * site;
1756  int ix1 = Nvc2 * ixzt;
1757  int ix2 = ix1 + NVC;
1758 
1759  double vt1[NVC], vt2[NVC];
1760  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
1761  for (int ic = 0; ic < m_Nc; ++ic) {
1762  int ic2 = 2 * ic;
1763  int icr = 2 * ic;
1764  int ici = 2 * ic + 1;
1765  vcp1_ym[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1766  vcp1_ym[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1767  vcp1_ym[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1768  vcp1_ym[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1769  }
1770  }
1771  }
1772 
1773 #pragma omp barrier
1774 
1775 #pragma omp master
1776  {
1777  const int Nv = m_Nvc * 2 * m_Nx * m_Nz * m_Nt;
1778  Communicator::exchange(Nv, vcp2_ym, vcp1_ym, 1, -1, 4);
1779  }
1780 #pragma omp barrier
1781 
1782  for (int site = is; site < ns; ++site) {
1783  int ix = site % m_Nx;
1784  int iyzt = site / m_Nx;
1785  int iy = iyzt % m_Ny;
1786  int izt = iyzt / m_Ny;
1787  int ixzt = ix + m_Nx * izt;
1788  int nei = ix + m_Nx * (iy - 1 + m_Ny * izt);
1789  int iv = Nvcd * site;
1790 
1791  if (iy > 0) {
1792  int ig = m_Ndf * nei;
1793  int in = Nvcd * nei;
1794  double vt1[NVC], vt2[NVC];
1795  set_sp2_ym(vt1, vt2, &wp[in], m_Nc);
1796  for (int ic = 0; ic < m_Nc; ++ic) {
1797  int ic2 = 2 * ic;
1798  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1799  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1800  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1801  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1802  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1803  }
1804  } else {
1805  int ix1 = Nvc2 * ixzt;
1806  int ix2 = ix1 + NVC;
1807  for (int ic = 0; ic < m_Nc; ++ic) {
1808  double wt1r = bc2 * vcp2_ym[2 * ic + ix1];
1809  double wt1i = bc2 * vcp2_ym[2 * ic + 1 + ix1];
1810  double wt2r = bc2 * vcp2_ym[2 * ic + ix2];
1811  double wt2i = bc2 * vcp2_ym[2 * ic + 1 + ix2];
1812  set_sp4_ym(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1813  }
1814  }
1815  }
1816 
1817 #pragma omp barrier
1818  }
1819 
1820 
1821 //====================================================================
1822  void Fopr_Wilson::mult_zp(Field& v, const Field& w)
1823  {
1824  int idir = 2;
1825 
1826  int Nvc2 = m_Nvc * 2;
1827  int Nvcd = m_Nvc * m_Nd;
1828 
1829  double bc2 = m_boundary_each_node[idir];
1830 
1831  double *vp = v.ptr(0);
1832  const double *wp = w.ptr(0);
1833  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1834 
1835  int ith, nth, is, ns;
1836  set_threadtask(ith, nth, is, ns, m_Nvol);
1837 
1838  int Nxy = m_Nx * m_Ny;
1839 
1840 #pragma omp barrier
1841 
1842  for (int site = is; site < ns; ++site) {
1843  int ixy = site % Nxy;
1844  int izt = site / Nxy;
1845  int iz = izt % m_Nz;
1846  int it = izt / m_Nz;
1847  int ixyt = ixy + Nxy * it;
1848  if (iz == 0) {
1849  int in = Nvcd * site;
1850  int ix1 = Nvc2 * ixyt;
1851  int ix2 = ix1 + NVC;
1852  double vt1[NVC], vt2[NVC];
1853  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
1854  for (int ivc = 0; ivc < NVC; ++ivc) {
1855  vcp1_zp[ivc + ix1] = bc2 * vt1[ivc];
1856  vcp1_zp[ivc + ix2] = bc2 * vt2[ivc];
1857  }
1858  }
1859  }
1860 
1861 #pragma omp barrier
1862 
1863 #pragma omp master
1864  {
1865  int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nt;
1866  Communicator::exchange(Nv, vcp2_zp, vcp1_zp, 2, 1, 5);
1867  }
1868 #pragma omp barrier
1869 
1870  for (int site = is; site < ns; ++site) {
1871  int ixy = site % Nxy;
1872  int izt = site / Nxy;
1873  int iz = izt % m_Nz;
1874  int it = izt / m_Nz;
1875  int ixyt = ixy + Nxy * it;
1876  int nei = ixy + Nxy * (iz + 1 + m_Nz * it);
1877  int iv = Nvcd * site;
1878  int ig = m_Ndf * site;
1879 
1880  if (iz < m_Nz - 1) {
1881  int in = Nvcd * nei;
1882  double vt1[NVC], vt2[NVC];
1883  set_sp2_zp(vt1, vt2, &wp[in], m_Nc);
1884  for (int ic = 0; ic < m_Nc; ++ic) {
1885  int ic2 = ic * NVC;
1886  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1887  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1888  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1889  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1890  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1891  }
1892  } else {
1893  int ix1 = Nvc2 * ixyt;
1894  int ix2 = ix1 + NVC;
1895  for (int ic = 0; ic < m_Nc; ++ic) {
1896  int ic2 = ic * NVC;
1897  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1898  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1899  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1900  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1901  set_sp4_zp(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1902  }
1903  }
1904  }
1905 
1906 #pragma omp barrier
1907  }
1908 
1909 
1910 //====================================================================
1911  void Fopr_Wilson::mult_zm(Field& v, const Field& w)
1912  {
1913  int idir = 2;
1914 
1915  int Nvc2 = m_Nvc * 2;
1916  int Nvcd = m_Nvc * m_Nd;
1917 
1918  double bc2 = m_boundary_each_node[idir];
1919 
1920  double *vp = v.ptr(0);
1921  const double *wp = w.ptr(0);
1922  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1923 
1924  int ith, nth, is, ns;
1925  set_threadtask(ith, nth, is, ns, m_Nvol);
1926 
1927  int Nxy = m_Nx * m_Ny;
1928 
1929 #pragma omp barrier
1930 
1931  for (int site = is; site < ns; ++site) {
1932  int ixy = site % Nxy;
1933  int izt = site / Nxy;
1934  int iz = izt % m_Nz;
1935  int it = izt / m_Nz;
1936  int ixyt = ixy + Nxy * it;
1937  if (iz == m_Nz - 1) {
1938  int in = Nvcd * site;
1939  int ig = m_Ndf * site;
1940  int ix1 = Nvc2 * ixyt;
1941  int ix2 = ix1 + NVC;
1942 
1943  double vt1[NVC], vt2[NVC];
1944  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
1945  for (int ic = 0; ic < m_Nc; ++ic) {
1946  int ic2 = 2 * ic;
1947  int icr = 2 * ic;
1948  int ici = 2 * ic + 1;
1949  vcp1_zm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1950  vcp1_zm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1951  vcp1_zm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1952  vcp1_zm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1953  }
1954  }
1955  }
1956 
1957 #pragma omp barrier
1958 
1959 #pragma omp master
1960  {
1961  const int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nt;
1962  Communicator::exchange(Nv, vcp2_zm, vcp1_zm, 2, -1, 6);
1963  }
1964 #pragma omp barrier
1965 
1966  for (int site = is; site < ns; ++site) {
1967  int ixy = site % Nxy;
1968  int izt = site / Nxy;
1969  int iz = izt % m_Nz;
1970  int it = izt / m_Nz;
1971  int ixyt = ixy + Nxy * it;
1972  int nei = ixy + Nxy * (iz - 1 + m_Nz * it);
1973  int iv = Nvcd * site;
1974 
1975  if (iz > 0) {
1976  int ig = m_Ndf * nei;
1977  int in = Nvcd * nei;
1978  double vt1[NVC], vt2[NVC];
1979  set_sp2_zm(vt1, vt2, &wp[in], m_Nc);
1980  for (int ic = 0; ic < m_Nc; ++ic) {
1981  int ic2 = 2 * ic;
1982  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1983  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1984  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1985  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1986  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1987  }
1988  } else {
1989  int ix1 = Nvc2 * ixyt;
1990  int ix2 = ix1 + NVC;
1991  for (int ic = 0; ic < m_Nc; ++ic) {
1992  double wt1r = bc2 * vcp2_zm[2 * ic + ix1];
1993  double wt1i = bc2 * vcp2_zm[2 * ic + 1 + ix1];
1994  double wt2r = bc2 * vcp2_zm[2 * ic + ix2];
1995  double wt2i = bc2 * vcp2_zm[2 * ic + 1 + ix2];
1996  set_sp4_zm(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1997  }
1998  }
1999  }
2000 
2001 #pragma omp barrier
2002  }
2003 
2004 
2005 //====================================================================
2006  void Fopr_Wilson::mult_tp_dirac(Field& v, const Field& w)
2007  {
2008  int idir = 3;
2009 
2010  int Nvc2 = m_Nvc * 2;
2011  int Nvcd = m_Nvc * m_Nd;
2012 
2013  double bc2 = m_boundary_each_node[idir];
2014 
2015  double *vp = v.ptr(0);
2016  const double *wp = w.ptr(0);
2017  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
2018 
2019  int ith, nth, is, ns;
2020  set_threadtask(ith, nth, is, ns, m_Nvol);
2021 
2022  int Nxyz = m_Nx * m_Ny * m_Nz;
2023 
2024 #pragma omp barrier
2025 
2026  for (int site = is; site < ns; ++site) {
2027  int ixyz = site % Nxyz;
2028  int it = site / Nxyz;
2029  if (it == 0) {
2030  int in = Nvcd * site;
2031  int ix1 = Nvc2 * ixyz;
2032  int ix2 = ix1 + NVC;
2033  double vt1[NVC], vt2[NVC];
2034  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
2035  for (int ivc = 0; ivc < NVC; ++ivc) {
2036  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
2037  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
2038  }
2039  }
2040  }
2041 
2042 #pragma omp barrier
2043 
2044 #pragma omp master
2045  {
2046  int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
2047  Communicator::exchange(Nv, vcp2_tp, vcp1_tp, 3, 1, 7);
2048  }
2049 #pragma omp barrier
2050 
2051  for (int site = is; site < ns; ++site) {
2052  int ixyz = site % Nxyz;
2053  int it = site / Nxyz;
2054  int nei = ixyz + Nxyz * (it + 1);
2055  int iv = Nvcd * site;
2056  int ig = m_Ndf * site;
2057 
2058  if (it < m_Nt - 1) {
2059  int in = Nvcd * nei;
2060  double vt1[NVC], vt2[NVC];
2061  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
2062  for (int ic = 0; ic < m_Nc; ++ic) {
2063  int ic2 = ic * NVC;
2064  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
2065  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
2066  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
2067  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
2068  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2069  }
2070  } else {
2071  int ix1 = Nvc2 * ixyz;
2072  int ix2 = ix1 + NVC;
2073  for (int ic = 0; ic < m_Nc; ++ic) {
2074  int ic2 = ic * NVC;
2075  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
2076  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
2077  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
2078  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
2079  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2080  }
2081  }
2082  }
2083 
2084 #pragma omp barrier
2085  }
2086 
2087 
2088 //====================================================================
2089  void Fopr_Wilson::mult_tm_dirac(Field& v, const Field& w)
2090  {
2091  int idir = 3;
2092 
2093  int Nvc2 = m_Nvc * 2;
2094  int Nvcd = m_Nvc * m_Nd;
2095 
2096  double bc2 = m_boundary_each_node[idir];
2097 
2098  double *vp = v.ptr(0);
2099  const double *wp = w.ptr(0);
2100  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
2101 
2102  int ith, nth, is, ns;
2103  set_threadtask(ith, nth, is, ns, m_Nvol);
2104 
2105  int Nxyz = m_Nx * m_Ny * m_Nz;
2106 
2107 #pragma omp barrier
2108 
2109  for (int site = is; site < ns; ++site) {
2110  int ixyz = site % Nxyz;
2111  int it = site / Nxyz;
2112  if (it == m_Nt - 1) {
2113  int in = Nvcd * site;
2114  int ig = m_Ndf * site;
2115  int ix1 = Nvc2 * ixyz;
2116  int ix2 = ix1 + NVC;
2117 
2118  double vt1[NVC], vt2[NVC];
2119  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
2120  for (int ic = 0; ic < m_Nc; ++ic) {
2121  int ic2 = 2 * ic;
2122  int icr = 2 * ic;
2123  int ici = 2 * ic + 1;
2124  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
2125  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
2126  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
2127  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
2128  }
2129  }
2130  }
2131 
2132 #pragma omp barrier
2133 
2134 #pragma omp master
2135  {
2136  const int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
2137  Communicator::exchange(Nv, vcp2_tm, vcp1_tm, 3, -1, 8);
2138  }
2139 #pragma omp barrier
2140 
2141  for (int site = is; site < ns; ++site) {
2142  int ixyz = site % Nxyz;
2143  int it = site / Nxyz;
2144  int nei = ixyz + Nxyz * (it - 1);
2145  int iv = Nvcd * site;
2146 
2147  if (it > 0) {
2148  int ig = m_Ndf * nei;
2149  int in = Nvcd * nei;
2150  double vt1[NVC], vt2[NVC];
2151  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
2152  for (int ic = 0; ic < m_Nc; ++ic) {
2153  int ic2 = 2 * ic;
2154  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
2155  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
2156  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
2157  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
2158  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2159  }
2160  } else {
2161  int ix1 = Nvc2 * ixyz;
2162  int ix2 = ix1 + NVC;
2163  for (int ic = 0; ic < m_Nc; ++ic) {
2164  int icr = 2 * ic;
2165  int ici = 2 * ic + 1;
2166  double wt1r = bc2 * vcp2_tm[icr + ix1];
2167  double wt1i = bc2 * vcp2_tm[ici + ix1];
2168  double wt2r = bc2 * vcp2_tm[icr + ix2];
2169  double wt2i = bc2 * vcp2_tm[ici + ix2];
2170  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2171  }
2172  }
2173  }
2174 
2175 #pragma omp barrier
2176  }
2177 
2178 
2179 //====================================================================
2180  void Fopr_Wilson::mult_tp_chiral(Field& v, const Field& w)
2181  {
2182  int idir = 3;
2183 
2184  int Nvc2 = m_Nvc * 2;
2185  int Nvcd = m_Nvc * m_Nd;
2186 
2187  double bc2 = m_boundary_each_node[idir];
2188 
2189  double *vp = v.ptr(0);
2190  const double *wp = w.ptr(0);
2191  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
2192 
2193  int ith, nth, is, ns;
2194  set_threadtask(ith, nth, is, ns, m_Nvol);
2195 
2196  int Nxyz = m_Nx * m_Ny * m_Nz;
2197 
2198 #pragma omp barrier
2199 
2200  for (int site = is; site < ns; ++site) {
2201  int ixyz = site % Nxyz;
2202  int it = site / Nxyz;
2203  if (it == 0) {
2204  int in = Nvcd * site;
2205  int ix1 = Nvc2 * ixyz;
2206  int ix2 = ix1 + NVC;
2207  double vt1[NVC], vt2[NVC];
2208  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
2209  for (int ivc = 0; ivc < NVC; ++ivc) {
2210  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
2211  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
2212  }
2213  }
2214  }
2215 
2216 #pragma omp barrier
2217 
2218 #pragma omp master
2219  {
2220  int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
2221  Communicator::exchange(Nv, vcp2_tp, vcp1_tp, 3, 1, 7);
2222  }
2223 #pragma omp barrier
2224 
2225  for (int site = is; site < ns; ++site) {
2226  int ixyz = site % Nxyz;
2227  int it = site / Nxyz;
2228  int nei = ixyz + Nxyz * (it + 1);
2229  int iv = Nvcd * site;
2230  int ig = m_Ndf * site;
2231 
2232  if (it < m_Nt - 1) {
2233  int in = Nvcd * nei;
2234  double vt1[NVC], vt2[NVC];
2235  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
2236  for (int ic = 0; ic < m_Nc; ++ic) {
2237  int ic2 = ic * NVC;
2238  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
2239  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
2240  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
2241  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
2242  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2243  }
2244  } else {
2245  int ix1 = Nvc2 * ixyz;
2246  int ix2 = ix1 + NVC;
2247  for (int ic = 0; ic < m_Nc; ++ic) {
2248  int ic2 = ic * NVC;
2249  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
2250  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
2251  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
2252  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
2253  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2254  }
2255  }
2256  }
2257 
2258 #pragma omp barrier
2259  }
2260 
2261 
2262 //====================================================================
2263  void Fopr_Wilson::mult_tm_chiral(Field& v, const Field& w)
2264  {
2265  int idir = 3;
2266 
2267  int Nvc2 = m_Nvc * 2;
2268  int Nvcd = m_Nvc * m_Nd;
2269 
2270  double bc2 = m_boundary_each_node[idir];
2271 
2272  double *vp = v.ptr(0);
2273  const double *wp = w.ptr(0);
2274  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
2275 
2276  int ith, nth, is, ns;
2277  set_threadtask(ith, nth, is, ns, m_Nvol);
2278 
2279  int Nxyz = m_Nx * m_Ny * m_Nz;
2280 
2281 #pragma omp barrier
2282 
2283  for (int site = is; site < ns; ++site) {
2284  int ixyz = site % Nxyz;
2285  int it = site / Nxyz;
2286  if (it == m_Nt - 1) {
2287  int in = Nvcd * site;
2288  int ig = m_Ndf * site;
2289  int ix1 = Nvc2 * ixyz;
2290  int ix2 = ix1 + NVC;
2291 
2292  double vt1[NVC], vt2[NVC];
2293  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
2294  for (int ic = 0; ic < m_Nc; ++ic) {
2295  int ic2 = 2 * ic;
2296  int icr = 2 * ic;
2297  int ici = 2 * ic + 1;
2298  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
2299  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
2300  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
2301  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
2302  }
2303  }
2304  }
2305 
2306 #pragma omp barrier
2307 
2308 #pragma omp master
2309  {
2310  const int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
2311  Communicator::exchange(Nv, vcp2_tm, vcp1_tm, 3, -1, 8);
2312  }
2313 #pragma omp barrier
2314 
2315  for (int site = is; site < ns; ++site) {
2316  int ixyz = site % Nxyz;
2317  int it = site / Nxyz;
2318  int nei = ixyz + Nxyz * (it - 1);
2319  int iv = Nvcd * site;
2320 
2321  if (it > 0) {
2322  int ig = m_Ndf * nei;
2323  int in = Nvcd * nei;
2324  double vt1[NVC], vt2[NVC];
2325  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
2326  for (int ic = 0; ic < m_Nc; ++ic) {
2327  int ic2 = 2 * ic;
2328  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
2329  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
2330  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
2331  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
2332  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2333  }
2334  } else {
2335  int ix1 = Nvc2 * ixyz;
2336  int ix2 = ix1 + NVC;
2337  for (int ic = 0; ic < m_Nc; ++ic) {
2338  double wt1r = bc2 * vcp2_tm[2 * ic + ix1];
2339  double wt1i = bc2 * vcp2_tm[2 * ic + 1 + ix1];
2340  double wt2r = bc2 * vcp2_tm[2 * ic + ix2];
2341  double wt2i = bc2 * vcp2_tm[2 * ic + 1 + ix2];
2342  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
2343  }
2344  }
2345  }
2346 
2347 #pragma omp barrier
2348  }
2349 
2350 
2351 //====================================================================
2352  double Fopr_Wilson::flop_count()
2353  {
2354  // Counting of floating point operations in giga unit.
2355  // The following counting explicitly depends on the implementation.
2356  // It will be recalculated when the code is modified.
2357  // The present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]
2358 
2359  const int Nvol = CommonParameters::Nvol();
2360  const int NPE = CommonParameters::NPE();
2361 
2362  int flop_site;
2363 
2364  if (m_repr == "Dirac") {
2365  flop_site = m_Nc * m_Nd * (4 + 6 * (4 * m_Nc + 2) + 2 * (4 * m_Nc + 1));
2366  } else if (m_repr == "Chiral") {
2367  flop_site = m_Nc * m_Nd * (4 + 8 * (4 * m_Nc + 2));
2368  } else {
2369  vout.crucial(m_vl, "Error at %s: input repr is undefined.\n",
2370  class_name.c_str());
2371  exit(EXIT_FAILURE);
2372  }
2373 
2374  double flop = flop_site * (Nvol * NPE);
2375 
2376  if ((m_mode == "DdagD") || (m_mode == "DDdag")) flop *= 2;
2377 
2378  double gflop = flop * 1.e-9;
2379  return gflop;
2380  }
2381 
2382 
2383 //====================================================================
2384 }
2385 //============================================================END=====
Imp::Fopr_Wilson::mult_ym
void mult_ym(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1729
Imp::Fopr_Wilson::mult_tp_dirac
void mult_tp_dirac(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:2006
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
Imp::Fopr_Wilson::class_name
static const std::string class_name
Definition: fopr_Wilson_impl.h:46
Imp::Fopr_Wilson::m_w2
Field m_w2
working fields
Definition: fopr_Wilson_impl.h:72
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
Imp::Fopr_Wilson::Ddag
void Ddag(Field &v, const Field &w)
Definition: fopr_Wilson_impl.cpp:433
Imp::Fopr_Wilson::flop_count
double flop_count()
returns the number of floating point operations.
Definition: fopr_Wilson_impl.cpp:2352
Imp::Fopr_Wilson::vcp1_zp
double * vcp1_zp
Definition: fopr_Wilson_impl.h:69
fopr_thread-inc.h
Imp::Fopr_Wilson::vcp1_xm
double * vcp1_xm
Definition: fopr_Wilson_impl.h:67
Imp::Fopr_Wilson::mult_gm5
void mult_gm5(Field &v, const Field &w)
multiplies gamma_5 matrix.
Definition: fopr_Wilson_impl.cpp:397
Imp::Fopr_Wilson::D_ex_dirac
void D_ex_dirac(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_Wilson_impl.cpp:470
Imp::Fopr_Wilson::mult
void mult(Field &v, const Field &w)
multiplies fermion operator to a given field.
Definition: fopr_Wilson_impl.cpp:265
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Imp::Fopr_Wilson::m_Ny
int m_Ny
Definition: fopr_Wilson_impl.h:59
fopr_Wilson_impl_SU_N-inc.h
fopr_Wilson_impl_SU3-inc.h
Imp::Fopr_Wilson::m_mode
std::string m_mode
mult mode
Definition: fopr_Wilson_impl.h:55
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Parameters
Class for parameters.
Definition: parameters.h:46
iy
int iy
Definition: mult_Wilson_xyz_openacc-inc.h:239
ixy
int ixy
Definition: mult_Domainwall_eo_xyz_openacc-inc.h:513
Imp::Fopr_Wilson::mult_zm
void mult_zm(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1911
Imp::Fopr_Wilson::vcp2_tp
double * vcp2_tp
Definition: fopr_Wilson_impl.h:70
izt
int izt
Definition: mult_Domainwall_eo_xyz_openacc-inc.h:262
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Imp::Fopr_Wilson::proj_chiral
void proj_chiral(Field &w, const int ex1, const Field &v, const int ex2, const int ipm)
Definition: fopr_Wilson_impl.cpp:1361
Imp::Fopr_Wilson::m_Nd
int m_Nd
Definition: fopr_Wilson_impl.h:58
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:518
Imp::Fopr_Wilson::mult_tm_dirac
void mult_tm_dirac(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:2089
Imp::Fopr_Wilson::mult_xp
void mult_xp(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1474
Imp::Fopr_Wilson::DdagD
void DdagD(Field &v, const Field &w)
Definition: fopr_Wilson_impl.cpp:442
Imp::Fopr_Wilson::DDdag
void DDdag(Field &v, const Field &w)
Definition: fopr_Wilson_impl.cpp:452
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Imp::Fopr_Wilson::vcp2_xm
double * vcp2_xm
Definition: fopr_Wilson_impl.h:67
Imp::Fopr_Wilson::mult_dn
void mult_dn(const int mu, Field &, const Field &)
downward nearest neighbor hopping term.
Definition: fopr_Wilson_impl.cpp:374
Imp::Fopr_Wilson::mult_xm
void mult_xm(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1555
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Imp::Fopr_Wilson::vcp2_tm
double * vcp2_tm
Definition: fopr_Wilson_impl.h:70
wt1i
wt1i
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:55
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
wt2r
wt2r
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:57
Imp::Fopr_Wilson::m_Ndf
int m_Ndf
Definition: fopr_Wilson_impl.h:58
Imp::Fopr_Wilson::vcp2_zp
double * vcp2_zp
Definition: fopr_Wilson_impl.h:69
Imp::Fopr_Wilson::mult_zp
void mult_zp(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1822
Imp::Fopr_Wilson::m_w1
Field m_w1
Definition: fopr_Wilson_impl.h:72
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
Imp::Fopr_Wilson::m_kappa
double m_kappa
hopping parameter
Definition: fopr_Wilson_impl.h:50
Imp::Fopr_Wilson::mult_gm5p
void mult_gm5p(const int mu, Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1352
iyzt
int iyzt
Definition: mult_Domainwall_eo_xyz_openacc-inc.h:14
Imp::Fopr_Wilson::m_Nz
int m_Nz
Definition: fopr_Wilson_impl.h:59
Imp::Fopr_Wilson::mult_yp
void mult_yp(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1642
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
ix
int ix
Definition: mult_Wilson_xyz_openacc-inc.h:16
Imp::Fopr_Wilson::D
void D(Field &v, const Field &w)
Definition: fopr_Wilson_impl.cpp:408
idir
idir
Definition: mult_Domainwall_eo_xyz_openacc-inc.h:264
Imp::Fopr_Wilson::setup
void setup()
initial setup main.
Definition: fopr_Wilson_impl.cpp:96
Imp::Fopr_Wilson::m_boundary_each_node
std::vector< double > m_boundary_each_node
b.c. on each node.
Definition: fopr_Wilson_impl.h:64
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Imp::Fopr_Wilson::mult_up
void mult_up(const int mu, Field &, const Field &)
upward nearest neighbor hopping term.
Definition: fopr_Wilson_impl.cpp:351
vt1
real_t vt1
Definition: mult_Staggered_uvdn1_openacc-inc.h:8
Imp::Fopr_Wilson::D_ex
void D_ex(Field &v, const int ex1, const Field &f, const int ex2)
Definition: fopr_Wilson_impl.cpp:421
Imp::Fopr_Wilson::tidyup
void tidyup()
final clean-up.
Definition: fopr_Wilson_impl.cpp:150
Imp::Fopr_Wilson::mult_gm5_chiral
void mult_gm5_chiral(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1453
bc2
int bc2
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:62
Imp::Fopr_Wilson::mult_gm5_dirac
void mult_gm5_dirac(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:1432
wt1r
wt1r
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:53
Imp::Fopr_Wilson::D_ex_dirac_alt
void D_ex_dirac_alt(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_Wilson_impl.cpp:1318
wt2i
wt2i
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:59
Imp::Fopr_Wilson::D_ex_chiral_alt
void D_ex_chiral_alt(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_Wilson_impl.cpp:1335
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
Imp::Fopr_Wilson::vcp2_yp
double * vcp2_yp
Definition: fopr_Wilson_impl.h:68
Imp::Fopr_Wilson::vcp2_ym
double * vcp2_ym
Definition: fopr_Wilson_impl.h:68
Imp::Fopr_Wilson::vcp1_xp
double * vcp1_xp
arrays for communication buffer.
Definition: fopr_Wilson_impl.h:67
Imp::Fopr_Wilson::m_U
const Field_G * m_U
gauge configuration.
Definition: fopr_Wilson_impl.h:62
Imp::Fopr_Wilson::vcp2_zm
double * vcp2_zm
Definition: fopr_Wilson_impl.h:69
Imp::Fopr_Wilson::vcp2_xp
double * vcp2_xp
Definition: fopr_Wilson_impl.h:67
Imp::Fopr_Wilson::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: fopr_Wilson_impl.cpp:246
Imp::Fopr_Wilson::vcp1_zm
double * vcp1_zm
Definition: fopr_Wilson_impl.h:69
Imp::Fopr_Wilson::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_Wilson_impl.cpp:175
fopr_Wilson_impl_common-inc.h
Imp::Fopr_Wilson::vcp1_ym
double * vcp1_ym
Definition: fopr_Wilson_impl.h:68
it
int it
Definition: mult_Wilson_xyz_openacc-inc.h:461
ix2
int ix2
Definition: mult_Domainwall_eo_xyz_openacc-inc.h:138
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Imp::Fopr_Wilson::D_ex_chiral
void D_ex_chiral(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_Wilson_impl.cpp:894
Field::reset
void reset(const int Nin, const int Nvol, const int Nex, const element_type cmpl=Element_type::COMPLEX)
Definition: field.h:95
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
Imp::Fopr_Wilson::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_Wilson_impl.cpp:235
fopr_Wilson_impl.h
Imp::Fopr_Wilson::daypx
void daypx(Field &, const double, const Field &)
Definition: fopr_Wilson_impl.cpp:1407
Imp::Fopr_Wilson::m_Nx
int m_Nx
Definition: fopr_Wilson_impl.h:59
Imp::Fopr_Wilson::m_boundary
std::vector< int > m_boundary
boundary condition
Definition: fopr_Wilson_impl.h:51
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
NVC
#define NVC
Definition: fopr_Wilson_impl_SU2-inc.h:15
Imp::Fopr_Wilson::vcp1_tm
double * vcp1_tm
Definition: fopr_Wilson_impl.h:70
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Imp::Fopr_Wilson::H
void H(Field &v, const Field &w)
Definition: fopr_Wilson_impl.cpp:462
fopr_Wilson_impl_SU2-inc.h
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
Imp::Fopr_Wilson::m_vl
Bridge::VerboseLevel m_vl
verbose level
Definition: fopr_Wilson_impl.h:53
Imp::Fopr_Wilson::m_repr
std::string m_repr
gamma-matrix representation
Definition: fopr_Wilson_impl.h:52
Imp::Fopr_Wilson::set_mode
void set_mode(const std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: fopr_Wilson_impl.cpp:255
Imp::Fopr_Wilson::vcp1_tp
double * vcp1_tp
Definition: fopr_Wilson_impl.h:70
Imp::Fopr_Wilson::vcp1_yp
double * vcp1_yp
Definition: fopr_Wilson_impl.h:68
iz
int iz
Definition: mult_Wilson_xyz_openacc-inc.h:462
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
ixyz
int ixyz
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:13
Imp::Fopr_Wilson::m_Nc
int m_Nc
Definition: fopr_Wilson_impl.h:58
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
Imp
Clover term operator.
Definition: fopr_CloverTerm_eo_impl.cpp:32
Imp::Fopr_Wilson::m_Nvc
int m_Nvc
Definition: fopr_Wilson_impl.h:58
Imp::Fopr_Wilson::mult_dag
void mult_dag(Field &v, const Field &w)
hermitian conjugate of mult.
Definition: fopr_Wilson_impl.cpp:286
Field
Container of Field-type object.
Definition: field.h:46
Communicator::exchange
static int exchange(int count, dcomplex *recv_buf, dcomplex *send_buf, int idir, int ipm, int tag)
receive array of dcomplex from upstream specified by idir and ipm, and send array to downstream.
Definition: communicator.cpp:207
Imp::Fopr_Wilson::m_Nvol
int m_Nvol
Definition: fopr_Wilson_impl.h:60
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
vt2
real_t vt2
Definition: mult_Staggered_uvdn1_openacc-inc.h:8
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Imp::Fopr_Wilson::clear
void clear(Field &)
Definition: fopr_Wilson_impl.cpp:1385
Imp::Fopr_Wilson::mult_tp_chiral
void mult_tp_chiral(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:2180
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
Imp::Fopr_Wilson::m_Nt
int m_Nt
Definition: fopr_Wilson_impl.h:59
bridge_defs.h
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Imp::Fopr_Wilson::mult_tm_chiral
void mult_tm_chiral(Field &, const Field &)
Definition: fopr_Wilson_impl.cpp:2263
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:572
Imp::Fopr_Wilson::init
void init()
to be discarded.
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:216
Imp::Fopr_Wilson::m_Ndim
int m_Ndim
Definition: fopr_Wilson_impl.h:60