Bridge++  Ver. 2.0.2
fopr_WilsonGeneral_impl.cpp
Go to the documentation of this file.
1 
14 #include "fopr_WilsonGeneral_impl.h"
15 
16 #include "Fopr/fopr_thread-inc.h"
17 
18 #if defined USE_GROUP_SU3
19 #include "fopr_Wilson_impl_SU3-inc.h"
20 #elif defined USE_GROUP_SU2
21 #include "fopr_Wilson_impl_SU2-inc.h"
22 #elif defined USE_GROUP_SU_N
23 #include "fopr_Wilson_impl_SU_N-inc.h"
24 #endif
25 
26 #include "Fopr/Imp/fopr_Wilson_impl_common-inc.h"
27 
28 namespace Imp {
29 #ifdef USE_FACTORY_AUTOREGISTER
30  namespace {
31  bool init = Fopr_WilsonGeneral::register_factory();
32  }
33 #endif
34 
35  const std::string Fopr_WilsonGeneral::class_name
36  = "Imp::Fopr_WilsonGeneral";
37 
38 //====================================================================
39  void Fopr_WilsonGeneral::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_WilsonGeneral::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_WilsonGeneral::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 * m_Nd * 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 * m_Nd * 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 * m_Nd * 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_WilsonGeneral::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_WilsonGeneral::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_s, kappa_t;
187  double nu_s, r_s;
188  std::vector<int> bc;
189 
190  int err = 0;
191  err += params.fetch_double("hopping_parameter_spatial", kappa_s);
192  err += params.fetch_double("hopping_parameter_temporal", kappa_t);
193  err += params.fetch_double("dispersion_parameter_spatial", nu_s);
194  err += params.fetch_double("Wilson_parameter_spatial", r_s);
195  err += params.fetch_int_vector("boundary_condition", bc);
196 
197  if (err) {
198  vout.crucial("Error at %s: input parameter not found.\n",
199  class_name.c_str());
200  exit(EXIT_FAILURE);
201  }
202 
203  set_parameters(kappa_s, kappa_t, nu_s, r_s, bc);
204  }
205 
206 
207 //====================================================================
208  void Fopr_WilsonGeneral::set_parameters(const double kappa_s,
209  const double kappa_t,
210  const double nu_s,
211  const double r_s,
212  const std::vector<int> bc)
213  {
214  assert(bc.size() == m_Ndim);
215 
216 #pragma omp barrier
217 
218  int ith = ThreadManager::get_thread_id();
219  if (ith == 0) {
220  m_kappa_s = kappa_s;
221  m_kappa_t = kappa_t;
222  m_nu_s = nu_s;
223  m_r_s = r_s;
224 
225  m_boundary = bc;
226 
227  for (int idir = 0; idir < m_Ndim; ++idir) {
228  m_boundary_each_node[idir] = 1.0;
229  if (Communicator::ipe(idir) == 0) {
230  m_boundary_each_node[idir] = m_boundary[idir];
231  }
232  }
233  }
234 #pragma omp barrier
235 
236  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
237  vout.general(m_vl, " gamma-matrix type = %s\n", m_repr.c_str());
238  vout.general(m_vl, " kappa_s = %12.8f\n", m_kappa_s);
239  vout.general(m_vl, " kappa_t = %12.8f\n", m_kappa_t);
240  vout.general(m_vl, " nu_s = %12.8f\n", m_nu_s);
241  vout.general(m_vl, " r_s = %12.8f\n", m_r_s);
242  for (int mu = 0; mu < m_Ndim; ++mu) {
243  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
244  }
245  }
246 
247 
248 //====================================================================
249  void Fopr_WilsonGeneral::get_parameters(Parameters& params) const
250  {
251  params.set_double("hopping_parameter_spatial", m_kappa_s);
252  params.set_double("hopping_parameter_temporal", m_kappa_t);
253  params.set_double("dispersion_parameter_spatial", m_nu_s);
254  params.set_double("Wilson_parameter_spatial", m_r_s);
255  params.set_int_vector("boundary_condition", m_boundary);
256  params.set_string("gamma_matrix_type", m_repr);
257 
258  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
259  }
260 
261 
262 //====================================================================
263  void Fopr_WilsonGeneral::set_config(Field *U)
264  {
265 #pragma omp barrier
266  m_U = (Field_G *)U;
267 #pragma omp barrier
268  }
269 
270 
271 //====================================================================
272  void Fopr_WilsonGeneral::set_mode(const std::string mode)
273  {
274 #pragma omp barrier
275  int ith = ThreadManager::get_thread_id();
276  if (ith == 0) m_mode = mode;
277 #pragma omp barrier
278  }
279 
280 
281 //====================================================================
282  void Fopr_WilsonGeneral::mult(Field& v, const Field& w)
283  {
284  if (m_mode == "D") {
285  D(v, w);
286  } else if (m_mode == "Ddag") {
287  Ddag(v, w);
288  } else if (m_mode == "DdagD") {
289  DdagD(v, w);
290  } else if (m_mode == "DDdag") {
291  DDdag(v, w);
292  } else if (m_mode == "H") {
293  H(v, w);
294  } else {
295  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
296  class_name.c_str(), m_mode.c_str());
297  exit(EXIT_FAILURE);
298  }
299  }
300 
301 
302 //====================================================================
303  void Fopr_WilsonGeneral::mult_dag(Field& v, const Field& w)
304  {
305  if (m_mode == "D") {
306  Ddag(v, w);
307  } else if (m_mode == "Ddag") {
308  D(v, w);
309  } else if (m_mode == "DdagD") {
310  DdagD(v, w);
311  } else if (m_mode == "DDdag") {
312  DDdag(v, w);
313  } else if (m_mode == "H") {
314  H(v, w);
315  } else {
316  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
317  class_name.c_str(), m_mode.c_str());
318  exit(EXIT_FAILURE);
319  }
320  }
321 
322 
323 //====================================================================
324  void Fopr_WilsonGeneral::mult(Field& v, const Field& w,
325  const std::string mode)
326  {
327  if (mode == "D") {
328  D(v, w);
329  } else if (mode == "Ddag") {
330  Ddag(v, w);
331  } else if (mode == "DdagD") {
332  DdagD(v, w);
333  } else if (mode == "DDdag") {
334  DDdag(v, w);
335  } else if (mode == "H") {
336  H(v, w);
337  } else {
338  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
339  class_name.c_str(), mode.c_str());
340  exit(EXIT_FAILURE);
341  }
342  }
343 
344 
345 //====================================================================
346  void Fopr_WilsonGeneral::mult_dag(Field& v, const Field& w,
347  const std::string mode)
348  {
349  if (mode == "D") {
350  Ddag(v, w);
351  } else if (mode == "Ddag") {
352  D(v, w);
353  } else if (mode == "DdagD") {
354  DdagD(v, w);
355  } else if (mode == "DDdag") {
356  DDdag(v, w);
357  } else if (mode == "H") {
358  H(v, w);
359  } else {
360  vout.crucial(m_vl, "Error at %s: undefined mode: %s\n",
361  class_name.c_str(), mode.c_str());
362  exit(EXIT_FAILURE);
363  }
364  }
365 
366 
367 //====================================================================
368  void Fopr_WilsonGeneral::mult_up(const int mu, Field& v, const Field& w)
369  {
370  if (mu == 0) {
371  mult_xp(v, w);
372  } else if (mu == 1) {
373  mult_yp(v, w);
374  } else if (mu == 2) {
375  mult_zp(v, w);
376  } else if (mu == 3) {
377  if (m_repr == "Dirac") {
378  mult_tp_dirac(v, w);
379  } else {
380  mult_tp_chiral(v, w);
381  }
382  } else {
383  vout.crucial(m_vl, "Error at %s::mult_up: illegal mu=%d\n",
384  class_name.c_str(), mu);
385  exit(EXIT_FAILURE);
386  }
387  }
388 
389 
390 //====================================================================
391  void Fopr_WilsonGeneral::mult_dn(const int mu, Field& v, const Field& w)
392  {
393  if (mu == 0) {
394  mult_xm(v, w);
395  } else if (mu == 1) {
396  mult_ym(v, w);
397  } else if (mu == 2) {
398  mult_zm(v, w);
399  } else if (mu == 3) {
400  if (m_repr == "Dirac") {
401  mult_tm_dirac(v, w);
402  } else {
403  mult_tm_chiral(v, w);
404  }
405  } else {
406  vout.crucial(m_vl, "Error at %s::mult_dn: illegal mu=%d\n",
407  class_name.c_str(), mu);
408  exit(EXIT_FAILURE);
409  }
410  }
411 
412 
413 //====================================================================
414  void Fopr_WilsonGeneral::mult_gm5(Field& v, const Field& w)
415  {
416  if (m_repr == "Dirac") {
417  mult_gm5_dirac(v, w);
418  } else if (m_repr == "Chiral") {
419  mult_gm5_chiral(v, w);
420  }
421  }
422 
423 
424 //====================================================================
425  void Fopr_WilsonGeneral::D(Field& v, const Field& w)
426  {
427  if (m_repr == "Dirac") {
428  D_dirac(v, 0, w, 0);
429  } else if (m_repr == "Chiral") {
430  D_chiral(v, 0, w, 0);
431  }
432  }
433 
434 
435 //====================================================================
436  void Fopr_WilsonGeneral::Ddag(Field& v, const Field& w)
437  {
438  mult_gm5(v, w);
439  D(m_w1, v);
440  mult_gm5(v, m_w1);
441  }
442 
443 
444 //====================================================================
445  void Fopr_WilsonGeneral::DdagD(Field& v, const Field& w)
446  {
447  D(m_w1, w);
448  mult_gm5(v, m_w1);
449  D(m_w1, v);
450  mult_gm5(v, m_w1);
451  }
452 
453 
454 //====================================================================
455  void Fopr_WilsonGeneral::DDdag(Field& v, const Field& w)
456  {
457  mult_gm5(m_w1, w);
458  D(v, m_w1);
459  mult_gm5(m_w1, v);
460  D(v, m_w1);
461  }
462 
463 
464 //====================================================================
465  void Fopr_WilsonGeneral::H(Field& v, const Field& w)
466  {
467  D(m_w1, w);
468  mult_gm5(v, m_w1);
469  }
470 
471 
472 //====================================================================
473  void Fopr_WilsonGeneral::D_dirac(Field& w, const int ex1,
474  const Field& f, const int ex2)
475  {
476  copy(w, f);
477 
478  clear(m_w2);
479  mult_xp(m_w2, f);
480  mult_xm(m_w2, f);
481  mult_yp(m_w2, f);
482  mult_ym(m_w2, f);
483  mult_zp(m_w2, f);
484  mult_zm(m_w2, f);
485  daxpy(w, -m_kappa_s, m_w2);
486 
487  clear(m_w2);
488  mult_tp_dirac(m_w2, f);
489  mult_tm_dirac(m_w2, f);
490  daxpy(w, -m_kappa_t, m_w2);
491  }
492 
493 
494 //====================================================================
495  void Fopr_WilsonGeneral::D_chiral(Field& w, const int ex1,
496  const Field& f, const int ex2)
497  {
498  copy(w, f);
499 
500  clear(m_w2);
501  mult_xp(m_w2, f);
502  mult_xm(m_w2, f);
503  mult_yp(m_w2, f);
504  mult_ym(m_w2, f);
505  mult_zp(m_w2, f);
506  mult_zm(m_w2, f);
507  daxpy(w, -m_kappa_s, m_w2);
508 
509  clear(m_w2);
510  mult_tp_chiral(m_w2, f);
511  mult_tm_chiral(m_w2, f);
512  daxpy(w, -m_kappa_t, m_w2);
513  }
514 
515 
516 //====================================================================
517  void Fopr_WilsonGeneral::clear(Field& w)
518  {
519 #pragma omp barrier
520 
521  int Nvcd = m_Nvc * m_Nd;
522 
523  double *wp = w.ptr(0);
524 
525  int ith, nth, is, ns;
526  set_threadtask(ith, nth, is, ns, m_Nvol);
527 
528  for (int site = is; site < ns; ++site) {
529  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
530  wp[ivcd + Nvcd * site] = 0.0;
531  }
532  }
533 
534 #pragma omp barrier
535  }
536 
537 
538 //====================================================================
539  void Fopr_WilsonGeneral::daxpy(Field& v,
540  const double fac, const Field& w)
541  {
542 #pragma omp barrier
543 
544  int Nvcd = m_Nvc * m_Nd;
545 
546  double *vp = v.ptr(0);
547  const double *wp = w.ptr(0);
548 
549  int ith, nth, is, ns;
550  set_threadtask(ith, nth, is, ns, m_Nvol);
551 
552  for (int site = is; site < ns; ++site) {
553  for (int ivcd = 0; ivcd < Nvcd; ++ivcd) {
554  vp[ivcd + Nvcd * site] += fac * wp[ivcd + Nvcd * site];
555  }
556  }
557 
558 #pragma omp barrier
559  }
560 
561 
562 //====================================================================
563  void Fopr_WilsonGeneral::mult_gm5_dirac(Field& v, const Field& w)
564  {
565 #pragma omp barrier
566 
567  int Nvcd = m_Nvc * m_Nd;
568 
569  double *vp = v.ptr(0);
570  const double *wp = w.ptr(0);
571 
572  int ith, nth, is, ns;
573  set_threadtask(ith, nth, is, ns, m_Nvol);
574 
575  for (int site = is; site < ns; ++site) {
576  mult_gamma5_dirac(&vp[Nvcd * site], &wp[Nvcd * site], m_Nc);
577  }
578 
579 #pragma omp barrier
580  }
581 
582 
583 //====================================================================
584  void Fopr_WilsonGeneral::mult_gm5_chiral(Field& v, const Field& w)
585  {
586 #pragma omp barrier
587 
588  int Nvcd = m_Nvc * m_Nd;
589 
590  double *vp = v.ptr(0);
591  const double *wp = w.ptr(0);
592 
593  int ith, nth, is, ns;
594  set_threadtask(ith, nth, is, ns, m_Nvol);
595 
596  for (int site = is; site < ns; ++site) {
597  mult_gamma5_chiral(&vp[Nvcd * site], &wp[Nvcd * site], m_Nc);
598  }
599 
600 #pragma omp barrier
601  }
602 
603 
604 //====================================================================
605  void Fopr_WilsonGeneral::mult_xp(Field& v, const Field& w)
606  {
607  const int idir = 0;
608 
609  const int Nvcd = m_Nvc * m_Nd;
610 
611  const int id1 = 0;
612  const int id2 = m_Nvc;
613  const int id3 = m_Nvc * 2;
614  const int id4 = m_Nvc * 3;
615 
616  double bc2 = m_boundary_each_node[idir];
617 
618  double *vp = v.ptr(0);
619  const double *wp = w.ptr(0);
620  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
621 
622  int ith, nth, is, ns;
623  set_threadtask(ith, nth, is, ns, m_Nvol);
624 
625 #pragma omp barrier
626 
627  for (int site = is; site < ns; ++site) {
628  int ix = site % m_Nx;
629  int iyzt = site / m_Nx;
630  if (ix == 0) {
631  int in = Nvcd * site;
632  int ix1 = Nvcd * iyzt;
633  int ix2 = ix1 + NVC;
634  int ix3 = ix1 + 2 * NVC;
635  int ix4 = ix1 + 3 * NVC;
636  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
637  for (int ic = 0; ic < m_Nc; ++ic) {
638  int icr = 2 * ic;
639  int ici = 2 * ic + 1;
640  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[ici + id4 + in];
641  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[icr + id4 + in];
642  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[ici + id3 + in];
643  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[icr + id3 + in];
644  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[ici + id2 + in];
645  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[icr + id2 + in];
646  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[ici + id1 + in];
647  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[icr + id1 + in];
648  }
649 
650  for (int ivc = 0; ivc < NVC; ++ivc) {
651  vcp1_xp[ivc + ix1] = bc2 * vt1[ivc];
652  vcp1_xp[ivc + ix2] = bc2 * vt2[ivc];
653  vcp1_xp[ivc + ix3] = bc2 * vt3[ivc];
654  vcp1_xp[ivc + ix4] = bc2 * vt4[ivc];
655  }
656  }
657  }
658 
659 #pragma omp barrier
660 
661 #pragma omp master
662  {
663  const int Nv = m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
664  Communicator::exchange(Nv, vcp2_xp, vcp1_xp, 0, 1, 1);
665  }
666 #pragma omp barrier
667 
668  for (int site = is; site < ns; ++site) {
669  int ix = site % m_Nx;
670  int iyzt = site / m_Nx;
671  int nei = ix + 1 + m_Nx * iyzt;
672  int iv = Nvcd * site;
673  int ig = m_Ndf * site;
674 
675  if (ix < m_Nx - 1) {
676  int in = Nvcd * nei;
677  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
678  for (int ic = 0; ic < m_Nc; ++ic) {
679  int icr = 2 * ic;
680  int ici = 2 * ic + 1;
681  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[ici + id4 + in];
682  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[icr + id4 + in];
683  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[ici + id3 + in];
684  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[icr + id3 + in];
685  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[ici + id2 + in];
686  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[icr + id2 + in];
687  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[ici + id1 + in];
688  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[icr + id1 + in];
689  }
690  for (int ic = 0; ic < m_Nc; ++ic) {
691  int ic2 = ic * NVC;
692  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
693  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
694  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
695  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
696  double wt3r = mult_uv_r(&up[ic2 + ig], vt3, m_Nc);
697  double wt3i = mult_uv_i(&up[ic2 + ig], vt3, m_Nc);
698  double wt4r = mult_uv_r(&up[ic2 + ig], vt4, m_Nc);
699  double wt4i = mult_uv_i(&up[ic2 + ig], vt4, m_Nc);
700 
701  int icr = 2 * ic;
702  int ici = 2 * ic + 1;
703  vp[icr + id1 + iv] += wt1r;
704  vp[ici + id1 + iv] += wt1i;
705  vp[icr + id2 + iv] += wt2r;
706  vp[ici + id2 + iv] += wt2i;
707  vp[icr + id3 + iv] += wt3r;
708  vp[ici + id3 + iv] += wt3i;
709  vp[icr + id4 + iv] += wt4r;
710  vp[ici + id4 + iv] += wt4i;
711  }
712  } else {
713  int ix1 = Nvcd * iyzt;
714  int ix2 = ix1 + NVC;
715  int ix3 = ix1 + 2 * NVC;
716  int ix4 = ix1 + 3 * NVC;
717  for (int ic = 0; ic < m_Nc; ++ic) {
718  int ic2 = ic * NVC;
719  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
720  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix1], m_Nc);
721  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
722  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix2], m_Nc);
723  double wt3r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix3], m_Nc);
724  double wt3i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix3], m_Nc);
725  double wt4r = mult_uv_r(&up[ic2 + ig], &vcp2_xp[ix4], m_Nc);
726  double wt4i = mult_uv_i(&up[ic2 + ig], &vcp2_xp[ix4], m_Nc);
727 
728  int icr = 2 * ic;
729  int ici = 2 * ic + 1;
730  vp[icr + id1 + iv] += wt1r;
731  vp[ici + id1 + iv] += wt1i;
732  vp[icr + id2 + iv] += wt2r;
733  vp[ici + id2 + iv] += wt2i;
734  vp[icr + id3 + iv] += wt3r;
735  vp[ici + id3 + iv] += wt3i;
736  vp[icr + id4 + iv] += wt4r;
737  vp[ici + id4 + iv] += wt4i;
738  }
739  }
740  }
741 
742 #pragma omp barrier
743  }
744 
745 
746 //====================================================================
747  void Fopr_WilsonGeneral::mult_xm(Field& v, const Field& w)
748  {
749  const int idir = 0;
750 
751  const int Nvcd = m_Nvc * m_Nd;
752 
753  const int id1 = 0;
754  const int id2 = m_Nvc;
755  const int id3 = m_Nvc * 2;
756  const int id4 = m_Nvc * 3;
757 
758  double bc2 = m_boundary_each_node[idir];
759 
760  double *vp = v.ptr(0);
761  const double *wp = w.ptr(0);
762  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
763 
764  int ith, nth, is, ns;
765  set_threadtask(ith, nth, is, ns, m_Nvol);
766 
767 #pragma omp barrier
768 
769  for (int site = is; site < ns; ++site) {
770  int ix = site % m_Nx;
771  int iyzt = site / m_Nx;
772  if (ix == m_Nx - 1) {
773  int in = Nvcd * site;
774  int ig = m_Ndf * site;
775  int ix1 = Nvcd * iyzt;
776  int ix2 = ix1 + NVC;
777  int ix3 = ix1 + 2 * NVC;
778  int ix4 = ix1 + 3 * NVC;
779 
780  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
781  for (int ic = 0; ic < m_Nc; ++ic) {
782  int icr = 2 * ic;
783  int ici = 2 * ic + 1;
784  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[ici + id4 + in];
785  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[icr + id4 + in];
786  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[ici + id3 + in];
787  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[icr + id3 + in];
788  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[ici + id2 + in];
789  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[icr + id2 + in];
790  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[ici + id1 + in];
791  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[icr + id1 + in];
792  }
793 
794  for (int ic = 0; ic < m_Nc; ++ic) {
795  int ic2 = 2 * ic;
796  int icr = 2 * ic;
797  int ici = 2 * ic + 1;
798  vcp1_xm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
799  vcp1_xm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
800  vcp1_xm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
801  vcp1_xm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
802  vcp1_xm[icr + ix3] = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
803  vcp1_xm[ici + ix3] = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
804  vcp1_xm[icr + ix4] = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
805  vcp1_xm[ici + ix4] = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
806  }
807  }
808  }
809 
810 #pragma omp barrier
811 
812 #pragma omp master
813  {
814  const int Nv = m_Nvc * m_Nd * m_Ny * m_Nz * m_Nt;
815  Communicator::exchange(Nv, vcp2_xm, vcp1_xm, 0, -1, 2);
816  }
817 #pragma omp barrier
818 
819  for (int site = is; site < ns; ++site) {
820  int ix = site % m_Nx;
821  int iyzt = site / m_Nx;
822  int nei = ix - 1 + m_Nx * iyzt;
823  int iv = Nvcd * site;
824 
825  if (ix > 0) {
826  int ig = m_Ndf * nei;
827  int in = Nvcd * nei;
828 
829  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
830  for (int ic = 0; ic < m_Nc; ++ic) {
831  int icr = 2 * ic;
832  int ici = 2 * ic + 1;
833  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[ici + id4 + in];
834  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[icr + id4 + in];
835  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[ici + id3 + in];
836  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[icr + id3 + in];
837  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[ici + id2 + in];
838  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[icr + id2 + in];
839  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[ici + id1 + in];
840  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[icr + id1 + in];
841  }
842 
843  for (int ic = 0; ic < m_Nc; ++ic) {
844  int ic2 = 2 * ic;
845  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
846  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
847  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
848  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
849  double wt3r = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
850  double wt3i = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
851  double wt4r = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
852  double wt4i = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
853 
854  int icr = 2 * ic;
855  int ici = 2 * ic + 1;
856  vp[icr + id1 + iv] += wt1r;
857  vp[ici + id1 + iv] += wt1i;
858  vp[icr + id2 + iv] += wt2r;
859  vp[ici + id2 + iv] += wt2i;
860  vp[icr + id3 + iv] += wt3r;
861  vp[ici + id3 + iv] += wt3i;
862  vp[icr + id4 + iv] += wt4r;
863  vp[ici + id4 + iv] += wt4i;
864  }
865  } else {
866  int ix1 = Nvcd * iyzt;
867  int ix2 = ix1 + NVC;
868  int ix3 = ix1 + 2 * NVC;
869  int ix4 = ix1 + 3 * NVC;
870  for (int ic = 0; ic < m_Nc; ++ic) {
871  double wt1r = bc2 * vcp2_xm[2 * ic + ix1];
872  double wt1i = bc2 * vcp2_xm[2 * ic + 1 + ix1];
873  double wt2r = bc2 * vcp2_xm[2 * ic + ix2];
874  double wt2i = bc2 * vcp2_xm[2 * ic + 1 + ix2];
875  double wt3r = bc2 * vcp2_xm[2 * ic + ix3];
876  double wt3i = bc2 * vcp2_xm[2 * ic + 1 + ix3];
877  double wt4r = bc2 * vcp2_xm[2 * ic + ix4];
878  double wt4i = bc2 * vcp2_xm[2 * ic + 1 + ix4];
879 
880  int icr = 2 * ic;
881  int ici = 2 * ic + 1;
882  vp[icr + id1 + iv] += wt1r;
883  vp[ici + id1 + iv] += wt1i;
884  vp[icr + id2 + iv] += wt2r;
885  vp[ici + id2 + iv] += wt2i;
886  vp[icr + id3 + iv] += wt3r;
887  vp[ici + id3 + iv] += wt3i;
888  vp[icr + id4 + iv] += wt4r;
889  vp[ici + id4 + iv] += wt4i;
890  }
891  }
892  }
893 
894 #pragma omp barrier
895  }
896 
897 
898 //====================================================================
899  void Fopr_WilsonGeneral::mult_yp(Field& v, const Field& w)
900  {
901  const int idir = 1;
902 
903  const int Nvcd = m_Nvc * m_Nd;
904 
905  const int id1 = 0;
906  const int id2 = m_Nvc;
907  const int id3 = m_Nvc * 2;
908  const int id4 = m_Nvc * 3;
909 
910  double bc2 = m_boundary_each_node[idir];
911 
912  double *vp = v.ptr(0);
913  const double *wp = w.ptr(0);
914  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
915 
916  int ith, nth, is, ns;
917  set_threadtask(ith, nth, is, ns, m_Nvol);
918 
919 #pragma omp barrier
920 
921  for (int site = is; site < ns; ++site) {
922  int ix = site % m_Nx;
923  int iyzt = site / m_Nx;
924  int iy = iyzt % m_Ny;
925  int izt = iyzt / m_Ny;
926  int ixzt = ix + m_Nx * izt;
927  if (iy == 0) {
928  int in = Nvcd * site;
929  int ix1 = Nvcd * ixzt;
930  int ix2 = ix1 + NVC;
931  int ix3 = ix1 + 2 * NVC;
932  int ix4 = ix1 + 3 * NVC;
933  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
934  for (int ic = 0; ic < m_Nc; ++ic) {
935  int icr = 2 * ic;
936  int ici = 2 * ic + 1;
937  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[icr + id4 + in];
938  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[ici + id4 + in];
939  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[icr + id3 + in];
940  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[ici + id3 + in];
941  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[icr + id2 + in];
942  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[ici + id2 + in];
943  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[icr + id1 + in];
944  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[ici + id1 + in];
945  }
946 
947  for (int ivc = 0; ivc < NVC; ++ivc) {
948  vcp1_yp[ivc + ix1] = bc2 * vt1[ivc];
949  vcp1_yp[ivc + ix2] = bc2 * vt2[ivc];
950  vcp1_yp[ivc + ix3] = bc2 * vt3[ivc];
951  vcp1_yp[ivc + ix4] = bc2 * vt4[ivc];
952  }
953  }
954  }
955 
956 #pragma omp barrier
957 
958 #pragma omp master
959  {
960  const int Nv = m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
961  Communicator::exchange(Nv, vcp2_yp, vcp1_yp, 1, 1, 3);
962  }
963 #pragma omp barrier
964 
965  for (int site = is; site < ns; ++site) {
966  int ix = site % m_Nx;
967  int iyzt = site / m_Nx;
968  int iy = iyzt % m_Ny;
969  int izt = iyzt / m_Ny;
970  int ixzt = ix + m_Nx * izt;
971  int nei = ix + m_Nx * (iy + 1 + m_Ny * izt);
972  int iv = Nvcd * site;
973  int ig = m_Ndf * site;
974 
975  if (iy < m_Ny - 1) {
976  int in = Nvcd * nei;
977 
978  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
979  for (int ic = 0; ic < m_Nc; ++ic) {
980  int icr = 2 * ic;
981  int ici = 2 * ic + 1;
982  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[icr + id4 + in];
983  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[ici + id4 + in];
984  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[icr + id3 + in];
985  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[ici + id3 + in];
986  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[icr + id2 + in];
987  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[ici + id2 + in];
988  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[icr + id1 + in];
989  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[ici + id1 + in];
990  }
991  for (int ic = 0; ic < m_Nc; ++ic) {
992  int ic2 = ic * NVC;
993  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
994  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
995  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
996  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
997  double wt3r = mult_uv_r(&up[ic2 + ig], vt3, m_Nc);
998  double wt3i = mult_uv_i(&up[ic2 + ig], vt3, m_Nc);
999  double wt4r = mult_uv_r(&up[ic2 + ig], vt4, m_Nc);
1000  double wt4i = mult_uv_i(&up[ic2 + ig], vt4, m_Nc);
1001 
1002  int icr = 2 * ic;
1003  int ici = 2 * ic + 1;
1004  vp[icr + id1 + iv] += wt1r;
1005  vp[ici + id1 + iv] += wt1i;
1006  vp[icr + id2 + iv] += wt2r;
1007  vp[ici + id2 + iv] += wt2i;
1008  vp[icr + id3 + iv] += wt3r;
1009  vp[ici + id3 + iv] += wt3i;
1010  vp[icr + id4 + iv] += wt4r;
1011  vp[ici + id4 + iv] += wt4i;
1012  }
1013  } else {
1014  int ix1 = Nvcd * ixzt;
1015  int ix2 = ix1 + NVC;
1016  int ix3 = ix1 + 2 * NVC;
1017  int ix4 = ix1 + 3 * NVC;
1018  for (int ic = 0; ic < m_Nc; ++ic) {
1019  int ic2 = ic * NVC;
1020  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1021  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix1], m_Nc);
1022  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1023  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix2], m_Nc);
1024  double wt3r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix3], m_Nc);
1025  double wt3i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix3], m_Nc);
1026  double wt4r = mult_uv_r(&up[ic2 + ig], &vcp2_yp[ix4], m_Nc);
1027  double wt4i = mult_uv_i(&up[ic2 + ig], &vcp2_yp[ix4], m_Nc);
1028 
1029  int icr = 2 * ic;
1030  int ici = 2 * ic + 1;
1031  vp[icr + id1 + iv] += wt1r;
1032  vp[ici + id1 + iv] += wt1i;
1033  vp[icr + id2 + iv] += wt2r;
1034  vp[ici + id2 + iv] += wt2i;
1035  vp[icr + id3 + iv] += wt3r;
1036  vp[ici + id3 + iv] += wt3i;
1037  vp[icr + id4 + iv] += wt4r;
1038  vp[ici + id4 + iv] += wt4i;
1039  }
1040  }
1041  }
1042 
1043 #pragma omp barrier
1044  }
1045 
1046 
1047 //====================================================================
1048  void Fopr_WilsonGeneral::mult_ym(Field& v, const Field& w)
1049  {
1050  const int idir = 1;
1051 
1052  const int Nvcd = m_Nvc * m_Nd;
1053 
1054  const int id1 = 0;
1055  const int id2 = m_Nvc;
1056  const int id3 = m_Nvc * 2;
1057  const int id4 = m_Nvc * 3;
1058 
1059  double bc2 = m_boundary_each_node[idir];
1060 
1061  double *vp = v.ptr(0);
1062  const double *wp = w.ptr(0);
1063  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1064 
1065  int ith, nth, is, ns;
1066  set_threadtask(ith, nth, is, ns, m_Nvol);
1067 
1068 #pragma omp barrier
1069 
1070  for (int site = is; site < ns; ++site) {
1071  int ix = site % m_Nx;
1072  int iyzt = site / m_Nx;
1073  int iy = iyzt % m_Ny;
1074  int izt = iyzt / m_Ny;
1075  int ixzt = ix + m_Nx * izt;
1076  if (iy == m_Ny - 1) {
1077  int in = Nvcd * site;
1078  int ig = m_Ndf * site;
1079  int ix1 = Nvcd * ixzt;
1080  int ix2 = ix1 + NVC;
1081  int ix3 = ix1 + 2 * NVC;
1082  int ix4 = ix1 + 3 * NVC;
1083 
1084  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1085  for (int ic = 0; ic < m_Nc; ++ic) {
1086  int icr = 2 * ic;
1087  int ici = 2 * ic + 1;
1088  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[icr + id4 + in];
1089  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[ici + id4 + in];
1090  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[icr + id3 + in];
1091  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[ici + id3 + in];
1092  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[icr + id2 + in];
1093  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[ici + id2 + in];
1094  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[icr + id1 + in];
1095  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[ici + id1 + in];
1096  }
1097 
1098  for (int ic = 0; ic < m_Nc; ++ic) {
1099  int ic2 = 2 * ic;
1100  int icr = 2 * ic;
1101  int ici = 2 * ic + 1;
1102  vcp1_ym[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1103  vcp1_ym[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1104  vcp1_ym[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1105  vcp1_ym[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1106  vcp1_ym[icr + ix3] = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
1107  vcp1_ym[ici + ix3] = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
1108  vcp1_ym[icr + ix4] = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
1109  vcp1_ym[ici + ix4] = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
1110  }
1111  }
1112  }
1113 
1114 #pragma omp barrier
1115 
1116 #pragma omp master
1117  {
1118  const int Nv = m_Nvc * m_Nd * m_Nx * m_Nz * m_Nt;
1119  Communicator::exchange(Nv, vcp2_ym, vcp1_ym, 1, -1, 4);
1120  }
1121 #pragma omp barrier
1122 
1123  for (int site = is; site < ns; ++site) {
1124  int ix = site % m_Nx;
1125  int iyzt = site / m_Nx;
1126  int iy = iyzt % m_Ny;
1127  int izt = iyzt / m_Ny;
1128  int ixzt = ix + m_Nx * izt;
1129  int nei = ix + m_Nx * (iy - 1 + m_Ny * izt);
1130  int iv = Nvcd * site;
1131 
1132  if (iy > 0) {
1133  int ig = m_Ndf * nei;
1134  int in = Nvcd * nei;
1135  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1136  for (int ic = 0; ic < m_Nc; ++ic) {
1137  int icr = 2 * ic;
1138  int ici = 2 * ic + 1;
1139  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[icr + id4 + in];
1140  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[ici + id4 + in];
1141  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[icr + id3 + in];
1142  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[ici + id3 + in];
1143  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[icr + id2 + in];
1144  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[ici + id2 + in];
1145  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[icr + id1 + in];
1146  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[ici + id1 + in];
1147  }
1148 
1149  for (int ic = 0; ic < m_Nc; ++ic) {
1150  int ic2 = 2 * ic;
1151  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1152  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1153  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1154  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1155  double wt3r = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
1156  double wt3i = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
1157  double wt4r = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
1158  double wt4i = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
1159 
1160  int icr = 2 * ic;
1161  int ici = 2 * ic + 1;
1162  vp[icr + id1 + iv] += wt1r;
1163  vp[ici + id1 + iv] += wt1i;
1164  vp[icr + id2 + iv] += wt2r;
1165  vp[ici + id2 + iv] += wt2i;
1166  vp[icr + id3 + iv] += wt3r;
1167  vp[ici + id3 + iv] += wt3i;
1168  vp[icr + id4 + iv] += wt4r;
1169  vp[ici + id4 + iv] += wt4i;
1170  }
1171  } else {
1172  int ix1 = Nvcd * ixzt;
1173  int ix2 = ix1 + NVC;
1174  int ix3 = ix1 + 2 * NVC;
1175  int ix4 = ix1 + 3 * NVC;
1176  for (int ic = 0; ic < m_Nc; ++ic) {
1177  double wt1r = bc2 * vcp2_ym[2 * ic + ix1];
1178  double wt1i = bc2 * vcp2_ym[2 * ic + 1 + ix1];
1179  double wt2r = bc2 * vcp2_ym[2 * ic + ix2];
1180  double wt2i = bc2 * vcp2_ym[2 * ic + 1 + ix2];
1181  double wt3r = bc2 * vcp2_ym[2 * ic + ix3];
1182  double wt3i = bc2 * vcp2_ym[2 * ic + 1 + ix3];
1183  double wt4r = bc2 * vcp2_ym[2 * ic + ix4];
1184  double wt4i = bc2 * vcp2_ym[2 * ic + 1 + ix4];
1185 
1186  int icr = 2 * ic;
1187  int ici = 2 * ic + 1;
1188  vp[icr + id1 + iv] += wt1r;
1189  vp[ici + id1 + iv] += wt1i;
1190  vp[icr + id2 + iv] += wt2r;
1191  vp[ici + id2 + iv] += wt2i;
1192  vp[icr + id3 + iv] += wt3r;
1193  vp[ici + id3 + iv] += wt3i;
1194  vp[icr + id4 + iv] += wt4r;
1195  vp[ici + id4 + iv] += wt4i;
1196  }
1197  }
1198  }
1199 
1200 #pragma omp barrier
1201  }
1202 
1203 
1204 //====================================================================
1205  void Fopr_WilsonGeneral::mult_zp(Field& v, const Field& w)
1206  {
1207  const int idir = 2;
1208 
1209  const int Nvcd = m_Nvc * m_Nd;
1210 
1211  const int id1 = 0;
1212  const int id2 = m_Nvc;
1213  const int id3 = m_Nvc * 2;
1214  const int id4 = m_Nvc * 3;
1215 
1216  double bc2 = m_boundary_each_node[idir];
1217 
1218  double *vp = v.ptr(0);
1219  const double *wp = w.ptr(0);
1220  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1221 
1222  int ith, nth, is, ns;
1223  set_threadtask(ith, nth, is, ns, m_Nvol);
1224 
1225  int Nxy = m_Nx * m_Ny;
1226 
1227 #pragma omp barrier
1228 
1229  for (int site = is; site < ns; ++site) {
1230  int ixy = site % Nxy;
1231  int izt = site / Nxy;
1232  int iz = izt % m_Nz;
1233  int it = izt / m_Nz;
1234  int ixyt = ixy + Nxy * it;
1235  if (iz == 0) {
1236  int in = Nvcd * site;
1237  int ix1 = Nvcd * ixyt;
1238  int ix2 = ix1 + NVC;
1239  int ix3 = ix1 + 2 * NVC;
1240  int ix4 = ix1 + 3 * NVC;
1241  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1242  for (int ic = 0; ic < m_Nc; ++ic) {
1243  int icr = 2 * ic;
1244  int ici = 2 * ic + 1;
1245  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[ici + id3 + in];
1246  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[icr + id3 + in];
1247  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[ici + id4 + in];
1248  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[icr + id4 + in];
1249  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[ici + id1 + in];
1250  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[icr + id1 + in];
1251  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[ici + id2 + in];
1252  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[icr + id2 + in];
1253  }
1254 
1255  for (int ivc = 0; ivc < NVC; ++ivc) {
1256  vcp1_zp[ivc + ix1] = bc2 * vt1[ivc];
1257  vcp1_zp[ivc + ix2] = bc2 * vt2[ivc];
1258  vcp1_zp[ivc + ix3] = bc2 * vt3[ivc];
1259  vcp1_zp[ivc + ix4] = bc2 * vt4[ivc];
1260  }
1261  }
1262  }
1263 
1264 #pragma omp barrier
1265 
1266 #pragma omp master
1267  {
1268  int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
1269  Communicator::exchange(Nv, vcp2_zp, vcp1_zp, 2, 1, 5);
1270  }
1271 #pragma omp barrier
1272 
1273  for (int site = is; site < ns; ++site) {
1274  int ixy = site % Nxy;
1275  int izt = site / Nxy;
1276  int iz = izt % m_Nz;
1277  int it = izt / m_Nz;
1278  int ixyt = ixy + Nxy * it;
1279  int nei = ixy + Nxy * (iz + 1 + m_Nz * it);
1280  int iv = Nvcd * site;
1281  int ig = m_Ndf * site;
1282 
1283  if (iz < m_Nz - 1) {
1284  int in = Nvcd * nei;
1285  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1286  for (int ic = 0; ic < m_Nc; ++ic) {
1287  int icr = 2 * ic;
1288  int ici = 2 * ic + 1;
1289  vt1[icr] = m_r_s * wp[icr + id1 + in] - m_nu_s * wp[ici + id3 + in];
1290  vt1[ici] = m_r_s * wp[ici + id1 + in] + m_nu_s * wp[icr + id3 + in];
1291  vt2[icr] = m_r_s * wp[icr + id2 + in] + m_nu_s * wp[ici + id4 + in];
1292  vt2[ici] = m_r_s * wp[ici + id2 + in] - m_nu_s * wp[icr + id4 + in];
1293  vt3[icr] = m_r_s * wp[icr + id3 + in] + m_nu_s * wp[ici + id1 + in];
1294  vt3[ici] = m_r_s * wp[ici + id3 + in] - m_nu_s * wp[icr + id1 + in];
1295  vt4[icr] = m_r_s * wp[icr + id4 + in] - m_nu_s * wp[ici + id2 + in];
1296  vt4[ici] = m_r_s * wp[ici + id4 + in] + m_nu_s * wp[icr + id2 + in];
1297  }
1298 
1299  for (int ic = 0; ic < m_Nc; ++ic) {
1300  int ic2 = ic * NVC;
1301  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1302  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1303  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1304  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1305  double wt3r = mult_uv_r(&up[ic2 + ig], vt3, m_Nc);
1306  double wt3i = mult_uv_i(&up[ic2 + ig], vt3, m_Nc);
1307  double wt4r = mult_uv_r(&up[ic2 + ig], vt4, m_Nc);
1308  double wt4i = mult_uv_i(&up[ic2 + ig], vt4, m_Nc);
1309 
1310  int icr = 2 * ic;
1311  int ici = 2 * ic + 1;
1312  vp[icr + id1 + iv] += wt1r;
1313  vp[ici + id1 + iv] += wt1i;
1314  vp[icr + id2 + iv] += wt2r;
1315  vp[ici + id2 + iv] += wt2i;
1316  vp[icr + id3 + iv] += wt3r;
1317  vp[ici + id3 + iv] += wt3i;
1318  vp[icr + id4 + iv] += wt4r;
1319  vp[ici + id4 + iv] += wt4i;
1320  }
1321  } else {
1322  int ix1 = Nvcd * ixyt;
1323  int ix2 = ix1 + NVC;
1324  int ix3 = ix1 + 2 * NVC;
1325  int ix4 = ix1 + 3 * NVC;
1326  for (int ic = 0; ic < m_Nc; ++ic) {
1327  int ic2 = ic * NVC;
1328  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1329  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix1], m_Nc);
1330  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1331  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix2], m_Nc);
1332  double wt3r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix3], m_Nc);
1333  double wt3i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix3], m_Nc);
1334  double wt4r = mult_uv_r(&up[ic2 + ig], &vcp2_zp[ix4], m_Nc);
1335  double wt4i = mult_uv_i(&up[ic2 + ig], &vcp2_zp[ix4], m_Nc);
1336 
1337  int icr = 2 * ic;
1338  int ici = 2 * ic + 1;
1339  vp[icr + id1 + iv] += wt1r;
1340  vp[ici + id1 + iv] += wt1i;
1341  vp[icr + id2 + iv] += wt2r;
1342  vp[ici + id2 + iv] += wt2i;
1343  vp[icr + id3 + iv] += wt3r;
1344  vp[ici + id3 + iv] += wt3i;
1345  vp[icr + id4 + iv] += wt4r;
1346  vp[ici + id4 + iv] += wt4i;
1347  }
1348  }
1349  }
1350 
1351 #pragma omp barrier
1352  }
1353 
1354 
1355 //====================================================================
1356  void Fopr_WilsonGeneral::mult_zm(Field& v, const Field& w)
1357  {
1358  const int idir = 2;
1359 
1360  const int Nvcd = m_Nvc * m_Nd;
1361 
1362  const int id1 = 0;
1363  const int id2 = m_Nvc;
1364  const int id3 = m_Nvc * 2;
1365  const int id4 = m_Nvc * 3;
1366 
1367  double bc2 = m_boundary_each_node[idir];
1368 
1369  double *vp = v.ptr(0);
1370  const double *wp = w.ptr(0);
1371  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1372 
1373  int ith, nth, is, ns;
1374  set_threadtask(ith, nth, is, ns, m_Nvol);
1375 
1376  int Nxy = m_Nx * m_Ny;
1377 
1378 #pragma omp barrier
1379 
1380  for (int site = is; site < ns; ++site) {
1381  int ixy = site % Nxy;
1382  int izt = site / Nxy;
1383  int iz = izt % m_Nz;
1384  int it = izt / m_Nz;
1385  int ixyt = ixy + Nxy * it;
1386  if (iz == m_Nz - 1) {
1387  int in = Nvcd * site;
1388  int ig = m_Ndf * site;
1389  int ix1 = Nvcd * ixyt;
1390  int ix2 = ix1 + NVC;
1391  int ix3 = ix1 + 2 * NVC;
1392  int ix4 = ix1 + 3 * NVC;
1393  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1394  for (int ic = 0; ic < m_Nc; ++ic) {
1395  int icr = 2 * ic;
1396  int ici = 2 * ic + 1;
1397  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[ici + id3 + in];
1398  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[icr + id3 + in];
1399  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[ici + id4 + in];
1400  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[icr + id4 + in];
1401  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[ici + id1 + in];
1402  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[icr + id1 + in];
1403  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[ici + id2 + in];
1404  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[icr + id2 + in];
1405  }
1406 
1407  for (int ic = 0; ic < m_Nc; ++ic) {
1408  int ic2 = 2 * ic;
1409  int icr = 2 * ic;
1410  int ici = 2 * ic + 1;
1411  vcp1_zm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1412  vcp1_zm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1413  vcp1_zm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1414  vcp1_zm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1415  vcp1_zm[icr + ix3] = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
1416  vcp1_zm[ici + ix3] = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
1417  vcp1_zm[icr + ix4] = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
1418  vcp1_zm[ici + ix4] = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
1419  }
1420  }
1421  }
1422 
1423 #pragma omp barrier
1424 
1425 #pragma omp master
1426  {
1427  const int Nv = m_Nvc * m_Nd * m_Nx * m_Ny * m_Nt;
1428  Communicator::exchange(Nv, vcp2_zm, vcp1_zm, 2, -1, 6);
1429  }
1430 #pragma omp barrier
1431 
1432  for (int site = is; site < ns; ++site) {
1433  int ixy = site % Nxy;
1434  int izt = site / Nxy;
1435  int iz = izt % m_Nz;
1436  int it = izt / m_Nz;
1437  int ixyt = ixy + Nxy * it;
1438  int nei = ixy + Nxy * (iz - 1 + m_Nz * it);
1439  int iv = Nvcd * site;
1440 
1441  if (iz > 0) {
1442  int ig = m_Ndf * nei;
1443  int in = Nvcd * nei;
1444  double vt1[NVC], vt2[NVC], vt3[NVC], vt4[NVC];
1445  for (int ic = 0; ic < m_Nc; ++ic) {
1446  int icr = 2 * ic;
1447  int ici = 2 * ic + 1;
1448  vt1[icr] = m_r_s * wp[icr + id1 + in] + m_nu_s * wp[ici + id3 + in];
1449  vt1[ici] = m_r_s * wp[ici + id1 + in] - m_nu_s * wp[icr + id3 + in];
1450  vt2[icr] = m_r_s * wp[icr + id2 + in] - m_nu_s * wp[ici + id4 + in];
1451  vt2[ici] = m_r_s * wp[ici + id2 + in] + m_nu_s * wp[icr + id4 + in];
1452  vt3[icr] = m_r_s * wp[icr + id3 + in] - m_nu_s * wp[ici + id1 + in];
1453  vt3[ici] = m_r_s * wp[ici + id3 + in] + m_nu_s * wp[icr + id1 + in];
1454  vt4[icr] = m_r_s * wp[icr + id4 + in] + m_nu_s * wp[ici + id2 + in];
1455  vt4[ici] = m_r_s * wp[ici + id4 + in] - m_nu_s * wp[icr + id2 + in];
1456  }
1457  for (int ic = 0; ic < m_Nc; ++ic) {
1458  int ic2 = 2 * ic;
1459  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1460  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1461  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1462  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1463  double wt3r = mult_udagv_r(&up[ic2 + ig], vt3, m_Nc);
1464  double wt3i = mult_udagv_i(&up[ic2 + ig], vt3, m_Nc);
1465  double wt4r = mult_udagv_r(&up[ic2 + ig], vt4, m_Nc);
1466  double wt4i = mult_udagv_i(&up[ic2 + ig], vt4, m_Nc);
1467 
1468  int icr = 2 * ic;
1469  int ici = 2 * ic + 1;
1470  vp[icr + id1 + iv] += wt1r;
1471  vp[ici + id1 + iv] += wt1i;
1472  vp[icr + id2 + iv] += wt2r;
1473  vp[ici + id2 + iv] += wt2i;
1474  vp[icr + id3 + iv] += wt3r;
1475  vp[ici + id3 + iv] += wt3i;
1476  vp[icr + id4 + iv] += wt4r;
1477  vp[ici + id4 + iv] += wt4i;
1478  }
1479  } else {
1480  int ix1 = Nvcd * ixyt;
1481  int ix2 = ix1 + NVC;
1482  int ix3 = ix1 + 2 * NVC;
1483  int ix4 = ix1 + 3 * NVC;
1484  for (int ic = 0; ic < m_Nc; ++ic) {
1485  double wt1r = bc2 * vcp2_zm[2 * ic + ix1];
1486  double wt1i = bc2 * vcp2_zm[2 * ic + 1 + ix1];
1487  double wt2r = bc2 * vcp2_zm[2 * ic + ix2];
1488  double wt2i = bc2 * vcp2_zm[2 * ic + 1 + ix2];
1489  double wt3r = bc2 * vcp2_zm[2 * ic + ix3];
1490  double wt3i = bc2 * vcp2_zm[2 * ic + 1 + ix3];
1491  double wt4r = bc2 * vcp2_zm[2 * ic + ix4];
1492  double wt4i = bc2 * vcp2_zm[2 * ic + 1 + ix4];
1493 
1494  int icr = 2 * ic;
1495  int ici = 2 * ic + 1;
1496  vp[icr + id1 + iv] += wt1r;
1497  vp[ici + id1 + iv] += wt1i;
1498  vp[icr + id2 + iv] += wt2r;
1499  vp[ici + id2 + iv] += wt2i;
1500  vp[icr + id3 + iv] += wt3r;
1501  vp[ici + id3 + iv] += wt3i;
1502  vp[icr + id4 + iv] += wt4r;
1503  vp[ici + id4 + iv] += wt4i;
1504  }
1505  }
1506  }
1507 
1508 #pragma omp barrier
1509  }
1510 
1511 
1512 //====================================================================
1513  void Fopr_WilsonGeneral::mult_tp_dirac(Field& v, const Field& w)
1514  {
1515  const int idir = 3;
1516 
1517  const int Nvc2 = m_Nvc * 2;
1518  const int Nvcd = m_Nvc * m_Nd;
1519 
1520  const int id1 = 0;
1521  const int id2 = m_Nvc;
1522  const int id3 = m_Nvc * 2;
1523  const int id4 = m_Nvc * 3;
1524 
1525  double bc2 = m_boundary_each_node[idir];
1526 
1527  double *vp = v.ptr(0);
1528  const double *wp = w.ptr(0);
1529  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1530 
1531  int ith, nth, is, ns;
1532  set_threadtask(ith, nth, is, ns, m_Nvol);
1533 
1534  int Nxyz = m_Nx * m_Ny * m_Nz;
1535 
1536 #pragma omp barrier
1537 
1538  for (int site = is; site < ns; ++site) {
1539  int ixyz = site % Nxyz;
1540  int it = site / Nxyz;
1541  if (it == 0) {
1542  int in = Nvcd * site;
1543  int ix1 = Nvc2 * ixyz;
1544  int ix2 = ix1 + NVC;
1545  double vt1[NVC], vt2[NVC];
1546  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
1547  for (int ivc = 0; ivc < NVC; ++ivc) {
1548  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
1549  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
1550  }
1551  }
1552  }
1553 
1554 #pragma omp barrier
1555 
1556 #pragma omp master
1557  {
1558  int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
1559  Communicator::exchange(Nv, vcp2_tp, vcp1_tp, 3, 1, 7);
1560  }
1561 #pragma omp barrier
1562 
1563  for (int site = is; site < ns; ++site) {
1564  int ixyz = site % Nxyz;
1565  int it = site / Nxyz;
1566  int nei = ixyz + Nxyz * (it + 1);
1567  int iv = Nvcd * site;
1568  int ig = m_Ndf * site;
1569 
1570  if (it < m_Nt - 1) {
1571  int in = Nvcd * nei;
1572  double vt1[NVC], vt2[NVC];
1573  set_sp2_tp_dirac(vt1, vt2, &wp[in], m_Nc);
1574  for (int ic = 0; ic < m_Nc; ++ic) {
1575  int ic2 = ic * NVC;
1576  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1577  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1578  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1579  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1580  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1581  }
1582  } else {
1583  int ix1 = Nvc2 * ixyz;
1584  int ix2 = ix1 + NVC;
1585  for (int ic = 0; ic < m_Nc; ++ic) {
1586  int ic2 = ic * NVC;
1587  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1588  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1589  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1590  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1591  set_sp4_tp_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1592  }
1593  }
1594  }
1595 
1596 #pragma omp barrier
1597  }
1598 
1599 
1600 //====================================================================
1601  void Fopr_WilsonGeneral::mult_tm_dirac(Field& v, const Field& w)
1602  {
1603  const int idir = 3;
1604 
1605  const int Nvc2 = m_Nvc * 2;
1606  const int Nvcd = m_Nvc * m_Nd;
1607 
1608  const int id1 = 0;
1609  const int id2 = m_Nvc;
1610  const int id3 = m_Nvc * 2;
1611  const int id4 = m_Nvc * 3;
1612 
1613  double bc2 = m_boundary_each_node[idir];
1614 
1615  double *vp = v.ptr(0);
1616  const double *wp = w.ptr(0);
1617  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1618 
1619  int ith, nth, is, ns;
1620  set_threadtask(ith, nth, is, ns, m_Nvol);
1621 
1622  int Nxyz = m_Nx * m_Ny * m_Nz;
1623 
1624 #pragma omp barrier
1625 
1626  for (int site = is; site < ns; ++site) {
1627  int ixyz = site % Nxyz;
1628  int it = site / Nxyz;
1629  if (it == m_Nt - 1) {
1630  int in = Nvcd * site;
1631  int ig = m_Ndf * site;
1632  int ix1 = Nvc2 * ixyz;
1633  int ix2 = ix1 + NVC;
1634 
1635  double vt1[NVC], vt2[NVC];
1636  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
1637  for (int ic = 0; ic < m_Nc; ++ic) {
1638  int ic2 = 2 * ic;
1639  int icr = 2 * ic;
1640  int ici = 2 * ic + 1;
1641  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1642  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1643  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1644  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1645  }
1646  }
1647  }
1648 
1649 #pragma omp barrier
1650 
1651 #pragma omp master
1652  {
1653  const int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
1654  Communicator::exchange(Nv, vcp2_tm, vcp1_tm, 3, -1, 8);
1655  }
1656 #pragma omp barrier
1657 
1658  for (int site = is; site < ns; ++site) {
1659  int ixyz = site % Nxyz;
1660  int it = site / Nxyz;
1661  int nei = ixyz + Nxyz * (it - 1);
1662  int iv = Nvcd * site;
1663 
1664  if (it > 0) {
1665  int ig = m_Ndf * nei;
1666  int in = Nvcd * nei;
1667  double vt1[NVC], vt2[NVC];
1668  set_sp2_tm_dirac(vt1, vt2, &wp[in], m_Nc);
1669  for (int ic = 0; ic < m_Nc; ++ic) {
1670  int ic2 = 2 * ic;
1671  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1672  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1673  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1674  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1675  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1676  }
1677  } else {
1678  int ix1 = Nvc2 * ixyz;
1679  int ix2 = ix1 + NVC;
1680  for (int ic = 0; ic < m_Nc; ++ic) {
1681  int icr = 2 * ic;
1682  int ici = 2 * ic + 1;
1683  double wt1r = bc2 * vcp2_tm[icr + ix1];
1684  double wt1i = bc2 * vcp2_tm[ici + ix1];
1685  double wt2r = bc2 * vcp2_tm[icr + ix2];
1686  double wt2i = bc2 * vcp2_tm[ici + ix2];
1687  set_sp4_tm_dirac(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1688  }
1689  }
1690  }
1691 
1692 #pragma omp barrier
1693  }
1694 
1695 
1696 //====================================================================
1697  void Fopr_WilsonGeneral::mult_tp_chiral(Field& v, const Field& w)
1698  {
1699  const int idir = 3;
1700 
1701  const int Nvc2 = m_Nvc * 2;
1702  const int Nvcd = m_Nvc * m_Nd;
1703 
1704  const int id1 = 0;
1705  const int id2 = m_Nvc;
1706  const int id3 = m_Nvc * 2;
1707  const int id4 = m_Nvc * 3;
1708 
1709  double bc2 = m_boundary_each_node[idir];
1710 
1711  double *vp = v.ptr(0);
1712  const double *wp = w.ptr(0);
1713  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1714 
1715  int ith, nth, is, ns;
1716  set_threadtask(ith, nth, is, ns, m_Nvol);
1717 
1718  int Nxyz = m_Nx * m_Ny * m_Nz;
1719 
1720 #pragma omp barrier
1721 
1722  for (int site = is; site < ns; ++site) {
1723  int ixyz = site % Nxyz;
1724  int it = site / Nxyz;
1725  if (it == 0) {
1726  int in = Nvcd * site;
1727  int ix1 = Nvc2 * ixyz;
1728  int ix2 = ix1 + NVC;
1729  double vt1[NVC], vt2[NVC];
1730  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
1731  for (int ivc = 0; ivc < NVC; ++ivc) {
1732  vcp1_tp[ivc + ix1] = bc2 * vt1[ivc];
1733  vcp1_tp[ivc + ix2] = bc2 * vt2[ivc];
1734  }
1735  }
1736  }
1737 
1738 #pragma omp barrier
1739 
1740 #pragma omp master
1741  {
1742  int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
1743  Communicator::exchange(Nv, vcp2_tp, vcp1_tp, 3, 1, 7);
1744  }
1745 #pragma omp barrier
1746 
1747  for (int site = is; site < ns; ++site) {
1748  int ixyz = site % Nxyz;
1749  int it = site / Nxyz;
1750  int nei = ixyz + Nxyz * (it + 1);
1751  int iv = Nvcd * site;
1752  int ig = m_Ndf * site;
1753 
1754  if (it < m_Nt - 1) {
1755  int in = Nvcd * nei;
1756  double vt1[NVC], vt2[NVC];
1757  set_sp2_tp_chiral(vt1, vt2, &wp[in], m_Nc);
1758  for (int ic = 0; ic < m_Nc; ++ic) {
1759  int ic2 = ic * NVC;
1760  double wt1r = mult_uv_r(&up[ic2 + ig], vt1, m_Nc);
1761  double wt1i = mult_uv_i(&up[ic2 + ig], vt1, m_Nc);
1762  double wt2r = mult_uv_r(&up[ic2 + ig], vt2, m_Nc);
1763  double wt2i = mult_uv_i(&up[ic2 + ig], vt2, m_Nc);
1764  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1765  }
1766  } else {
1767  int ix1 = Nvc2 * ixyz;
1768  int ix2 = ix1 + NVC;
1769  for (int ic = 0; ic < m_Nc; ++ic) {
1770  int ic2 = ic * NVC;
1771  double wt1r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1772  double wt1i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix1], m_Nc);
1773  double wt2r = mult_uv_r(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1774  double wt2i = mult_uv_i(&up[ic2 + ig], &vcp2_tp[ix2], m_Nc);
1775  set_sp4_tp_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1776  }
1777  }
1778  }
1779 
1780 #pragma omp barrier
1781  }
1782 
1783 
1784 //====================================================================
1785  void Fopr_WilsonGeneral::mult_tm_chiral(Field& v, const Field& w)
1786  {
1787  const int idir = 3;
1788 
1789  const int Nvc2 = m_Nvc * 2;
1790  const int Nvcd = m_Nvc * m_Nd;
1791 
1792  const int id1 = 0;
1793  const int id2 = m_Nvc;
1794  const int id3 = m_Nvc * 2;
1795  const int id4 = m_Nvc * 3;
1796 
1797  double bc2 = m_boundary_each_node[idir];
1798 
1799  double *vp = v.ptr(0);
1800  const double *wp = w.ptr(0);
1801  const double *up = m_U->ptr(m_Ndf * m_Nvol * idir);
1802 
1803  int ith, nth, is, ns;
1804  set_threadtask(ith, nth, is, ns, m_Nvol);
1805 
1806  int Nxyz = m_Nx * m_Ny * m_Nz;
1807 
1808 #pragma omp barrier
1809 
1810  for (int site = is; site < ns; ++site) {
1811  int ixyz = site % Nxyz;
1812  int it = site / Nxyz;
1813  if (it == m_Nt - 1) {
1814  int in = Nvcd * site;
1815  int ig = m_Ndf * site;
1816  int ix1 = Nvc2 * ixyz;
1817  int ix2 = ix1 + NVC;
1818 
1819  double vt1[NVC], vt2[NVC];
1820  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
1821  for (int ic = 0; ic < m_Nc; ++ic) {
1822  int ic2 = 2 * ic;
1823  int icr = 2 * ic;
1824  int ici = 2 * ic + 1;
1825  vcp1_tm[icr + ix1] = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1826  vcp1_tm[ici + ix1] = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1827  vcp1_tm[icr + ix2] = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1828  vcp1_tm[ici + ix2] = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1829  }
1830  }
1831  }
1832 
1833 #pragma omp barrier
1834 
1835 #pragma omp master
1836  {
1837  const int Nv = m_Nvc * 2 * m_Nx * m_Ny * m_Nz;
1838  Communicator::exchange(Nv, vcp2_tm, vcp1_tm, 3, -1, 8);
1839  }
1840 #pragma omp barrier
1841 
1842  for (int site = is; site < ns; ++site) {
1843  int ixyz = site % Nxyz;
1844  int it = site / Nxyz;
1845  int nei = ixyz + Nxyz * (it - 1);
1846  int iv = Nvcd * site;
1847 
1848  if (it > 0) {
1849  int ig = m_Ndf * nei;
1850  int in = Nvcd * nei;
1851  double vt1[NVC], vt2[NVC];
1852  set_sp2_tm_chiral(vt1, vt2, &wp[in], m_Nc);
1853  for (int ic = 0; ic < m_Nc; ++ic) {
1854  int ic2 = 2 * ic;
1855  double wt1r = mult_udagv_r(&up[ic2 + ig], vt1, m_Nc);
1856  double wt1i = mult_udagv_i(&up[ic2 + ig], vt1, m_Nc);
1857  double wt2r = mult_udagv_r(&up[ic2 + ig], vt2, m_Nc);
1858  double wt2i = mult_udagv_i(&up[ic2 + ig], vt2, m_Nc);
1859  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1860  }
1861  } else {
1862  int ix1 = Nvc2 * ixyz;
1863  int ix2 = ix1 + NVC;
1864  for (int ic = 0; ic < m_Nc; ++ic) {
1865  double wt1r = bc2 * vcp2_tm[2 * ic + ix1];
1866  double wt1i = bc2 * vcp2_tm[2 * ic + 1 + ix1];
1867  double wt2r = bc2 * vcp2_tm[2 * ic + ix2];
1868  double wt2i = bc2 * vcp2_tm[2 * ic + 1 + ix2];
1869  set_sp4_tm_chiral(&vp[2 * ic + iv], wt1r, wt1i, wt2r, wt2i, m_Nc);
1870  }
1871  }
1872  }
1873 
1874 #pragma omp barrier
1875  }
1876 
1877 
1878 //====================================================================
1879  double Fopr_WilsonGeneral::flop_count()
1880  {
1881  // Counting of floating point operations in giga unit.
1882  // The following counting explicitly depends on the implementation.
1883  // It will be recalculated when the code is modified.
1884  // The present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]
1885 
1886  const int Nvol = CommonParameters::Nvol();
1887  const int NPE = CommonParameters::NPE();
1888 
1889  int flop_site;
1890 
1891  if (m_repr == "Dirac") {
1892  flop_site = m_Nc * m_Nd * (4 + 6 * (4 * m_Nc + 2) + 2 * (4 * m_Nc + 1));
1893  } else if (m_repr == "Chiral") {
1894  flop_site = m_Nc * m_Nd * (4 + 8 * (4 * m_Nc + 2));
1895  } else {
1896  vout.crucial(m_vl, "Error at %s: input repr is undefined.\n",
1897  class_name.c_str());
1898  exit(EXIT_FAILURE);
1899  }
1900 
1901  double gflop = flop_site * (Nvol * (NPE / 1.0e+9));
1902 
1903  if ((m_mode == "DdagD") || (m_mode == "DDdag")) gflop *= 2;
1904 
1905  return gflop;
1906  }
1907 
1908 
1909 //====================================================================
1910 }
1911 //============================================================END=====
Imp::Fopr_WilsonGeneral::tidyup
void tidyup()
final clean-up.
Definition: fopr_WilsonGeneral_impl.cpp:150
Imp::Fopr_WilsonGeneral::vcp2_ym
double * vcp2_ym
Definition: fopr_WilsonGeneral_impl.h:64
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
fopr_WilsonGeneral_impl.h
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
Imp::Fopr_WilsonGeneral::m_Ndf
int m_Ndf
Definition: fopr_WilsonGeneral_impl.h:54
fopr_thread-inc.h
Imp::Fopr_WilsonGeneral::m_boundary_each_node
std::vector< double > m_boundary_each_node
b.c. on each node.
Definition: fopr_WilsonGeneral_impl.h:60
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Imp::Fopr_WilsonGeneral::H
void H(Field &v, const Field &w)
Definition: fopr_WilsonGeneral_impl.cpp:465
fopr_Wilson_impl_SU_N-inc.h
Imp::Fopr_WilsonGeneral::mult_tp_chiral
void mult_tp_chiral(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1697
fopr_Wilson_impl_SU3-inc.h
Imp::Fopr_WilsonGeneral::m_kappa_s
double m_kappa_s
Definition: fopr_WilsonGeneral_impl.h:45
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Imp::Fopr_WilsonGeneral::m_U
const Field_G * m_U
gauge configuration.
Definition: fopr_WilsonGeneral_impl.h:58
Imp::Fopr_WilsonGeneral::DDdag
void DDdag(Field &v, const Field &w)
Definition: fopr_WilsonGeneral_impl.cpp:455
Imp::Fopr_WilsonGeneral::m_Nz
int m_Nz
Definition: fopr_WilsonGeneral_impl.h:55
Parameters
Class for parameters.
Definition: parameters.h:46
Imp::Fopr_WilsonGeneral::m_Nvc
int m_Nvc
Definition: fopr_WilsonGeneral_impl.h:54
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
Imp::Fopr_WilsonGeneral::m_r_s
double m_r_s
Definition: fopr_WilsonGeneral_impl.h:46
Imp::Fopr_WilsonGeneral::m_Nvol
int m_Nvol
Definition: fopr_WilsonGeneral_impl.h:56
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
Imp::Fopr_WilsonGeneral::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_WilsonGeneral_impl.cpp:249
Imp::Fopr_WilsonGeneral::vcp2_xp
double * vcp2_xp
Definition: fopr_WilsonGeneral_impl.h:63
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Imp::Fopr_WilsonGeneral::vcp1_xp
double * vcp1_xp
arrays for communication buffer.
Definition: fopr_WilsonGeneral_impl.h:63
Imp::Fopr_WilsonGeneral::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_WilsonGeneral_impl.cpp:175
Imp::Fopr_WilsonGeneral::m_vl
Bridge::VerboseLevel m_vl
verbose level
Definition: fopr_WilsonGeneral_impl.h:49
Imp::Fopr_WilsonGeneral::mult_dag
void mult_dag(Field &v, const Field &w)
hermitian conjugate of mult.
Definition: fopr_WilsonGeneral_impl.cpp:303
Imp::Fopr_WilsonGeneral::mult_tm_dirac
void mult_tm_dirac(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1601
Imp::Fopr_WilsonGeneral::m_nu_s
double m_nu_s
Definition: fopr_WilsonGeneral_impl.h:46
Imp::Fopr_WilsonGeneral::class_name
static const std::string class_name
Definition: fopr_WilsonGeneral_impl.h:41
Imp::Fopr_WilsonGeneral::vcp2_zp
double * vcp2_zp
Definition: fopr_WilsonGeneral_impl.h:65
Imp::Fopr_WilsonGeneral::mult_dn
void mult_dn(const int mu, Field &, const Field &)
downward nearest neighbor hopping term.
Definition: fopr_WilsonGeneral_impl.cpp:391
Imp::Fopr_WilsonGeneral::D_dirac
void D_dirac(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_WilsonGeneral_impl.cpp:473
Imp::Fopr_WilsonGeneral::mult_xm
void mult_xm(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:747
Imp::Fopr_WilsonGeneral::setup
void setup()
initial setup main.
Definition: fopr_WilsonGeneral_impl.cpp:96
Imp::Fopr_WilsonGeneral::vcp2_xm
double * vcp2_xm
Definition: fopr_WilsonGeneral_impl.h:63
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:212
Imp::Fopr_WilsonGeneral::m_repr
std::string m_repr
gamma-matrix representation
Definition: fopr_WilsonGeneral_impl.h:48
Imp::Fopr_WilsonGeneral::m_w1
Field m_w1
Definition: fopr_WilsonGeneral_impl.h:68
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
Imp::Fopr_WilsonGeneral::vcp1_tp
double * vcp1_tp
Definition: fopr_WilsonGeneral_impl.h:66
Imp::Fopr_WilsonGeneral::mult_ym
void mult_ym(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1048
Imp::Fopr_WilsonGeneral::vcp1_zp
double * vcp1_zp
Definition: fopr_WilsonGeneral_impl.h:65
Imp::Fopr_WilsonGeneral::vcp1_tm
double * vcp1_tm
Definition: fopr_WilsonGeneral_impl.h:66
Imp::Fopr_WilsonGeneral::m_Nd
int m_Nd
Definition: fopr_WilsonGeneral_impl.h:54
Imp::Fopr_WilsonGeneral::Ddag
void Ddag(Field &v, const Field &w)
Definition: fopr_WilsonGeneral_impl.cpp:436
Imp::Fopr_WilsonGeneral::m_Nt
int m_Nt
Definition: fopr_WilsonGeneral_impl.h:55
Imp::Fopr_WilsonGeneral::D_chiral
void D_chiral(Field &, const int ex1, const Field &, const int ex2)
Definition: fopr_WilsonGeneral_impl.cpp:495
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Imp::Fopr_WilsonGeneral::m_kappa_t
double m_kappa_t
Definition: fopr_WilsonGeneral_impl.h:45
Imp::Fopr_WilsonGeneral::m_Nx
int m_Nx
Definition: fopr_WilsonGeneral_impl.h:55
Imp::Fopr_WilsonGeneral::mult_zm
void mult_zm(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1356
Imp::Fopr_WilsonGeneral::m_w2
Field m_w2
working fields
Definition: fopr_WilsonGeneral_impl.h:68
Imp::Fopr_WilsonGeneral::D
void D(Field &v, const Field &w)
Definition: fopr_WilsonGeneral_impl.cpp:425
Imp::Fopr_WilsonGeneral::mult_xp
void mult_xp(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:605
Imp::Fopr_WilsonGeneral::m_mode
std::string m_mode
mult mode
Definition: fopr_WilsonGeneral_impl.h:51
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Imp::Fopr_WilsonGeneral::DdagD
void DdagD(Field &v, const Field &w)
Definition: fopr_WilsonGeneral_impl.cpp:445
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
Imp::Fopr_WilsonGeneral::m_Ndim
int m_Ndim
Definition: fopr_WilsonGeneral_impl.h:56
Imp::Fopr_WilsonGeneral::vcp2_yp
double * vcp2_yp
Definition: fopr_WilsonGeneral_impl.h:64
fopr_Wilson_impl_common-inc.h
Imp::Fopr_WilsonGeneral::vcp1_xm
double * vcp1_xm
Definition: fopr_WilsonGeneral_impl.h:63
Imp::Fopr_WilsonGeneral::mult_yp
void mult_yp(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:899
Imp::Fopr_WilsonGeneral::m_Ny
int m_Ny
Definition: fopr_WilsonGeneral_impl.h:55
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
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_WilsonGeneral::mult_tm_chiral
void mult_tm_chiral(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1785
Imp::Fopr_WilsonGeneral::vcp2_tp
double * vcp2_tp
Definition: fopr_WilsonGeneral_impl.h:66
Imp::Fopr_WilsonGeneral::mult_gm5_chiral
void mult_gm5_chiral(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:584
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
Imp::Fopr_WilsonGeneral::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_WilsonGeneral_impl.cpp:272
Imp::Fopr_WilsonGeneral::mult_gm5_dirac
void mult_gm5_dirac(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:563
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
NVC
#define NVC
Definition: fopr_Wilson_impl_SU2-inc.h:15
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Imp::Fopr_WilsonGeneral::mult_gm5
void mult_gm5(Field &v, const Field &w)
multiplies gamma_5 matrix.
Definition: fopr_WilsonGeneral_impl.cpp:414
fopr_Wilson_impl_SU2-inc.h
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Imp::Fopr_WilsonGeneral::vcp1_yp
double * vcp1_yp
Definition: fopr_WilsonGeneral_impl.h:64
Imp::Fopr_WilsonGeneral::vcp1_ym
double * vcp1_ym
Definition: fopr_WilsonGeneral_impl.h:64
Imp::Fopr_WilsonGeneral::vcp2_tm
double * vcp2_tm
Definition: fopr_WilsonGeneral_impl.h:66
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Imp::Fopr_WilsonGeneral::daxpy
void daxpy(Field &, const double, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:539
Imp::Fopr_WilsonGeneral::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: fopr_WilsonGeneral_impl.cpp:263
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
Imp::Fopr_WilsonGeneral::flop_count
double flop_count()
returns the number of floating point operations.
Definition: fopr_WilsonGeneral_impl.cpp:1879
Imp::Fopr_WilsonGeneral::mult_zp
void mult_zp(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1205
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Imp
Clover term operator.
Definition: fopr_CloverTerm_eo_impl.cpp:32
Imp::Fopr_WilsonGeneral::mult_up
void mult_up(const int mu, Field &, const Field &)
upward nearest neighbor hopping term.
Definition: fopr_WilsonGeneral_impl.cpp:368
Imp::Fopr_WilsonGeneral::mult_tp_dirac
void mult_tp_dirac(Field &, const Field &)
Definition: fopr_WilsonGeneral_impl.cpp:1513
Field
Container of Field-type object.
Definition: field.h:46
Imp::Fopr_WilsonGeneral::clear
void clear(Field &)
Definition: fopr_WilsonGeneral_impl.cpp:517
Imp::Fopr_WilsonGeneral::vcp1_zm
double * vcp1_zm
Definition: fopr_WilsonGeneral_impl.h:65
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_WilsonGeneral::mult
void mult(Field &v, const Field &w)
multiplies fermion operator to a given field.
Definition: fopr_WilsonGeneral_impl.cpp:282
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Imp::Fopr_WilsonGeneral::m_boundary
std::vector< int > m_boundary
boundary condition
Definition: fopr_WilsonGeneral_impl.h:47
Imp::Fopr_WilsonGeneral::init
void init(const std::string repr)
obsolete initial setup.
Definition: fopr_WilsonGeneral_impl.cpp:73
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
Imp::Fopr_WilsonGeneral::vcp2_zm
double * vcp2_zm
Definition: fopr_WilsonGeneral_impl.h:65
Imp::Fopr_WilsonGeneral::m_Nc
int m_Nc
Definition: fopr_WilsonGeneral_impl.h:54
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