Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_Wilson_impl.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Wilson_impl.h"
15 
16 using std::valarray;
17 using std::string;
18 
19 //====================================================================
20 namespace {
21  inline double mult_uv_r(double *u, double *v)
22  {
23  return
24  u[0] * v[0] - u[1] * v[1]
25  + u[2] * v[2] - u[3] * v[3]
26  + u[4] * v[4] - u[5] * v[5];
27  }
28 
29 
30  inline double mult_uv_i(double *u, double *v)
31  {
32  return
33  u[0] * v[1] + u[1] * v[0]
34  + u[2] * v[3] + u[3] * v[2]
35  + u[4] * v[5] + u[5] * v[4];
36  }
37 
38 
39  inline double mult_udagv_r(double *u, double *v)
40  {
41  return
42  u[0] * v[0] + u[1] * v[1]
43  + u[6] * v[2] + u[7] * v[3]
44  + u[12] * v[4] + u[13] * v[5];
45  }
46 
47 
48  inline double mult_udagv_i(double *u, double *v)
49  {
50  return
51  u[0] * v[1] - u[1] * v[0]
52  + u[6] * v[3] - u[7] * v[2]
53  + u[12] * v[5] - u[13] * v[4];
54  }
55 }
56 
57 //====================================================================
58 void Fopr_Wilson::Fopr_Wilson_impl::set_parameters(const double kappa, const valarray<int> bc)
59 {
60  //- print input parameters
61  vout.general(m_vl, "Parameters of Fopr_Wilson_impl:\n");
62  vout.general(m_vl, " kappa = %8.4f\n", kappa);
63  for (int mu = 0; mu < m_Ndim; ++mu) {
64  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
65  }
66 
67  //- range check
68  // NB. kappa = 0 is allowed.
69  assert(bc.size() == m_Ndim);
70 
71  //- store values
72  m_kappa = kappa;
73 
74  // m_boundary.resize(m_Ndim); // already resized in init.
75  for (int mu = 0; mu < m_Ndim; ++mu) {
76  m_boundary[mu] = bc[mu];
77  }
78 }
79 
80 
81 //====================================================================
82 void Fopr_Wilson::Fopr_Wilson_impl::init(string repr)
83 {
84  m_vl = CommonParameters::Vlevel();
85 
86  m_Nvol = CommonParameters::Nvol();
87  m_Ndim = CommonParameters::Ndim();
88  m_boundary.resize(m_Ndim);
89 
90  m_U = 0;
91 
92  m_repr = repr;
93 
94  m_GM.resize(m_Ndim + 1);
95 
96  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);
97 
98  m_GM[0] = gmset->get_GM(GammaMatrixSet::GAMMA1);
99  m_GM[1] = gmset->get_GM(GammaMatrixSet::GAMMA2);
100  m_GM[2] = gmset->get_GM(GammaMatrixSet::GAMMA3);
101  m_GM[3] = gmset->get_GM(GammaMatrixSet::GAMMA4);
102  m_GM[4] = gmset->get_GM(GammaMatrixSet::GAMMA5);
103 
104 
105  m_mult = &Fopr_Wilson::Fopr_Wilson_impl::mult_undef;
106  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::mult_undef;
107 
108  if (m_repr == "Dirac") {
109  m_D = &Fopr_Wilson::Fopr_Wilson_impl::D_dirac;
110  m_gm5 = &Fopr_Wilson::Fopr_Wilson_impl::gm5_dirac;
111  m_mult_tp = &Fopr_Wilson::Fopr_Wilson_impl::mult_tp_dirac;
112  m_mult_tm = &Fopr_Wilson::Fopr_Wilson_impl::mult_tm_dirac;
113  } else if (m_repr == "Chiral") {
114  m_D = &Fopr_Wilson::Fopr_Wilson_impl::D_chiral;
115  m_gm5 = &Fopr_Wilson::Fopr_Wilson_impl::gm5_chiral;
116  m_mult_tp = &Fopr_Wilson::Fopr_Wilson_impl::mult_tp_chiral;
117  m_mult_tm = &Fopr_Wilson::Fopr_Wilson_impl::mult_tm_chiral;
118  } else {
119  vout.crucial(m_vl, "Fopr_Wilson: input repr is undefined.\n");
120  abort();
121  }
122 
123  delete gmset;
124 }
125 
126 
127 //====================================================================
128 void Fopr_Wilson::Fopr_Wilson_impl::set_mode(std::string mode)
129 {
130  m_mode = mode;
131 
132  if (m_mode == "D") {
133  m_mult = &Fopr_Wilson::Fopr_Wilson_impl::D;
134  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::Ddag;
135  } else if (m_mode == "Ddag") {
136  m_mult = &Fopr_Wilson::Fopr_Wilson_impl::Ddag;
137  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::D;
138  } else if (m_mode == "DdagD") {
139  m_mult = &Fopr_Wilson::Fopr_Wilson_impl::DdagD;
140  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::DdagD;
141  } else if (m_mode == "H") {
142  m_mult = &Fopr_Wilson::Fopr_Wilson_impl::H;
143  m_mult_dag = &Fopr_Wilson::Fopr_Wilson_impl::H;
144  } else {
145  vout.crucial(m_vl, "Fopr_Wilson: input mode is undefined.\n");
146  abort();
147  }
148 }
149 
150 
151 //====================================================================
152 std::string Fopr_Wilson::Fopr_Wilson_impl::get_mode() const
153 {
154  return m_mode;
155 }
156 
157 
158 //====================================================================
159 void Fopr_Wilson::Fopr_Wilson_impl::D_dirac(Field& w, const Field& f)
160 {
161  w = 0.0;
162 
163  mult_xp(w, f);
164  mult_xm(w, f);
165 
166  mult_yp(w, f);
167  mult_ym(w, f);
168 
169  mult_zp(w, f);
170  mult_zm(w, f);
171 
172  mult_tp_dirac(w, f);
173  mult_tm_dirac(w, f);
174 
175  w *= -m_kappa;
176  w += f;
177 }
178 
179 
180 //====================================================================
181 void Fopr_Wilson::Fopr_Wilson_impl::D_chiral(Field& w, const Field& f)
182 {
183  mult_xp(w, f);
184  mult_xm(w, f);
185 
186  mult_yp(w, f);
187  mult_ym(w, f);
188 
189  mult_zp(w, f);
190  mult_zm(w, f);
191 
192  mult_tp_chiral(w, f);
193  mult_tm_chiral(w, f);
194 
195  w *= -m_kappa;
196  w += f;
197 }
198 
199 
200 //====================================================================
201 void Fopr_Wilson::Fopr_Wilson_impl::mult_up(int mu, Field& w, const Field& f)
202 {
203  if (mu == 0) {
204  mult_xp(w, f);
205  } else if (mu == 1) {
206  mult_yp(w, f);
207  } else if (mu == 2) {
208  mult_zp(w, f);
209  } else if (mu == 3) {
210  (this->*m_mult_tp)(w, f);
211  } else {
212  vout.crucial(m_vl, "Fopr_Wilson: illegal parameter mu in mult_up.\n");
213  abort();
214  }
215 }
216 
217 
218 //====================================================================
219 void Fopr_Wilson::Fopr_Wilson_impl::mult_dn(int mu, Field& w, const Field& f)
220 {
221  if (mu == 0) {
222  mult_xm(w, f);
223  } else if (mu == 1) {
224  mult_ym(w, f);
225  } else if (mu == 2) {
226  mult_zm(w, f);
227  } else if (mu == 3) {
228  (this->*m_mult_tm)(w, f);
229  } else {
230  vout.crucial(m_vl, "Fopr_Wilson: illegal parameter mu in mult_dn.\n");
231  abort();
232  }
233 }
234 
235 
236 //====================================================================
237 void Fopr_Wilson::Fopr_Wilson_impl::gm5_dirac(Field& w, const Field& f)
238 {
239  int Nvc = 2 * CommonParameters::Nc();
240  int Nd = CommonParameters::Nd();
241 
242  double *v1;
243  double *v2;
244 
245  v1 = const_cast<Field *>(&f)->ptr(0);
246  v2 = w.ptr(0);
247 
248  int id1 = 0;
249  int id2 = Nvc;
250  int id3 = Nvc * 2;
251  int id4 = Nvc * 3;
252 
253  for (int site = 0; site < m_Nvol; ++site) {
254  for (int icc = 0; icc < Nvc; icc++) {
255  int in = Nvc * Nd * site;
256  v2[icc + id1 + in] = v1[icc + id3 + in];
257  v2[icc + id2 + in] = v1[icc + id4 + in];
258  v2[icc + id3 + in] = v1[icc + id1 + in];
259  v2[icc + id4 + in] = v1[icc + id2 + in];
260  }
261  }
262 }
263 
264 
265 //====================================================================
266 void Fopr_Wilson::Fopr_Wilson_impl::gm5_chiral(Field& w, const Field& f)
267 {
268  int Nvc = 2 * CommonParameters::Nc();
269  int Nd = CommonParameters::Nd();
270 
271  double *v1;
272  double *v2;
273 
274  v1 = const_cast<Field *>(&f)->ptr(0);
275  v2 = w.ptr(0);
276 
277  int id1 = 0;
278  int id2 = Nvc;
279  int id3 = Nvc * 2;
280  int id4 = Nvc * 3;
281 
282  for (int site = 0; site < m_Nvol; ++site) {
283  for (int icc = 0; icc < Nvc; icc++) {
284  int in = Nvc * Nd * site;
285  v2[icc + id1 + in] = v1[icc + id1 + in];
286  v2[icc + id2 + in] = v1[icc + id2 + in];
287  v2[icc + id3 + in] = -v1[icc + id3 + in];
288  v2[icc + id4 + in] = -v1[icc + id4 + in];
289  }
290  }
291 }
292 
293 
294 //====================================================================
295 const Field_F Fopr_Wilson::Fopr_Wilson_impl::mult_gm5p(int mu, const Field_F& w)
296 {
297  Field vt(w.nin(), w.nvol(), w.nex());
298 
299  vt = 0.0;
300 
301  if (mu == 0) {
302  mult_xp(vt, (Field)w);
303  } else if (mu == 1) {
304  mult_yp(vt, (Field)w);
305  } else if (mu == 2) {
306  mult_zp(vt, (Field)w);
307  } else if (mu == 3) {
308  if (m_repr == "Dirac") {
309  mult_tp_dirac(vt, w);
310  } else {
311  mult_tp_chiral(vt, w);
312  }
313  } else {
314  abort();
315  }
316 
317  // return (Field_F) mult_gm5(vt);
318 
319  Field vt2(w.nin(), w.nvol(), w.nex());
320  mult_gm5(vt2, vt);
321 
322  return vt2;
323 }
324 
325 
326 //====================================================================
327 void Fopr_Wilson::Fopr_Wilson_impl::mult_xp(Field& w, const Field& f)
328 {
329  int Ncol = CommonParameters::Nc();
330  int Nvc = 2 * Ncol;
331  int Ndf = 2 * Ncol * Ncol;
332  int ND = CommonParameters::Nd();
333  int Nx = CommonParameters::Nx();
334  int Ny = CommonParameters::Ny();
335  int Nz = CommonParameters::Nz();
336  int Nt = CommonParameters::Nt();
337  int Nst = Nx * Ny * Nz * Nt;
338 
339  double *v1;
340  double *v2;
341  double *u;
342 
343  v2 = w.ptr(0);
344  v1 = const_cast<Field *>(&f)->ptr(0);
345  u = const_cast<Field_G *>(m_U)->ptr(0);
346 
347  int idir = 0;
348 
349  int id1 = 0;
350  int id2 = Nvc;
351  int id3 = Nvc * 2;
352  int id4 = Nvc * 3;
353 
354  int ix, nn;
355 
356  double vt1[Nvc], vt2[Nvc];
357  double wt1r, wt1i, wt2r, wt2i;
358  double bc2 = 1.0;
359  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
360 
361  double vcp1[Nvc * 2 * Ny * Nz * Nt], vcp2[Nvc * 2 * Ny * Nz * Nt];
362 
363  //- boundary part
364  ix = Nx - 1;
365  nn = 0;
366  for (int iyzt = 0; iyzt < (Ny * Nz * Nt); iyzt++) {
367  int in = Nvc * ND * (nn + iyzt * Nx);
368  int ix1 = Nvc * 2 * iyzt;
369  int ix2 = ix1 + Nvc;
370 
371  for (int ic = 0; ic < Ncol; ic++) {
372  vcp1[2 * ic + ix1] = bc2 * (v1[2 * ic + id1 + in] - v1[2 * ic + 1 + id4 + in]);
373  vcp1[2 * ic + 1 + ix1] = bc2 * (v1[2 * ic + 1 + id1 + in] + v1[2 * ic + id4 + in]);
374  vcp1[2 * ic + ix2] = bc2 * (v1[2 * ic + id2 + in] - v1[2 * ic + 1 + id3 + in]);
375  vcp1[2 * ic + 1 + ix2] = bc2 * (v1[2 * ic + 1 + id2 + in] + v1[2 * ic + id3 + in]);
376  }
377  }
378 
379  Communicator::exchange(Nvc * 2 * Ny * Nz * Nt, vcp2, vcp1, 0, 1, 1);
380 
381  ix = Nx - 1;
382  nn = 0;
383  for (int iyzt = 0; iyzt < (Ny * Nz * Nt); iyzt++) {
384  int iv = Nvc * ND * (ix + iyzt * Nx);
385  int ig = Ndf * (ix + iyzt * Nx + idir * Nst);
386  int ix1 = Nvc * 2 * iyzt;
387  int ix2 = ix1 + Nvc;
388 
389  for (int ic = 0; ic < Ncol; ic++) {
390  int ic2 = ic * Nvc;
391 
392  wt1r = mult_uv_r(&u[ic2 + ig], &vcp2[ix1]);
393  wt1i = mult_uv_i(&u[ic2 + ig], &vcp2[ix1]);
394  wt2r = mult_uv_r(&u[ic2 + ig], &vcp2[ix2]);
395  wt2i = mult_uv_i(&u[ic2 + ig], &vcp2[ix2]);
396 
397  v2[2 * ic + id1 + iv] = wt1r;
398  v2[2 * ic + 1 + id1 + iv] = wt1i;
399  v2[2 * ic + id2 + iv] = wt2r;
400  v2[2 * ic + 1 + id2 + iv] = wt2i;
401  v2[2 * ic + id3 + iv] = wt2i;
402  v2[2 * ic + 1 + id3 + iv] = -wt2r;
403  v2[2 * ic + id4 + iv] = wt1i;
404  v2[2 * ic + 1 + id4 + iv] = -wt1r;
405  }
406  }
407 
408  //- bulk part
409  for (int iyzt = 0; iyzt < (Ny * Nz * Nt); iyzt++) {
410  for (int ix = 0; ix < (Nx - 1); ix++) {
411  nn = ix + 1;
412 
413  int iv = Nvc * ND * (ix + iyzt * Nx);
414  int in = Nvc * ND * (nn + iyzt * Nx);
415  int ig = Ndf * (ix + iyzt * Nx + idir * Nst);
416 
417  for (int ic = 0; ic < Ncol; ic++) {
418  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + 1 + id4 + in];
419  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] + v1[2 * ic + id4 + in];
420  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + 1 + id3 + in];
421  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + id3 + in];
422  }
423 
424  for (int ic = 0; ic < Ncol; ic++) {
425  int ic2 = ic * Nvc;
426 
427  wt1r = mult_uv_r(&u[ic2 + ig], vt1);
428  wt1i = mult_uv_i(&u[ic2 + ig], vt1);
429  wt2r = mult_uv_r(&u[ic2 + ig], vt2);
430  wt2i = mult_uv_i(&u[ic2 + ig], vt2);
431 
432  v2[2 * ic + id1 + iv] = wt1r;
433  v2[2 * ic + 1 + id1 + iv] = wt1i;
434  v2[2 * ic + id2 + iv] = wt2r;
435  v2[2 * ic + 1 + id2 + iv] = wt2i;
436  v2[2 * ic + id3 + iv] = wt2i;
437  v2[2 * ic + 1 + id3 + iv] = -wt2r;
438  v2[2 * ic + id4 + iv] = wt1i;
439  v2[2 * ic + 1 + id4 + iv] = -wt1r;
440  }
441  }
442  }
443 }
444 
445 
446 //====================================================================
447 void Fopr_Wilson::Fopr_Wilson_impl::mult_xm(Field& w, const Field& f)
448 {
449  int Ncol = CommonParameters::Nc();
450  int Nvc = 2 * Ncol;
451  int Ndf = 2 * Ncol * Ncol;
452  int ND = CommonParameters::Nd();
453  int Nx = CommonParameters::Nx();
454  int Ny = CommonParameters::Ny();
455  int Nz = CommonParameters::Nz();
456  int Nt = CommonParameters::Nt();
457  int Nst = Nx * Ny * Nz * Nt;
458 
459  double *v1;
460  double *v2;
461  double *u;
462 
463  v2 = w.ptr(0);
464  v1 = const_cast<Field *>(&f)->ptr(0);
465  u = const_cast<Field_G *>(m_U)->ptr(0);
466 
467  int idir = 0;
468 
469  int id1 = 0;
470  int id2 = Nvc;
471  int id3 = Nvc * 2;
472  int id4 = Nvc * 3;
473 
474  int ix, nn;
475 
476  double vt1[Nvc], vt2[Nvc];
477  double wt1r, wt1i, wt2r, wt2i;
478  double bc2 = 1.0;
479  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
480 
481  double vcp1[Nvc * 2 * Ny * Nz * Nt], vcp2[Nvc * 2 * Ny * Nz * Nt];
482 
483  //- boundary part
484  ix = 0;
485  nn = Nx - 1;
486  for (int iyzt = 0; iyzt < (Ny * Nz * Nt); iyzt++) {
487  int in = Nvc * ND * (nn + iyzt * Nx);
488  int ig = Ndf * (nn + iyzt * Nx + idir * Nst);
489  int ix1 = Nvc * 2 * iyzt;
490  int ix2 = ix1 + Nvc;
491 
492  for (int ic = 0; ic < Ncol; ic++) {
493  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + 1 + id4 + in];
494  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + id4 + in];
495 
496  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + 1 + id3 + in];
497  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + id3 + in];
498  }
499 
500  for (int ic = 0; ic < Ncol; ic++) {
501  int icr = 2 * ic;
502  vcp1[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1);
503  vcp1[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1);
504  vcp1[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2);
505  vcp1[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2);
506  }
507  }
508 
509  Communicator::exchange(Nvc * 2 * Ny * Nz * Nt, vcp2, vcp1, 0, -1, 2);
510 
511  ix = 0;
512  nn = Nx - 1;
513  for (int iyzt = 0; iyzt < (Ny * Nz * Nt); iyzt++) {
514  int iv = Nvc * ND * (ix + iyzt * Nx);
515  int ix1 = Nvc * 2 * iyzt;
516  int ix2 = ix1 + Nvc;
517 
518  for (int ic = 0; ic < Ncol; ic++) {
519  int icr = 2 * ic;
520  int ici = 2 * ic + 1;
521  v2[icr + id1 + iv] += bc2 * vcp2[icr + ix1];
522  v2[ici + id1 + iv] += bc2 * vcp2[ici + ix1];
523  v2[icr + id2 + iv] += bc2 * vcp2[icr + ix2];
524  v2[ici + id2 + iv] += bc2 * vcp2[ici + ix2];
525  v2[icr + id3 + iv] += -bc2 * vcp2[ici + ix2];
526  v2[ici + id3 + iv] += +bc2 * vcp2[icr + ix2];
527  v2[icr + id4 + iv] += -bc2 * vcp2[ici + ix1];
528  v2[ici + id4 + iv] += +bc2 * vcp2[icr + ix1];
529  }
530  }
531 
532  //- bulk part
533  for (int iyzt = 0; iyzt < Ny * Nz * Nt; iyzt++) {
534  for (int ix = 1; ix < Nx; ix++) {
535  nn = ix - 1;
536 
537  int iv = Nvc * ND * (ix + (iyzt) * Nx);
538  int in = Nvc * ND * (nn + (iyzt) * Nx);
539  int ig = Ndf * (nn + iyzt * Nx + idir * Nst);
540 
541  for (int ic = 0; ic < Ncol; ic++) {
542  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + 1 + id4 + in];
543  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + id4 + in];
544 
545  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + 1 + id3 + in];
546  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + id3 + in];
547  }
548 
549  for (int ic = 0; ic < Ncol; ic++) {
550  int icr = 2 * ic;
551  int ici = 2 * ic + 1;
552  int ic1 = 0;
553  int ic2 = Nvc;
554  int ic3 = 2 * Nvc;
555 
556  wt1r = mult_udagv_r(&u[icr + ig], vt1);
557  wt1i = mult_udagv_i(&u[icr + ig], vt1);
558  wt2r = mult_udagv_r(&u[icr + ig], vt2);
559  wt2i = mult_udagv_i(&u[icr + ig], vt2);
560 
561  v2[icr + id1 + iv] += wt1r;
562  v2[ici + id1 + iv] += wt1i;
563  v2[icr + id2 + iv] += wt2r;
564  v2[ici + id2 + iv] += wt2i;
565  v2[icr + id3 + iv] += -wt2i;
566  v2[ici + id3 + iv] += +wt2r;
567  v2[icr + id4 + iv] += -wt1i;
568  v2[ici + id4 + iv] += +wt1r;
569  }
570  }
571  }
572 }
573 
574 
575 //====================================================================
576 void Fopr_Wilson::Fopr_Wilson_impl::mult_yp(Field& w, const Field& f)
577 {
578  int Ncol = CommonParameters::Nc();
579  int Nvc = 2 * Ncol;
580  int Ndf = 2 * Ncol * Ncol;
581  int ND = CommonParameters::Nd();
582  int Nx = CommonParameters::Nx();
583  int Ny = CommonParameters::Ny();
584  int Nz = CommonParameters::Nz();
585  int Nt = CommonParameters::Nt();
586  int Nst = Nx * Ny * Nz * Nt;
587 
588  double *v1;
589  double *v2;
590  double *u;
591 
592  v2 = w.ptr(0);
593  v1 = const_cast<Field *>(&f)->ptr(0);
594  u = const_cast<Field_G *>(m_U)->ptr(0);
595 
596  int idir = 1;
597 
598  int id1 = 0;
599  int id2 = Nvc;
600  int id3 = Nvc * 2;
601  int id4 = Nvc * 3;
602 
603  double vt1[Nvc], vt2[Nvc];
604  double wt1r, wt1i, wt2r, wt2i;
605  double bc2 = 1.0;
606  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
607 
608  int iy, nn;
609  double vcp1[Nvc * 2 * Nx * Nz * Nt], vcp2[Nvc * 2 * Nx * Nz * Nt];
610 
611  //- boundary part
612  iy = Ny - 1;
613  nn = 0;
614  for (int izt = 0; izt < (Nz * Nt); izt++) {
615  for (int ix = 0; ix < Nx; ix++) {
616  int in = Nvc * ND * (ix + nn * Nx + izt * Nx * Ny);
617  int ix1 = Nvc * 2 * (ix + izt * Nx);
618  int ix2 = ix1 + Nvc;
619 
620  for (int ic = 0; ic < Ncol; ic++) {
621  vcp1[2 * ic + ix1] = bc2 * (v1[2 * ic + id1 + in] + v1[2 * ic + id4 + in]);
622  vcp1[2 * ic + 1 + ix1] = bc2 * (v1[2 * ic + 1 + id1 + in] + v1[2 * ic + 1 + id4 + in]);
623  vcp1[2 * ic + ix2] = bc2 * (v1[2 * ic + id2 + in] - v1[2 * ic + id3 + in]);
624  vcp1[2 * ic + 1 + ix2] = bc2 * (v1[2 * ic + 1 + id2 + in] - v1[2 * ic + 1 + id3 + in]);
625  }
626  }
627  }
628 
629  Communicator::exchange(Nvc * 2 * Nx * Nz * Nt, vcp2, vcp1, 1, 1, 3);
630 
631  iy = Ny - 1;
632  nn = 0;
633  for (int izt = 0; izt < (Nz * Nt); izt++) {
634  for (int ix = 0; ix < Nx; ix++) {
635  int iv = Nvc * ND * (ix + iy * Nx + izt * Nx * Ny);
636  int ig = Ndf * (ix + iy * Nx + izt * Nx * Ny + idir * Nst);
637  int ix1 = Nvc * 2 * (ix + izt * Nx);
638  int ix2 = ix1 + Nvc;
639 
640  for (int ic = 0; ic < Ncol; ic++) {
641  int ic2 = ic * Nvc;
642 
643  wt1r = mult_uv_r(&u[ic2 + ig], &vcp2[ix1]);
644  wt1i = mult_uv_i(&u[ic2 + ig], &vcp2[ix1]);
645  wt2r = mult_uv_r(&u[ic2 + ig], &vcp2[ix2]);
646  wt2i = mult_uv_i(&u[ic2 + ig], &vcp2[ix2]);
647 
648  v2[2 * ic + id1 + iv] += wt1r;
649  v2[2 * ic + 1 + id1 + iv] += wt1i;
650  v2[2 * ic + id2 + iv] += wt2r;
651  v2[2 * ic + 1 + id2 + iv] += wt2i;
652  v2[2 * ic + id3 + iv] += -wt2r;
653  v2[2 * ic + 1 + id3 + iv] += -wt2i;
654  v2[2 * ic + id4 + iv] += wt1r;
655  v2[2 * ic + 1 + id4 + iv] += wt1i;
656  }
657  }
658  }
659 
660  //- bulk part
661  for (int izt = 0; izt < (Nz * Nt); izt++) {
662  for (int iy = 0; iy < (Ny - 1); iy++) {
663  int nn = iy + 1;
664  for (int ix = 0; ix < Nx; ix++) {
665  int iv = Nvc * ND * (ix + iy * Nx + izt * Nx * Ny);
666  int in = Nvc * ND * (ix + nn * Nx + izt * Nx * Ny);
667  int ig = Ndf * (ix + iy * Nx + izt * Nx * Ny + idir * Nst);
668 
669  for (int ic = 0; ic < Ncol; ic++) {
670  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + id4 + in];
671  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] + v1[2 * ic + 1 + id4 + in];
672  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + id3 + in];
673  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + 1 + id3 + in];
674  }
675 
676  for (int ic = 0; ic < Ncol; ic++) {
677  int ic2 = ic * Nvc;
678 
679  wt1r = mult_uv_r(&u[ic2 + ig], vt1);
680  wt1i = mult_uv_i(&u[ic2 + ig], vt1);
681  wt2r = mult_uv_r(&u[ic2 + ig], vt2);
682  wt2i = mult_uv_i(&u[ic2 + ig], vt2);
683 
684  v2[2 * ic + id1 + iv] += wt1r;
685  v2[2 * ic + 1 + id1 + iv] += wt1i;
686  v2[2 * ic + id2 + iv] += wt2r;
687  v2[2 * ic + 1 + id2 + iv] += wt2i;
688  v2[2 * ic + id3 + iv] += -wt2r;
689  v2[2 * ic + 1 + id3 + iv] += -wt2i;
690  v2[2 * ic + id4 + iv] += wt1r;
691  v2[2 * ic + 1 + id4 + iv] += wt1i;
692  }
693  }
694  }
695  }
696 }
697 
698 
699 //====================================================================
700 void Fopr_Wilson::Fopr_Wilson_impl::mult_ym(Field& w, const Field& f)
701 {
702  int Ncol = CommonParameters::Nc();
703  int Nvc = 2 * Ncol;
704  int Ndf = 2 * Ncol * Ncol;
705  int ND = CommonParameters::Nd();
706  int Nx = CommonParameters::Nx();
707  int Ny = CommonParameters::Ny();
708  int Nz = CommonParameters::Nz();
709  int Nt = CommonParameters::Nt();
710  int Nst = Nx * Ny * Nz * Nt;
711 
712  double *v1;
713  double *v2;
714  double *u;
715 
716  v2 = w.ptr(0);
717  v1 = const_cast<Field *>(&f)->ptr(0);
718  u = const_cast<Field_G *>(m_U)->ptr(0);
719 
720  int idir = 1;
721 
722  int id1 = 0;
723  int id2 = Nvc;
724  int id3 = Nvc * 2;
725  int id4 = Nvc * 3;
726 
727  double vt1[Nvc], vt2[Nvc];
728  double wt1r, wt1i, wt2r, wt2i;
729  double bc2 = 1.0;
730  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
731 
732  int iy, nn;
733  double vcp1[Nvc * 2 * Nx * Nz * Nt], vcp2[Nvc * 2 * Nx * Nz * Nt];
734 
735  //- boundary part
736  iy = 0;
737  nn = Ny - 1;
738  for (int izt = 0; izt < (Nz * Nt); izt++) {
739  for (int ix = 0; ix < Nx; ix++) {
740  int in = Nvc * ND * (ix + nn * Nx + izt * Nx * Ny);
741  int ig = Ndf * (ix + nn * Nx + izt * Nx * Ny + idir * Nst);
742  int ix1 = Nvc * 2 * (ix + izt * Nx);
743  int ix2 = ix1 + Nvc;
744 
745  for (int ic = 0; ic < Ncol; ic++) {
746  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + id4 + in];
747  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + 1 + id4 + in];
748 
749  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + id3 + in];
750  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + 1 + id3 + in];
751  }
752 
753  for (int ic = 0; ic < Ncol; ic++) {
754  int icr = 2 * ic;
755  int ici = 2 * ic + 1;
756  vcp1[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1);
757  vcp1[ici + ix1] = mult_udagv_i(&u[icr + ig], vt1);
758  vcp1[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2);
759  vcp1[ici + ix2] = mult_udagv_i(&u[icr + ig], vt2);
760  }
761  }
762  }
763 
764  Communicator::exchange(Nvc * 2 * Nx * Nz * Nt, vcp2, vcp1, 1, -1, 4);
765 
766  iy = 0;
767  nn = Ny - 1;
768  for (int izt = 0; izt < (Nz * Nt); izt++) {
769  for (int ix = 0; ix < Nx; ix++) {
770  int iv = Nvc * ND * (ix + iy * Nx + izt * Nx * Ny);
771  int ix1 = Nvc * 2 * (ix + izt * Nx);
772  int ix2 = ix1 + Nvc;
773 
774  for (int ic = 0; ic < Ncol; ic++) {
775  int icr = 2 * ic;
776  int ici = 2 * ic + 1;
777  v2[icr + id1 + iv] += bc2 * vcp2[icr + ix1];
778  v2[ici + id1 + iv] += bc2 * vcp2[ici + ix1];
779  v2[icr + id2 + iv] += bc2 * vcp2[icr + ix2];
780  v2[ici + id2 + iv] += bc2 * vcp2[ici + ix2];
781  v2[icr + id3 + iv] += bc2 * vcp2[icr + ix2];
782  v2[ici + id3 + iv] += bc2 * vcp2[ici + ix2];
783  v2[icr + id4 + iv] += -bc2 * vcp2[icr + ix1];
784  v2[ici + id4 + iv] += -bc2 * vcp2[ici + ix1];
785  }
786  }
787  }
788 
789  //- bulk part
790  for (int izt = 0; izt < (Nz * Nt); izt++) {
791  for (int iy = 1; iy < Ny; iy++) {
792  int nn = iy - 1;
793  for (int ix = 0; ix < Nx; ix++) {
794  int iv = Nvc * ND * (ix + iy * Nx + izt * Nx * Ny);
795  int in = Nvc * ND * (ix + nn * Nx + izt * Nx * Ny);
796  int ig = Ndf * (ix + nn * Nx + izt * Nx * Ny + idir * Nst);
797 
798  for (int ic = 0; ic < Ncol; ic++) {
799  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + id4 + in];
800  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + 1 + id4 + in];
801 
802  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + id3 + in];
803  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + 1 + id3 + in];
804  }
805 
806  for (int ic = 0; ic < Ncol; ic++) {
807  int icr = 2 * ic;
808  int ici = 2 * ic + 1;
809  int ic1 = 0;
810  int ic2 = Nvc;
811  int ic3 = 2 * Nvc;
812 
813  wt1r = mult_udagv_r(&u[icr + ig], vt1);
814  wt1i = mult_udagv_i(&u[icr + ig], vt1);
815  wt2r = mult_udagv_r(&u[icr + ig], vt2);
816  wt2i = mult_udagv_i(&u[icr + ig], vt2);
817 
818  v2[icr + id1 + iv] += wt1r;
819  v2[ici + id1 + iv] += wt1i;
820  v2[icr + id2 + iv] += wt2r;
821  v2[ici + id2 + iv] += wt2i;
822  v2[icr + id3 + iv] += wt2r;
823  v2[ici + id3 + iv] += wt2i;
824  v2[icr + id4 + iv] += -wt1r;
825  v2[ici + id4 + iv] += -wt1i;
826  }
827  }
828  }
829  }
830 }
831 
832 
833 //====================================================================
834 void Fopr_Wilson::Fopr_Wilson_impl::mult_zp(Field& w, const Field& f)
835 {
836  int Ncol = CommonParameters::Nc();
837  int Nvc = 2 * Ncol;
838  int Ndf = 2 * Ncol * Ncol;
839  int ND = CommonParameters::Nd();
840  int Nx = CommonParameters::Nx();
841  int Ny = CommonParameters::Ny();
842  int Nz = CommonParameters::Nz();
843  int Nt = CommonParameters::Nt();
844  int Nst = Nx * Ny * Nz * Nt;
845 
846  double *v1;
847  double *v2;
848  double *u;
849 
850  v2 = w.ptr(0);
851  v1 = const_cast<Field *>(&f)->ptr(0);
852  u = const_cast<Field_G *>(m_U)->ptr(0);
853 
854  int idir = 2;
855 
856  int id1 = 0;
857  int id2 = Nvc;
858  int id3 = Nvc * 2;
859  int id4 = Nvc * 3;
860 
861  double vt1[Nvc], vt2[Nvc];
862  double wt1r, wt1i, wt2r, wt2i;
863  double bc2 = 1.0;
864  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
865 
866  int iz, nn;
867  double vcp1[Nvc * 2 * Nx * Ny * Nt], vcp2[Nvc * 2 * Nx * Ny * Nt];
868 
869  //- boundary part
870  iz = Nz - 1;
871  nn = 0;
872  for (int it = 0; it < Nt; it++) {
873  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
874  int in = Nvc * ND * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz);
875  int ix1 = Nvc * 2 * (ixy + it * Nx * Ny);
876  int ix2 = ix1 + Nvc;
877 
878  for (int ic = 0; ic < Ncol; ic++) {
879  vcp1[2 * ic + ix1] = bc2 * (v1[2 * ic + id1 + in] - v1[2 * ic + 1 + id3 + in]);
880  vcp1[2 * ic + 1 + ix1] = bc2 * (v1[2 * ic + 1 + id1 + in] + v1[2 * ic + id3 + in]);
881  vcp1[2 * ic + ix2] = bc2 * (v1[2 * ic + id2 + in] + v1[2 * ic + 1 + id4 + in]);
882  vcp1[2 * ic + 1 + ix2] = bc2 * (v1[2 * ic + 1 + id2 + in] - v1[2 * ic + id4 + in]);
883  }
884  }
885  }
886 
887  Communicator::exchange(Nvc * 2 * Nx * Ny * Nt, vcp2, vcp1, 2, 1, 5);
888 
889  iz = Nz - 1;
890  nn = 0;
891  for (int it = 0; it < Nt; it++) {
892  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
893  int iv = Nvc * ND * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz);
894  int ig = Ndf * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz + idir * Nst);
895  int ix1 = Nvc * 2 * (ixy + it * Nx * Ny);
896  int ix2 = ix1 + Nvc;
897 
898  for (int ic = 0; ic < Ncol; ic++) {
899  int ic2 = ic * Nvc;
900 
901  wt1r = mult_uv_r(&u[ic2 + ig], &vcp2[ix1]);
902  wt1i = mult_uv_i(&u[ic2 + ig], &vcp2[ix1]);
903  wt2r = mult_uv_r(&u[ic2 + ig], &vcp2[ix2]);
904  wt2i = mult_uv_i(&u[ic2 + ig], &vcp2[ix2]);
905 
906  v2[2 * ic + id1 + iv] += wt1r;
907  v2[2 * ic + 1 + id1 + iv] += wt1i;
908  v2[2 * ic + id2 + iv] += wt2r;
909  v2[2 * ic + 1 + id2 + iv] += wt2i;
910  v2[2 * ic + id3 + iv] += wt1i;
911  v2[2 * ic + 1 + id3 + iv] += -wt1r;
912  v2[2 * ic + id4 + iv] += -wt2i;
913  v2[2 * ic + 1 + id4 + iv] += wt2r;
914  }
915  }
916  }
917 
918  //- bulk part
919  for (int it = 0; it < Nt; it++) {
920  for (int iz = 0; iz < (Nz - 1); iz++) {
921  int nn = iz + 1;
922  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
923  int iv = Nvc * ND * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz);
924  int in = Nvc * ND * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz);
925  int ig = Ndf * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz + idir * Nst);
926 
927  for (int ic = 0; ic < Ncol; ic++) {
928  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + 1 + id3 + in];
929  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] + v1[2 * ic + id3 + in];
930  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + 1 + id4 + in];
931  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + id4 + in];
932  }
933 
934  for (int ic = 0; ic < Ncol; ic++) {
935  int ic2 = ic * Nvc;
936 
937  wt1r = mult_uv_r(&u[ic2 + ig], vt1);
938  wt1i = mult_uv_i(&u[ic2 + ig], vt1);
939  wt2r = mult_uv_r(&u[ic2 + ig], vt2);
940  wt2i = mult_uv_i(&u[ic2 + ig], vt2);
941 
942  v2[2 * ic + id1 + iv] += wt1r;
943  v2[2 * ic + 1 + id1 + iv] += wt1i;
944  v2[2 * ic + id2 + iv] += wt2r;
945  v2[2 * ic + 1 + id2 + iv] += wt2i;
946  v2[2 * ic + id3 + iv] += wt1i;
947  v2[2 * ic + 1 + id3 + iv] += -wt1r;
948  v2[2 * ic + id4 + iv] += -wt2i;
949  v2[2 * ic + 1 + id4 + iv] += wt2r;
950  }
951  }
952  }
953  }
954 }
955 
956 
957 //====================================================================
958 void Fopr_Wilson::Fopr_Wilson_impl::mult_zm(Field& w, const Field& f)
959 {
960  int Ncol = CommonParameters::Nc();
961  int Nvc = 2 * Ncol;
962  int Ndf = 2 * Ncol * Ncol;
963  int ND = CommonParameters::Nd();
964  int Nx = CommonParameters::Nx();
965  int Ny = CommonParameters::Ny();
966  int Nz = CommonParameters::Nz();
967  int Nt = CommonParameters::Nt();
968  int Nst = Nx * Ny * Nz * Nt;
969 
970  double *v1;
971  double *v2;
972  double *u;
973 
974  v2 = w.ptr(0);
975  v1 = const_cast<Field *>(&f)->ptr(0);
976  u = const_cast<Field_G *>(m_U)->ptr(0);
977 
978  int idir = 2;
979 
980  int id1 = 0;
981  int id2 = Nvc;
982  int id3 = Nvc * 2;
983  int id4 = Nvc * 3;
984 
985  double vt1[Nvc], vt2[Nvc];
986  double wt1r, wt1i, wt2r, wt2i;
987  double bc2 = 1.0;
988  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
989 
990  int iz, nn;
991  double vcp1[Nvc * 2 * Nx * Ny * Nt], vcp2[Nvc * 2 * Nx * Ny * Nt];
992 
993  //- boundary part
994  iz = 0;
995  nn = Nz - 1;
996  for (int it = 0; it < Nt; it++) {
997  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
998  int in = Nvc * ND * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz);
999  int ig = Ndf * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz + idir * Nst);
1000  int ix1 = Nvc * 2 * (ixy + it * Nx * Ny);
1001  int ix2 = ix1 + Nvc;
1002 
1003  for (int ic = 0; ic < Ncol; ic++) {
1004  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + 1 + id3 + in];
1005  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + id3 + in];
1006  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + 1 + id4 + in];
1007  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + id4 + in];
1008  }
1009 
1010  for (int ic = 0; ic < Ncol; ic++) {
1011  int icr = 2 * ic;
1012  int ici = 2 * ic + 1;
1013  vcp1[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1);
1014  vcp1[ici + ix1] = mult_udagv_i(&u[icr + ig], vt1);
1015  vcp1[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2);
1016  vcp1[ici + ix2] = mult_udagv_i(&u[icr + ig], vt2);
1017  }
1018  }
1019  }
1020 
1021  Communicator::exchange(Nvc * 2 * Nx * Ny * Nt, vcp2, vcp1, 2, -1, 6);
1022 
1023  iz = 0;
1024  nn = Nz - 1;
1025  for (int it = 0; it < Nt; it++) {
1026  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
1027  int iv = Nvc * ND * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz);
1028  int ix1 = Nvc * 2 * (ixy + it * Nx * Ny);
1029  int ix2 = ix1 + Nvc;
1030 
1031  for (int ic = 0; ic < Ncol; ic++) {
1032  int icr = 2 * ic;
1033  int ici = 2 * ic + 1;
1034  v2[icr + id1 + iv] += bc2 * vcp2[icr + ix1];
1035  v2[ici + id1 + iv] += bc2 * vcp2[ici + ix1];
1036  v2[icr + id2 + iv] += bc2 * vcp2[icr + ix2];
1037  v2[ici + id2 + iv] += bc2 * vcp2[ici + ix2];
1038  v2[icr + id3 + iv] += -bc2 * vcp2[ici + ix1];
1039  v2[ici + id3 + iv] += bc2 * vcp2[icr + ix1];
1040  v2[icr + id4 + iv] += bc2 * vcp2[ici + ix2];
1041  v2[ici + id4 + iv] += -bc2 * vcp2[icr + ix2];
1042  }
1043  }
1044  }
1045 
1046  //- bulk part
1047  for (int it = 0; it < Nt; it++) {
1048  for (int iz = 1; iz < Nz; iz++) {
1049  int nn = iz - 1;
1050  for (int ixy = 0; ixy < (Nx * Ny); ixy++) {
1051  int iv = Nvc * ND * (ixy + iz * Nx * Ny + it * Nx * Ny * Nz);
1052  int in = Nvc * ND * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz);
1053  int ig = Ndf * (ixy + nn * Nx * Ny + it * Nx * Ny * Nz + idir * Nst);
1054 
1055  for (int ic = 0; ic < Ncol; ic++) {
1056  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + 1 + id3 + in];
1057  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + id3 + in];
1058  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + 1 + id4 + in];
1059  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + id4 + in];
1060  }
1061 
1062  for (int ic = 0; ic < Ncol; ic++) {
1063  int icr = 2 * ic;
1064  int ici = 2 * ic + 1;
1065  int ic1 = 0;
1066  int ic2 = Nvc;
1067  int ic3 = 2 * Nvc;
1068 
1069  wt1r = mult_udagv_r(&u[icr + ig], vt1);
1070  wt1i = mult_udagv_i(&u[icr + ig], vt1);
1071  wt2r = mult_udagv_r(&u[icr + ig], vt2);
1072  wt2i = mult_udagv_i(&u[icr + ig], vt2);
1073 
1074  v2[icr + id1 + iv] += wt1r;
1075  v2[ici + id1 + iv] += wt1i;
1076  v2[icr + id2 + iv] += wt2r;
1077  v2[ici + id2 + iv] += wt2i;
1078  v2[icr + id3 + iv] += -wt1i;
1079  v2[ici + id3 + iv] += wt1r;
1080  v2[icr + id4 + iv] += wt2i;
1081  v2[ici + id4 + iv] += -wt2r;
1082  }
1083  }
1084  }
1085  }
1086 }
1087 
1088 
1089 //====================================================================
1090 void Fopr_Wilson::Fopr_Wilson_impl::mult_tp_dirac(Field& w, const Field& f)
1091 {
1092  int Ncol = CommonParameters::Nc();
1093  int Nvc = 2 * Ncol;
1094  int Ndf = 2 * Ncol * Ncol;
1095  int ND = CommonParameters::Nd();
1096  int Nx = CommonParameters::Nx();
1097  int Ny = CommonParameters::Ny();
1098  int Nz = CommonParameters::Nz();
1099  int Nt = CommonParameters::Nt();
1100  int Nst = Nx * Ny * Nz * Nt;
1101 
1102  double *v1;
1103  double *v2;
1104  double *u;
1105 
1106  v2 = w.ptr(0);
1107  v1 = const_cast<Field *>(&f)->ptr(0);
1108  u = const_cast<Field_G *>(m_U)->ptr(0);
1109 
1110  int idir = 3;
1111 
1112  int id1 = 0;
1113  int id2 = Nvc;
1114  int id3 = Nvc * 2;
1115  int id4 = Nvc * 3;
1116 
1117  double vt1[Nvc], vt2[Nvc];
1118  double wt1r, wt1i, wt2r, wt2i;
1119  double bc2 = 1.0;
1120  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1121 
1122  int it, nn;
1123  double vcp1[Nvc * 2 * Nx * Ny * Nz], vcp2[Nvc * 2 * Nx * Ny * Nz];
1124 
1125  //- boundary part
1126  it = Nt - 1;
1127  nn = 0;
1128  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1129  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1130  int ix1 = Nvc * 2 * ixyz;
1131  int ix2 = ix1 + Nvc;
1132 
1133  for (int ic = 0; ic < Ncol; ic++) {
1134  vcp1[2 * ic + ix1] = 2.0 * bc2 * v1[2 * ic + id3 + in];
1135  vcp1[2 * ic + 1 + ix1] = 2.0 * bc2 * v1[2 * ic + 1 + id3 + in];
1136  vcp1[2 * ic + ix2] = 2.0 * bc2 * v1[2 * ic + id4 + in];
1137  vcp1[2 * ic + 1 + ix2] = 2.0 * bc2 * v1[2 * ic + 1 + id4 + in];
1138  }
1139  }
1140 
1141  Communicator::exchange(Nvc * 2 * Nx * Ny * Nz, vcp2, vcp1, 3, 1, 7);
1142 
1143  it = Nt - 1;
1144  nn = 0;
1145  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1146  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1147  int ig = Ndf * (ixyz + it * Nx * Ny * Nz + idir * Nst);
1148  int ix1 = Nvc * 2 * ixyz;
1149  int ix2 = ix1 + Nvc;
1150 
1151  for (int ic = 0; ic < Ncol; ic++) {
1152  int ic2 = ic * Nvc;
1153 
1154  wt1r = mult_uv_r(&u[ic2 + ig], &vcp2[ix1]);
1155  wt1i = mult_uv_i(&u[ic2 + ig], &vcp2[ix1]);
1156  wt2r = mult_uv_r(&u[ic2 + ig], &vcp2[ix2]);
1157  wt2i = mult_uv_i(&u[ic2 + ig], &vcp2[ix2]);
1158 
1159  v2[2 * ic + id3 + iv] += wt1r;
1160  v2[2 * ic + 1 + id3 + iv] += wt1i;
1161  v2[2 * ic + id4 + iv] += wt2r;
1162  v2[2 * ic + 1 + id4 + iv] += wt2i;
1163  }
1164  }
1165 
1166  //- bulk part
1167  for (int it = 0; it < (Nt - 1); it++) {
1168  int nn = it + 1;
1169  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1170  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1171  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1172  int ig = Ndf * (ixyz + it * Nx * Ny * Nz + idir * Nst);
1173 
1174  for (int ic = 0; ic < Ncol; ic++) {
1175  vt1[2 * ic] = 2.0 * v1[2 * ic + id3 + in];
1176  vt1[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id3 + in];
1177  vt2[2 * ic] = 2.0 * v1[2 * ic + id4 + in];
1178  vt2[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id4 + in];
1179  }
1180 
1181  for (int ic = 0; ic < Ncol; ic++) {
1182  int ic2 = ic * Nvc;
1183 
1184  wt1r = mult_uv_r(&u[ic2 + ig], vt1);
1185  wt1i = mult_uv_i(&u[ic2 + ig], vt1);
1186  wt2r = mult_uv_r(&u[ic2 + ig], vt2);
1187  wt2i = mult_uv_i(&u[ic2 + ig], vt2);
1188 
1189  v2[2 * ic + id3 + iv] += wt1r;
1190  v2[2 * ic + 1 + id3 + iv] += wt1i;
1191  v2[2 * ic + id4 + iv] += wt2r;
1192  v2[2 * ic + 1 + id4 + iv] += wt2i;
1193  }
1194  }
1195  }
1196 }
1197 
1198 
1199 //====================================================================
1200 void Fopr_Wilson::Fopr_Wilson_impl::mult_tp_chiral(Field& w, const Field& f)
1201 {
1202  int Ncol = CommonParameters::Nc();
1203  int Nvc = 2 * Ncol;
1204  int Ndf = 2 * Ncol * Ncol;
1205  int ND = CommonParameters::Nd();
1206  int Nx = CommonParameters::Nx();
1207  int Ny = CommonParameters::Ny();
1208  int Nz = CommonParameters::Nz();
1209  int Nt = CommonParameters::Nt();
1210  int Nst = Nx * Ny * Nz * Nt;
1211 
1212  double *v1;
1213  double *v2;
1214  double *u;
1215 
1216  v2 = w.ptr(0);
1217  v1 = const_cast<Field *>(&f)->ptr(0);
1218  u = const_cast<Field_G *>(m_U)->ptr(0);
1219 
1220  int idir = 3;
1221 
1222  int id1 = 0;
1223  int id2 = Nvc;
1224  int id3 = Nvc * 2;
1225  int id4 = Nvc * 3;
1226 
1227  double vt1[Nvc], vt2[Nvc];
1228  double wt1r, wt1i, wt2r, wt2i;
1229  double bc2 = 1.0;
1230  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1231 
1232  int it, nn;
1233  double vcp1[Nvc * 2 * Nx * Ny * Nz], vcp2[Nvc * 2 * Nx * Ny * Nz];
1234 
1235  //- boundary part
1236  it = Nt - 1;
1237  nn = 0;
1238  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1239  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1240  int ix1 = Nvc * 2 * ixyz;
1241  int ix2 = ix1 + Nvc;
1242 
1243  for (int ic = 0; ic < Ncol; ic++) {
1244  vcp1[2 * ic + ix1] = bc2 * (v1[2 * ic + id1 + in] + v1[2 * ic + id3 + in]);
1245  vcp1[2 * ic + 1 + ix1] = bc2 * (v1[2 * ic + 1 + id1 + in] + v1[2 * ic + 1 + id3 + in]);
1246  vcp1[2 * ic + ix2] = bc2 * (v1[2 * ic + id2 + in] + v1[2 * ic + id4 + in]);
1247  vcp1[2 * ic + 1 + ix2] = bc2 * (v1[2 * ic + 1 + id2 + in] + v1[2 * ic + 1 + id4 + in]);
1248  }
1249  }
1250 
1251  Communicator::exchange(Nvc * 2 * Nx * Ny * Nz, vcp2, vcp1, 3, 1, 7);
1252 
1253  it = Nt - 1;
1254  nn = 0;
1255  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1256  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1257  int ig = Ndf * (ixyz + it * Nx * Ny * Nz + idir * Nst);
1258  int ix1 = Nvc * 2 * ixyz;
1259  int ix2 = ix1 + Nvc;
1260 
1261  for (int ic = 0; ic < Ncol; ic++) {
1262  int ic2 = ic * Nvc;
1263 
1264  wt1r = mult_uv_r(&u[ic2 + ig], &vcp2[ix1]);
1265  wt1i = mult_uv_i(&u[ic2 + ig], &vcp2[ix1]);
1266  wt2r = mult_uv_r(&u[ic2 + ig], &vcp2[ix2]);
1267  wt2i = mult_uv_i(&u[ic2 + ig], &vcp2[ix2]);
1268 
1269  v2[2 * ic + id1 + iv] += wt1r;
1270  v2[2 * ic + 1 + id1 + iv] += wt1i;
1271  v2[2 * ic + id2 + iv] += wt2r;
1272  v2[2 * ic + 1 + id2 + iv] += wt2i;
1273  v2[2 * ic + id3 + iv] += wt1r;
1274  v2[2 * ic + 1 + id3 + iv] += wt1i;
1275  v2[2 * ic + id4 + iv] += wt2r;
1276  v2[2 * ic + 1 + id4 + iv] += wt2i;
1277  }
1278  }
1279 
1280  //- bulk part
1281  for (int it = 0; it < (Nt - 1); it++) {
1282  int nn = it + 1;
1283  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1284  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1285  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1286  int ig = Ndf * (ixyz + it * Nx * Ny * Nz + idir * Nst);
1287 
1288  for (int ic = 0; ic < Ncol; ic++) {
1289  vt1[2 * ic] = v1[2 * ic + id1 + in] + v1[2 * ic + id3 + in];
1290  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] + v1[2 * ic + 1 + id3 + in];
1291  vt2[2 * ic] = v1[2 * ic + id2 + in] + v1[2 * ic + id4 + in];
1292  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] + v1[2 * ic + 1 + id4 + in];
1293  }
1294 
1295  for (int ic = 0; ic < Ncol; ic++) {
1296  int ic2 = ic * Nvc;
1297 
1298  wt1r = mult_uv_r(&u[ic2 + ig], vt1);
1299  wt1i = mult_uv_i(&u[ic2 + ig], vt1);
1300  wt2r = mult_uv_r(&u[ic2 + ig], vt2);
1301  wt2i = mult_uv_i(&u[ic2 + ig], vt2);
1302 
1303  v2[2 * ic + id1 + iv] += wt1r;
1304  v2[2 * ic + 1 + id1 + iv] += wt1i;
1305  v2[2 * ic + id2 + iv] += wt2r;
1306  v2[2 * ic + 1 + id2 + iv] += wt2i;
1307  v2[2 * ic + id3 + iv] += wt1r;
1308  v2[2 * ic + 1 + id3 + iv] += wt1i;
1309  v2[2 * ic + id4 + iv] += wt2r;
1310  v2[2 * ic + 1 + id4 + iv] += wt2i;
1311  }
1312  }
1313  }
1314 }
1315 
1316 
1317 //====================================================================
1318 void Fopr_Wilson::Fopr_Wilson_impl::mult_tm_dirac(Field& w, const Field& f)
1319 {
1320  int Ncol = CommonParameters::Nc();
1321  int Nvc = 2 * Ncol;
1322  int Ndf = 2 * Ncol * Ncol;
1323  int ND = CommonParameters::Nd();
1324  int Nx = CommonParameters::Nx();
1325  int Ny = CommonParameters::Ny();
1326  int Nz = CommonParameters::Nz();
1327  int Nt = CommonParameters::Nt();
1328  int Nst = Nx * Ny * Nz * Nt;
1329 
1330  double *v1;
1331  double *v2;
1332  double *u;
1333 
1334  v2 = w.ptr(0);
1335  v1 = const_cast<Field *>(&f)->ptr(0);
1336  u = const_cast<Field_G *>(m_U)->ptr(0);
1337 
1338  int idir = 3;
1339 
1340  int id1 = 0;
1341  int id2 = Nvc;
1342  int id3 = Nvc * 2;
1343  int id4 = Nvc * 3;
1344 
1345  double vt1[Nvc], vt2[Nvc];
1346  double wt1r, wt1i, wt2r, wt2i;
1347  double bc2 = 1.0;
1348  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1349 
1350  int it, nn;
1351  double vcp1[Nvc * 2 * Nx * Ny * Nz], vcp2[Nvc * 2 * Nx * Ny * Nz];
1352 
1353  //- boundary part
1354  it = 0;
1355  nn = Nt - 1;
1356  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1357  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1358  int ig = Ndf * (ixyz + nn * Nx * Ny * Nz + idir * Nst);
1359  int ix1 = Nvc * 2 * ixyz;
1360  int ix2 = ix1 + Nvc;
1361 
1362  for (int ic = 0; ic < Ncol; ic++) {
1363  vt1[2 * ic] = 2.0 * v1[2 * ic + id1 + in];
1364  vt1[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id1 + in];
1365  vt2[2 * ic] = 2.0 * v1[2 * ic + id2 + in];
1366  vt2[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id2 + in];
1367  }
1368 
1369  for (int ic = 0; ic < Ncol; ic++) {
1370  int icr = 2 * ic;
1371  int ici = 2 * ic + 1;
1372  vcp1[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1);
1373  vcp1[ici + ix1] = mult_udagv_i(&u[icr + ig], vt1);
1374  vcp1[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2);
1375  vcp1[ici + ix2] = mult_udagv_i(&u[icr + ig], vt2);
1376  }
1377  }
1378 
1379  Communicator::exchange(Nvc * 2 * Nx * Ny * Nz, vcp2, vcp1, 3, -1, 8);
1380 
1381  it = 0;
1382  nn = Nt - 1;
1383  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1384  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1385  int ix1 = Nvc * 2 * ixyz;
1386  int ix2 = ix1 + Nvc;
1387 
1388  for (int ic = 0; ic < Ncol; ic++) {
1389  int icr = 2 * ic;
1390  int ici = 2 * ic + 1;
1391  v2[icr + id1 + iv] += bc2 * vcp2[icr + ix1];
1392  v2[ici + id1 + iv] += bc2 * vcp2[ici + ix1];
1393  v2[icr + id2 + iv] += bc2 * vcp2[icr + ix2];
1394  v2[ici + id2 + iv] += bc2 * vcp2[ici + ix2];
1395  }
1396  }
1397 
1398  //- bulk part
1399  for (int it = 1; it < Nt; it++) {
1400  int nn = it - 1;
1401  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1402  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1403  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1404  int ig = Ndf * (ixyz + nn * Nx * Ny * Nz + idir * Nst);
1405 
1406  for (int ic = 0; ic < Ncol; ic++) {
1407  vt1[2 * ic] = 2.0 * v1[2 * ic + id1 + in];
1408  vt1[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id1 + in];
1409  vt2[2 * ic] = 2.0 * v1[2 * ic + id2 + in];
1410  vt2[2 * ic + 1] = 2.0 * v1[2 * ic + 1 + id2 + in];
1411  }
1412 
1413  for (int ic = 0; ic < Ncol; ic++) {
1414  int icr = 2 * ic;
1415  int ici = 2 * ic + 1;
1416  int ic1 = 0;
1417  int ic2 = Nvc;
1418  int ic3 = 2 * Nvc;
1419 
1420  wt1r = mult_udagv_r(&u[icr + ig], vt1);
1421  wt1i = mult_udagv_i(&u[icr + ig], vt1);
1422  wt2r = mult_udagv_r(&u[icr + ig], vt2);
1423  wt2i = mult_udagv_i(&u[icr + ig], vt2);
1424 
1425  v2[icr + id1 + iv] += wt1r;
1426  v2[ici + id1 + iv] += wt1i;
1427  v2[icr + id2 + iv] += wt2r;
1428  v2[ici + id2 + iv] += wt2i;
1429  }
1430  }
1431  }
1432 }
1433 
1434 
1435 //====================================================================
1436 void Fopr_Wilson::Fopr_Wilson_impl::mult_tm_chiral(Field& w, const Field& f)
1437 {
1438  int Ncol = CommonParameters::Nc();
1439  int Nvc = 2 * Ncol;
1440  int Ndf = 2 * Ncol * Ncol;
1441  int ND = CommonParameters::Nd();
1442  int Nx = CommonParameters::Nx();
1443  int Ny = CommonParameters::Ny();
1444  int Nz = CommonParameters::Nz();
1445  int Nt = CommonParameters::Nt();
1446  int Nst = Nx * Ny * Nz * Nt;
1447 
1448  double *v1;
1449  double *v2;
1450  double *u;
1451 
1452  v2 = w.ptr(0);
1453  v1 = const_cast<Field *>(&f)->ptr(0);
1454  u = const_cast<Field_G *>(m_U)->ptr(0);
1455 
1456  int idir = 3;
1457 
1458  int id1 = 0;
1459  int id2 = Nvc;
1460  int id3 = Nvc * 2;
1461  int id4 = Nvc * 3;
1462 
1463  double vt1[Nvc], vt2[Nvc];
1464  double wt1r, wt1i, wt2r, wt2i;
1465  double bc2 = 1.0;
1466  if (Communicator::ipe(idir) == 0) bc2 = m_boundary[idir];
1467 
1468  int it, nn;
1469  double vcp1[Nvc * 2 * Nx * Ny * Nz], vcp2[Nvc * 2 * Nx * Ny * Nz];
1470 
1471  //- boundary part
1472  it = 0;
1473  nn = Nt - 1;
1474  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1475  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1476  int ig = Ndf * (ixyz + nn * Nx * Ny * Nz + idir * Nst);
1477  int ix1 = Nvc * 2 * ixyz;
1478  int ix2 = ix1 + Nvc;
1479 
1480  for (int ic = 0; ic < Ncol; ic++) {
1481  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + id3 + in];
1482  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + 1 + id3 + in];
1483  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + id4 + in];
1484  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + 1 + id4 + in];
1485  }
1486 
1487  for (int ic = 0; ic < Ncol; ic++) {
1488  int icr = 2 * ic;
1489  int ici = 2 * ic + 1;
1490  vcp1[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1);
1491  vcp1[ici + ix1] = mult_udagv_i(&u[icr + ig], vt1);
1492  vcp1[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2);
1493  vcp1[ici + ix2] = mult_udagv_i(&u[icr + ig], vt2);
1494  }
1495  }
1496 
1497  Communicator::exchange(Nvc * 2 * Nx * Ny * Nz, vcp2, vcp1, 3, -1, 8);
1498 
1499  it = 0;
1500  nn = Nt - 1;
1501  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1502  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1503  int ix1 = Nvc * 2 * ixyz;
1504  int ix2 = ix1 + Nvc;
1505 
1506  for (int ic = 0; ic < Ncol; ic++) {
1507  int icr = 2 * ic;
1508  int ici = 2 * ic + 1;
1509  v2[icr + id1 + iv] += bc2 * vcp2[icr + ix1];
1510  v2[ici + id1 + iv] += bc2 * vcp2[ici + ix1];
1511  v2[icr + id2 + iv] += bc2 * vcp2[icr + ix2];
1512  v2[ici + id2 + iv] += bc2 * vcp2[ici + ix2];
1513  v2[icr + id3 + iv] -= bc2 * vcp2[icr + ix1];
1514  v2[ici + id3 + iv] -= bc2 * vcp2[ici + ix1];
1515  v2[icr + id4 + iv] -= bc2 * vcp2[icr + ix2];
1516  v2[ici + id4 + iv] -= bc2 * vcp2[ici + ix2];
1517  }
1518  }
1519 
1520  //- bulk part
1521  for (int it = 1; it < Nt; it++) {
1522  int nn = it - 1;
1523  for (int ixyz = 0; ixyz < (Nx * Ny * Nz); ixyz++) {
1524  int iv = Nvc * ND * (ixyz + it * Nx * Ny * Nz);
1525  int in = Nvc * ND * (ixyz + nn * Nx * Ny * Nz);
1526  int ig = Ndf * (ixyz + nn * Nx * Ny * Nz + idir * Nst);
1527 
1528  for (int ic = 0; ic < Ncol; ic++) {
1529  vt1[2 * ic] = v1[2 * ic + id1 + in] - v1[2 * ic + id3 + in];
1530  vt1[2 * ic + 1] = v1[2 * ic + 1 + id1 + in] - v1[2 * ic + 1 + id3 + in];
1531  vt2[2 * ic] = v1[2 * ic + id2 + in] - v1[2 * ic + id4 + in];
1532  vt2[2 * ic + 1] = v1[2 * ic + 1 + id2 + in] - v1[2 * ic + 1 + id4 + in];
1533  }
1534 
1535  for (int ic = 0; ic < Ncol; ic++) {
1536  int icr = 2 * ic;
1537  int ici = 2 * ic + 1;
1538  int ic1 = 0;
1539  int ic2 = Nvc;
1540  int ic3 = 2 * Nvc;
1541 
1542  wt1r = mult_udagv_r(&u[icr + ig], vt1);
1543  wt1i = mult_udagv_i(&u[icr + ig], vt1);
1544  wt2r = mult_udagv_r(&u[icr + ig], vt2);
1545  wt2i = mult_udagv_i(&u[icr + ig], vt2);
1546 
1547  v2[icr + id1 + iv] += wt1r;
1548  v2[ici + id1 + iv] += wt1i;
1549  v2[icr + id2 + iv] += wt2r;
1550  v2[ici + id2 + iv] += wt2i;
1551  v2[icr + id3 + iv] -= wt1r;
1552  v2[ici + id3 + iv] -= wt1i;
1553  v2[icr + id4 + iv] -= wt2r;
1554  v2[ici + id4 + iv] -= wt2i;
1555  }
1556  }
1557  }
1558 }
1559 
1560 
1561 //====================================================================
1562 //============================================================END=====