Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_Wilson_eo_impl_thread.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Wilson_eo_impl.h"
15 
16 #include "field_G.h"
17 #include "bridgeIO.h"
18 using Bridge::vout;
19 
20 #include "threadManager_OpenMP.h"
21 
22 
23 #if defined USE_GROUP_SU3
24 #include "fopr_Wilson_impl_SU3.inc"
25 #elif defined USE_GROUP_SU2
26 #include "fopr_Wilson_impl_SU2.inc"
27 #elif defined USE_GROUP_SU_N
28 #include "fopr_Wilson_impl_SU_N.inc"
29 #endif
30 
31 //====================================================================
33 {
35 
36  // The following setup corresponds to uniform division of volume.
37  if (m_Nthread <= m_Nt) {
39  } else if (m_Nthread <= m_Nz * m_Nt) {
40  m_Ntask_t = m_Nt;
41  } else {
42  vout.crucial(m_vl, "%s: Too large Nthread = %d\n", class_name.c_str(), m_Nthread);
43  abort();
44  }
45 
47 
48  if (m_Ntask_z * m_Ntask_t != m_Nthread) {
49  vout.crucial(m_vl, "%s: Nz = %d and Nt = %d do not match Nthread = %d\n",
50  class_name.c_str(), m_Nz, m_Nt, m_Nthread);
51  abort();
52  }
53 
55  m_Mz = m_Nz / m_Ntask_z;
56  m_Mt = m_Nt / m_Ntask_t;
57 
58  if (m_Mz * m_Ntask_z != m_Nz) {
59  vout.crucial(m_vl, "%s: Mz = %d and Ntask_z = %d do not match Nz = %d\n",
60  class_name.c_str(), m_Mz, m_Ntask_z, m_Nz);
61  abort();
62  }
63 
64  if (m_Mt * m_Ntask_t != m_Nt) {
65  vout.crucial(m_vl, "%s: Mt = %d and Ntask_t = %d do not match Nt = %d\n",
66  class_name.c_str(), m_Mt, m_Ntask_t, m_Nt);
67  abort();
68  }
69 
70  vout.general(m_vl, " Nthread = %d\n", m_Nthread);
71  vout.general(m_vl, " Ntask = %d\n", m_Ntask);
72  vout.general(m_vl, " Ntask_z = %d Ntask_t = %d\n", m_Ntask_z, m_Ntask_t);
73  vout.general(m_vl, " Mz = %d Mt = %d\n", m_Mz, m_Mt);
74 
75  // setup of arguments
76  int Nxy2 = m_Nx2 * m_Ny;
77  m_arg.resize(m_Ntask);
78  for (int ith_t = 0; ith_t < m_Ntask_t; ++ith_t) {
79  for (int ith_z = 0; ith_z < m_Ntask_z; ++ith_z) {
80  int itask = ith_z + m_Ntask_z * ith_t;
81 
82  m_arg[itask].isite = (ith_z * m_Mz + ith_t * (m_Nz * m_Mt)) * Nxy2;
83 
84  m_arg[itask].kt0 = 0;
85  m_arg[itask].kt1 = 0;
86  m_arg[itask].kz0 = 0;
87  m_arg[itask].kz1 = 0;
88  if (ith_t == 0) m_arg[itask].kt0 = 1;
89  if (ith_z == 0) m_arg[itask].kz0 = 1;
90  if (ith_t == m_Ntask_t - 1) m_arg[itask].kt1 = 1;
91  if (ith_z == m_Ntask_z - 1) m_arg[itask].kz1 = 1;
92 
93  m_arg[itask].isite_cpx = itask * m_Mz * m_Mt * (m_Ny / 2);
94  m_arg[itask].isite_cpy = itask * m_Mz * m_Mt * m_Nx2;
95  m_arg[itask].isite_cpz = ith_t * m_Mt * Nxy2;
96  m_arg[itask].isite_cpt = ith_z * m_Mz * Nxy2;
97  }
98  }
99 
100 }
101 
102 //====================================================================
104  double *w, double fac)
105 {
106  int Nvcd = m_Nvc * m_Nd;
107  int Nvxy = Nvcd * m_Nx2 * m_Ny;
108 
109  int isite = m_arg[itask].isite;
110  double *wp = &w[Nvcd * isite];
111 
112  for (int it = 0; it < m_Mt; ++it) {
113  for (int iz = 0; iz < m_Mz; ++iz) {
114  for (int ivxy = 0; ivxy < Nvxy; ++ivxy) {
115  int iv = ivxy + Nvxy * (iz + m_Nz * it);
116  wp[iv] = fac * wp[iv];
117  }
118  }
119  }
120 
121 }
122 
123 //====================================================================
125  double *v)
126 {
127  int Nvcd = m_Nvc * m_Nd;
128  int Nvxy = Nvcd * m_Nx2 * m_Ny;
129 
130  int isite = m_arg[itask].isite;
131  double *wp = &v[Nvcd * isite];
132 
133  for (int it = 0; it < m_Mt; ++it) {
134  for (int iz = 0; iz < m_Mz; ++iz) {
135  for (int ivxy = 0; ivxy < Nvxy; ++ivxy) {
136  int iv = ivxy + Nvxy * (iz + m_Nz * it);
137  wp[iv] = 0.0;
138  }
139  }
140  }
141 
142 }
143 
144 //====================================================================
146  double *vcp1, double *v1, int ieo)
147 {
148  int Nvc2 = 2 * m_Nvc;
149  int Nvcd = m_Nvc * m_Nd;
150  int Nvcd2 = Nvcd / 2;
151 
152  int id1 = 0;
153  int id2 = m_Nvc;
154  int id3 = m_Nvc * 2;
155  int id4 = m_Nvc * 3;
156 
157  int isite = m_arg[itask].isite;
158  int isite_cp = m_arg[itask].isite_cpx;
159  int iyzt0 = isite / m_Nx2;
160 
161  double *w2 = &vcp1[Nvcd2 * isite_cp];
162  double *w1 = &v1[Nvcd * isite];
163 
164  int idir = 0;
165  double bc2 = m_boundary2[idir];
166 
167  int ix = 0;
168  int ibf = 0;
169 
170  for (int it = 0; it < m_Mt; ++it) {
171  for (int iz = 0; iz < m_Mz; ++iz) {
172  for (int iy = 0; iy < m_Ny; ++iy) {
173  int iyzt = iy + m_Ny * (iz + m_Nz * it);
174  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
175  if (Leo == 1) {
176  int is = ix + m_Nx2 * iyzt;
177  int in = Nvcd * is;
178 
179  int ix1 = Nvc2 * ibf;
180  int ix2 = ix1 + m_Nvc;
181 
182  for (int ic = 0; ic < m_Nc; ++ic) {
183  w2[2 * ic + ix1] = bc2 * (w1[2 * ic + id1 + in] - w1[2 * ic + 1 + id4 + in]);
184  w2[2 * ic + 1 + ix1] = bc2 * (w1[2 * ic + 1 + id1 + in] + w1[2 * ic + id4 + in]);
185  w2[2 * ic + ix2] = bc2 * (w1[2 * ic + id2 + in] - w1[2 * ic + 1 + id3 + in]);
186  w2[2 * ic + 1 + ix2] = bc2 * (w1[2 * ic + 1 + id2 + in] + w1[2 * ic + id3 + in]);
187  }
188  ++ibf;
189  }
190  }
191  }
192  }
193 
194 }
195 
196 //====================================================================
198  double *v2, double *vcp2, int ieo)
199 {
200  int Nvc2 = 2 * m_Nvc;
201  int Nvcd = m_Nvc * m_Nd;
202  int Nvcd2 = Nvcd / 2;
203 
204  int id1 = 0;
205  int id2 = m_Nvc;
206  int id3 = m_Nvc * 2;
207  int id4 = m_Nvc * 3;
208 
209  int idir = 0;
210 
211  double wt1r, wt1i, wt2r, wt2i;
212 
213  int isite = m_arg[itask].isite;
214  int isite_cp = m_arg[itask].isite_cpx;
215  int iyzt0 = isite / m_Nx2;
216 
217  double *w2 = &v2[Nvcd * isite];
218  double *w1 = &vcp2[Nvcd2 * isite_cp];
219  double *u = const_cast<Field_G *>(m_U)->ptr(
220  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
221 
222  int ix = m_Nx2 - 1;
223  int ibf = 0;
224 
225  for (int it = 0; it < m_Mt; ++it) {
226  for (int iz = 0; iz < m_Mz; ++iz) {
227  for (int iy = 0; iy < m_Ny; ++iy) {
228  int iyzt = iy + m_Ny * (iz + m_Nz * it);
229  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
230 
231  if (Leo == 1) {
232  int is = ix + m_Nx2 * iyzt;
233  int iv = Nvcd * is;
234  int ig = m_Ndf * is;
235  int ix1 = Nvc2 * ibf;
236  int ix2 = ix1 + m_Nvc;
237 
238  for (int ic = 0; ic < m_Nc; ++ic) {
239  int ic2 = ic * m_Nvc;
240  wt1r = mult_uv_r(&u[ic2 + ig], &w1[ix1], m_Nc);
241  wt1i = mult_uv_i(&u[ic2 + ig], &w1[ix1], m_Nc);
242  wt2r = mult_uv_r(&u[ic2 + ig], &w1[ix2], m_Nc);
243  wt2i = mult_uv_i(&u[ic2 + ig], &w1[ix2], m_Nc);
244  w2[2 * ic + id1 + iv] += wt1r;
245  w2[2 * ic + 1 + id1 + iv] += wt1i;
246  w2[2 * ic + id2 + iv] += wt2r;
247  w2[2 * ic + 1 + id2 + iv] += wt2i;
248  w2[2 * ic + id3 + iv] += wt2i;
249  w2[2 * ic + 1 + id3 + iv] += -wt2r;
250  w2[2 * ic + id4 + iv] += wt1i;
251  w2[2 * ic + 1 + id4 + iv] += -wt1r;
252  }
253  ++ibf;
254  }
255  }
256  }
257  }
258 
259 }
260 
261 //====================================================================
263  double *v2, double *v1, int ieo)
264 {
265  int Nvcd = m_Nvc * m_Nd;
266 
267  int id1 = 0;
268  int id2 = m_Nvc;
269  int id3 = m_Nvc * 2;
270  int id4 = m_Nvc * 3;
271 
272  int idir = 0;
273 
274  double vt1[m_Nvc], vt2[m_Nvc];
275  double wt1r, wt1i, wt2r, wt2i;
276 
277  int isite = m_arg[itask].isite;
278  int iyzt0 = isite / m_Nx2;
279 
280  double *w2 = &v2[Nvcd * isite];
281  double *w1 = &v1[Nvcd * isite];
282  double *u = const_cast<Field_G *>(m_U)->ptr(
283  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
284 
285  for (int it = 0; it < m_Mt; ++it) {
286  for (int iz = 0; iz < m_Mz; ++iz) {
287  for (int iy = 0; iy < m_Ny; ++iy) {
288  int iyzt = iy + m_Ny * (iz + m_Nz * it);
289  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
290  for (int ix = 0; ix < m_Nx2 - Leo; ++ix) {
291  int is = ix + m_Nx2 * iyzt;
292  int iv = Nvcd * is;
293  int in = Nvcd * (is + Leo);
294  int ig = m_Ndf * is;
295 
296  for (int ic = 0; ic < m_Nc; ++ic) {
297  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + 1 + id4 + in];
298  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] + w1[2 * ic + id4 + in];
299  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + 1 + id3 + in];
300  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + id3 + in];
301  }
302 
303  for (int ic = 0; ic < m_Nc; ++ic) {
304  int ic2 = ic * m_Nvc;
305 
306  wt1r = mult_uv_r(&u[ic2 + ig], vt1, m_Nc);
307  wt1i = mult_uv_i(&u[ic2 + ig], vt1, m_Nc);
308  wt2r = mult_uv_r(&u[ic2 + ig], vt2, m_Nc);
309  wt2i = mult_uv_i(&u[ic2 + ig], vt2, m_Nc);
310 
311  w2[2 * ic + id1 + iv] += wt1r;
312  w2[2 * ic + 1 + id1 + iv] += wt1i;
313  w2[2 * ic + id2 + iv] += wt2r;
314  w2[2 * ic + 1 + id2 + iv] += wt2i;
315  w2[2 * ic + id3 + iv] += wt2i;
316  w2[2 * ic + 1 + id3 + iv] += -wt2r;
317  w2[2 * ic + id4 + iv] += wt1i;
318  w2[2 * ic + 1 + id4 + iv] += -wt1r;
319  }
320  }
321  }
322  }
323  }
324 
325 }
326 
327 //====================================================================
329  double *vcp1, double *v1, int ieo)
330 {
331  int Nvc2 = 2 * m_Nvc;
332  int Nvcd = m_Nvc * m_Nd;
333  int Nvcd2 = Nvcd / 2;
334 
335  int id1 = 0;
336  int id2 = m_Nvc;
337  int id3 = m_Nvc * 2;
338  int id4 = m_Nvc * 3;
339 
340  int idir = 0;
341 
342  int isite = m_arg[itask].isite;
343  int isite_cp = m_arg[itask].isite_cpx;
344  int iyzt0 = isite / m_Nx2;
345 
346  double *w2 = &vcp1[Nvcd2 * isite_cp];
347  double *w1 = &v1[Nvcd * isite];
348  double *u = const_cast<Field_G *>(m_U)->ptr(
349  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
350 
351  double vt1[m_Nvc], vt2[m_Nvc];
352 
353  int ix = m_Nx2 - 1;
354  int ibf = 0;
355 
356  for (int it = 0; it < m_Mt; ++it) {
357  for (int iz = 0; iz < m_Mz; ++iz) {
358  for (int iy = 0; iy < m_Ny; ++iy) {
359  int iyzt = iy + m_Ny * (iz + m_Nz * it);
360  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
361  if (Leo == 0) {
362  int is = ix + m_Nx2 * iyzt;
363  int in = Nvcd * is;
364  int ig = m_Ndf * is;
365 
366  int ix1 = Nvc2 * ibf;
367  int ix2 = ix1 + m_Nvc;
368 
369  for (int ic = 0; ic < m_Nc; ++ic) {
370  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + 1 + id4 + in];
371  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + id4 + in];
372  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + 1 + id3 + in];
373  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + id3 + in];
374  }
375 
376  for (int ic = 0; ic < m_Nc; ++ic) {
377  int icr = 2 * ic;
378  w2[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1, m_Nc);
379  w2[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1, m_Nc);
380  w2[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2, m_Nc);
381  w2[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2, m_Nc);
382  }
383  ++ibf;
384  }
385  }
386  }
387  }
388 
389 }
390 
391 //====================================================================
393  double *v2, double *vcp2, int ieo)
394 {
395  int Nvc2 = 2 * m_Nvc;
396  int Nvcd = m_Nvc * m_Nd;
397  int Nvcd2 = Nvcd / 2;
398 
399  int id1 = 0;
400  int id2 = m_Nvc;
401  int id3 = m_Nvc * 2;
402  int id4 = m_Nvc * 3;
403 
404  int idir = 0;
405  double bc2 = m_boundary2[idir];
406 
407  // double wt1r, wt1i, wt2r, wt2i;
408 
409  int isite = m_arg[itask].isite;
410  int isite_cp = m_arg[itask].isite_cpx;
411  int iyzt0 = isite / m_Nx2;
412 
413  double *w2 = &v2[Nvcd * isite];
414  double *w1 = &vcp2[Nvcd2 * isite_cp];
415 
416  int ix = 0;
417  int ibf = 0;
418 
419  for (int it = 0; it < m_Mt; ++it) {
420  for (int iz = 0; iz < m_Mz; ++iz) {
421  for (int iy = 0; iy < m_Ny; ++iy) {
422  int iyzt = iy + m_Ny * (iz + m_Nz * it);
423  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
424  if (Leo == 0) {
425  int is = ix + m_Nx2 * iyzt;
426  int iv = Nvcd * is;
427 
428  int ix1 = Nvc2 * ibf;
429  int ix2 = ix1 + m_Nvc;
430 
431  for (int ic = 0; ic < m_Nc; ++ic) {
432  int icr = 2 * ic;
433  int ici = 2 * ic + 1;
434  w2[icr + id1 + iv] += bc2 * w1[icr + ix1];
435  w2[ici + id1 + iv] += bc2 * w1[ici + ix1];
436  w2[icr + id2 + iv] += bc2 * w1[icr + ix2];
437  w2[ici + id2 + iv] += bc2 * w1[ici + ix2];
438  w2[icr + id3 + iv] += -bc2 * w1[ici + ix2];
439  w2[ici + id3 + iv] += +bc2 * w1[icr + ix2];
440  w2[icr + id4 + iv] += -bc2 * w1[ici + ix1];
441  w2[ici + id4 + iv] += +bc2 * w1[icr + ix1];
442  }
443  ++ibf;
444  }
445  }
446  }
447  }
448 
449 }
450 
451 //====================================================================
453  double *v2, double *v1, int ieo)
454 {
455  int Nvcd = m_Nvc * m_Nd;
456 
457  int id1 = 0;
458  int id2 = m_Nvc;
459  int id3 = m_Nvc * 2;
460  int id4 = m_Nvc * 3;
461 
462  int idir = 0;
463 
464  double vt1[m_Nvc], vt2[m_Nvc];
465  double wt1r, wt1i, wt2r, wt2i;
466 
467  int isite = m_arg[itask].isite;
468  int iyzt0 = isite / m_Nx2;
469 
470  double *w2 = &v2[Nvcd * isite];
471  double *w1 = &v1[Nvcd * isite];
472  double *u = const_cast<Field_G *>(m_U)->ptr(
473  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
474 
475  for (int it = 0; it < m_Mt; ++it) {
476  for (int iz = 0; iz < m_Mz; ++iz) {
477  for (int iy = 0; iy < m_Ny; ++iy) {
478  int iyzt = iy + m_Ny * (iz + m_Nz * it);
479  int Leo = ieo + (1 - 2 * ieo) * m_Leo[iyzt0 + iyzt];
480  int Meo = 1 - Leo;
481  for (int ix = Meo; ix < m_Nx2; ++ix) {
482  int is = ix + m_Nx2 * iyzt;
483  int iv = Nvcd * is;
484  int in = Nvcd * (is - Meo);
485  int ig = m_Ndf * (is - Meo);
486 
487  for (int ic = 0; ic < m_Nc; ++ic) {
488  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + 1 + id4 + in];
489  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + id4 + in];
490  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + 1 + id3 + in];
491  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + id3 + in];
492  }
493 
494  for (int ic = 0; ic < m_Nc; ++ic) {
495  int ic2 = 2 * ic;
496 
497  wt1r = mult_udagv_r(&u[ic2 + ig], vt1, m_Nc);
498  wt1i = mult_udagv_i(&u[ic2 + ig], vt1, m_Nc);
499  wt2r = mult_udagv_r(&u[ic2 + ig], vt2, m_Nc);
500  wt2i = mult_udagv_i(&u[ic2 + ig], vt2, m_Nc);
501 
502  w2[2 * ic + id1 + iv] += wt1r;
503  w2[2 * ic + 1 + id1 + iv] += wt1i;
504  w2[2 * ic + id2 + iv] += wt2r;
505  w2[2 * ic + 1 + id2 + iv] += wt2i;
506  w2[2 * ic + id3 + iv] += -wt2i;
507  w2[2 * ic + 1 + id3 + iv] += +wt2r;
508  w2[2 * ic + id4 + iv] += -wt1i;
509  w2[2 * ic + 1 + id4 + iv] += +wt1r;
510  }
511  }
512  }
513  }
514  }
515 
516 }
517 
518 //====================================================================
520  double *vcp1, double *v1, int ieo)
521 {
522  int Nvc2 = 2 * m_Nvc;
523  int Nvcd = m_Nvc * m_Nd;
524  int Nvcd2 = Nvcd / 2;
525 
526  int id1 = 0;
527  int id2 = m_Nvc;
528  int id3 = m_Nvc * 2;
529  int id4 = m_Nvc * 3;
530 
531  int isite = m_arg[itask].isite;
532  int isite_cp = m_arg[itask].isite_cpy;
533 
534  double *w2 = &vcp1[Nvcd2 * isite_cp];
535  double *w1 = &v1[Nvcd * isite];
536 
537  int idir = 1;
538  double bc2 = m_boundary2[idir];
539 
540  int iy = 0;
541 
542  for (int it = 0; it < m_Mt; ++it) {
543  for (int iz = 0; iz < m_Mz; ++iz) {
544  for (int ix = 0; ix < m_Nx2; ++ix) {
545  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
546  int is2 = ix + m_Nx2 * (iz + m_Mz * it);
547  int in = Nvcd * is;
548  int ix1 = Nvc2 * is2;
549  int ix2 = ix1 + m_Nvc;
550 
551  for (int ic = 0; ic < m_Nc; ++ic) {
552  w2[2 * ic + ix1] = bc2 * (w1[2 * ic + id1 + in] + w1[2 * ic + id4 + in]);
553  w2[2 * ic + 1 + ix1] = bc2 * (w1[2 * ic + 1 + id1 + in] + w1[2 * ic + 1 + id4 + in]);
554  w2[2 * ic + ix2] = bc2 * (w1[2 * ic + id2 + in] - w1[2 * ic + id3 + in]);
555  w2[2 * ic + 1 + ix2] = bc2 * (w1[2 * ic + 1 + id2 + in] - w1[2 * ic + 1 + id3 + in]);
556  }
557  }
558  }
559  }
560 
561 }
562 
563 //====================================================================
565  double *v2, double *vcp2, int ieo)
566 {
567  int Nvc2 = 2 * m_Nvc;
568  int Nvcd = m_Nvc * m_Nd;
569  int Nvcd2 = Nvcd / 2;
570 
571  int id1 = 0;
572  int id2 = m_Nvc;
573  int id3 = m_Nvc * 2;
574  int id4 = m_Nvc * 3;
575 
576  int idir = 1;
577 
578  double wt1r, wt1i, wt2r, wt2i;
579 
580  int isite = m_arg[itask].isite;
581  int isite_cp = m_arg[itask].isite_cpy;
582 
583  double *w2 = &v2[Nvcd * isite];
584  double *w1 = &vcp2[Nvcd2 * isite_cp];
585  double *u = const_cast<Field_G *>(m_U)->ptr(
586  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
587 
588  int iy = m_Ny - 1;
589 
590  for (int it = 0; it < m_Mt; ++it) {
591  for (int iz = 0; iz < m_Mz; ++iz) {
592  for (int ix = 0; ix < m_Nx2; ++ix) {
593  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
594  int is2 = ix + m_Nx2 * (iz + m_Mz * it);
595  int iv = Nvcd * is;
596  int ig = m_Ndf * is;
597  int ix1 = Nvc2 * is2;
598  int ix2 = ix1 + m_Nvc;
599 
600  for (int ic = 0; ic < m_Nc; ++ic) {
601  int ic2 = ic * m_Nvc;
602 
603  wt1r = mult_uv_r(&u[ic2 + ig], &w1[ix1], m_Nc);
604  wt1i = mult_uv_i(&u[ic2 + ig], &w1[ix1], m_Nc);
605  wt2r = mult_uv_r(&u[ic2 + ig], &w1[ix2], m_Nc);
606  wt2i = mult_uv_i(&u[ic2 + ig], &w1[ix2], m_Nc);
607 
608  w2[2 * ic + id1 + iv] += wt1r;
609  w2[2 * ic + 1 + id1 + iv] += wt1i;
610  w2[2 * ic + id2 + iv] += wt2r;
611  w2[2 * ic + 1 + id2 + iv] += wt2i;
612  w2[2 * ic + id3 + iv] += -wt2r;
613  w2[2 * ic + 1 + id3 + iv] += -wt2i;
614  w2[2 * ic + id4 + iv] += wt1r;
615  w2[2 * ic + 1 + id4 + iv] += wt1i;
616  }
617  }
618  }
619  }
620 
621 }
622 
623 //====================================================================
625  double *v2, double *v1, int ieo)
626 {
627  int Nvcd = m_Nvc * m_Nd;
628 
629  int id1 = 0;
630  int id2 = m_Nvc;
631  int id3 = m_Nvc * 2;
632  int id4 = m_Nvc * 3;
633 
634  int idir = 1;
635 
636  double vt1[m_Nvc], vt2[m_Nvc];
637  double wt1r, wt1i, wt2r, wt2i;
638 
639  int isite = m_arg[itask].isite;
640 
641  double *w2 = &v2[Nvcd * isite];
642  double *w1 = &v1[Nvcd * isite];
643  double *u = const_cast<Field_G *>(m_U)->ptr(
644  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
645 
646  for (int it = 0; it < m_Mt; ++it) {
647  for (int iz = 0; iz < m_Mz; ++iz) {
648  for (int iy = 0; iy < m_Ny - 1; ++iy) {
649  for (int ix = 0; ix < m_Nx2; ++ix) {
650  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
651  int iv = Nvcd * is;
652  int in = Nvcd * (is + m_Nx2);
653  int ig = m_Ndf * is;
654 
655  for (int ic = 0; ic < m_Nc; ++ic) {
656  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + id4 + in];
657  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] + w1[2 * ic + 1 + id4 + in];
658  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + id3 + in];
659  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + 1 + id3 + in];
660  }
661 
662  for (int ic = 0; ic < m_Nc; ++ic) {
663  int ic2 = ic * m_Nvc;
664 
665  wt1r = mult_uv_r(&u[ic2 + ig], vt1, m_Nc);
666  wt1i = mult_uv_i(&u[ic2 + ig], vt1, m_Nc);
667  wt2r = mult_uv_r(&u[ic2 + ig], vt2, m_Nc);
668  wt2i = mult_uv_i(&u[ic2 + ig], vt2, m_Nc);
669 
670  w2[2 * ic + id1 + iv] += wt1r;
671  w2[2 * ic + 1 + id1 + iv] += wt1i;
672  w2[2 * ic + id2 + iv] += wt2r;
673  w2[2 * ic + 1 + id2 + iv] += wt2i;
674  w2[2 * ic + id3 + iv] += -wt2r;
675  w2[2 * ic + 1 + id3 + iv] += -wt2i;
676  w2[2 * ic + id4 + iv] += wt1r;
677  w2[2 * ic + 1 + id4 + iv] += wt1i;
678  }
679  }
680  }
681  }
682  }
683 
684 }
685 
686 //====================================================================
688  double *vcp1, double *v1, int ieo)
689 {
690  int Nvc2 = 2 * m_Nvc;
691  int Nvcd = m_Nvc * m_Nd;
692  int Nvcd2 = Nvcd / 2;
693 
694  int id1 = 0;
695  int id2 = m_Nvc;
696  int id3 = m_Nvc * 2;
697  int id4 = m_Nvc * 3;
698 
699  int idir = 1;
700 
701  int isite = m_arg[itask].isite;
702  int isite_cp = m_arg[itask].isite_cpy;
703 
704  double *w2 = &vcp1[Nvcd2 * isite_cp];
705  double *w1 = &v1[Nvcd * isite];
706  double *u = const_cast<Field_G *>(m_U)->ptr(
707  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
708 
709  double vt1[m_Nvc], vt2[m_Nvc];
710 
711  int iy = m_Ny - 1;
712 
713  for (int it = 0; it < m_Mt; ++it) {
714  for (int iz = 0; iz < m_Mz; ++iz) {
715  for (int ix = 0; ix < m_Nx2; ++ix) {
716  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
717  int is2 = ix + m_Nx2 * (iz + m_Mz * it);
718  int in = Nvcd * is;
719  int ig = m_Ndf * is;
720  int ix1 = Nvc2 * is2;
721  int ix2 = ix1 + m_Nvc;
722 
723  for (int ic = 0; ic < m_Nc; ++ic) {
724  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + id4 + in];
725  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + 1 + id4 + in];
726  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + id3 + in];
727  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + 1 + id3 + in];
728  }
729 
730  for (int ic = 0; ic < m_Nc; ++ic) {
731  int icr = 2 * ic;
732  w2[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1, m_Nc);
733  w2[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1, m_Nc);
734  w2[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2, m_Nc);
735  w2[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2, m_Nc);
736  }
737  }
738  }
739  }
740 
741 }
742 
743 //====================================================================
745  double *v2, double *vcp2, int ieo)
746 {
747  int Nvc2 = 2 * m_Nvc;
748  int Nvcd = m_Nvc * m_Nd;
749  int Nvcd2 = Nvcd / 2;
750 
751  int id1 = 0;
752  int id2 = m_Nvc;
753  int id3 = m_Nvc * 2;
754  int id4 = m_Nvc * 3;
755 
756  int idir = 1;
757  double bc2 = m_boundary2[idir];
758 
759  // double wt1r, wt1i, wt2r, wt2i;
760 
761  int isite = m_arg[itask].isite;
762  int isite_cp = m_arg[itask].isite_cpy;
763 
764  double *w2 = &v2[Nvcd * isite];
765  double *w1 = &vcp2[Nvcd2 * isite_cp];
766 
767  int iy = 0;
768 
769  for (int it = 0; it < m_Mt; ++it) {
770  for (int iz = 0; iz < m_Mz; ++iz) {
771  for (int ix = 0; ix < m_Nx2; ++ix) {
772  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
773  int is2 = ix + m_Nx2 * (iz + m_Mz * it);
774  int iv = Nvcd * is;
775  int ix1 = Nvc2 * is2;
776  int ix2 = ix1 + m_Nvc;
777 
778  for (int ic = 0; ic < m_Nc; ++ic) {
779  int icr = 2 * ic;
780  int ici = 2 * ic + 1;
781  w2[icr + id1 + iv] += bc2 * w1[icr + ix1];
782  w2[ici + id1 + iv] += bc2 * w1[ici + ix1];
783  w2[icr + id2 + iv] += bc2 * w1[icr + ix2];
784  w2[ici + id2 + iv] += bc2 * w1[ici + ix2];
785  w2[icr + id3 + iv] += bc2 * w1[icr + ix2];
786  w2[ici + id3 + iv] += bc2 * w1[ici + ix2];
787  w2[icr + id4 + iv] += -bc2 * w1[icr + ix1];
788  w2[ici + id4 + iv] += -bc2 * w1[ici + ix1];
789  }
790  }
791  }
792  }
793 
794 }
795 
796 //====================================================================
798  double *v2, double *v1, int ieo)
799 {
800  int Nvcd = m_Nvc * m_Nd;
801 
802  int id1 = 0;
803  int id2 = m_Nvc;
804  int id3 = m_Nvc * 2;
805  int id4 = m_Nvc * 3;
806 
807  int idir = 1;
808 
809  double vt1[m_Nvc], vt2[m_Nvc];
810  double wt1r, wt1i, wt2r, wt2i;
811 
812  int isite = m_arg[itask].isite;
813 
814  double *w2 = &v2[Nvcd * isite];
815  double *w1 = &v1[Nvcd * isite];
816  double *u = const_cast<Field_G *>(m_U)->ptr(
817  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
818 
819  for (int it = 0; it < m_Mt; ++it) {
820  for (int iz = 0; iz < m_Mz; ++iz) {
821  for (int iy = 1; iy < m_Ny; ++iy) {
822  for (int ix = 0; ix < m_Nx2; ++ix) {
823  int is = ix + m_Nx2 * (iy + m_Ny * (iz + m_Nz * it));
824  int iv = Nvcd * is;
825  int in = Nvcd * (is - m_Nx2);
826  int ig = m_Ndf * (is - m_Nx2);
827 
828  for (int ic = 0; ic < m_Nc; ++ic) {
829  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + id4 + in];
830  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + 1 + id4 + in];
831  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + id3 + in];
832  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + 1 + id3 + in];
833  }
834 
835  for (int ic = 0; ic < m_Nc; ++ic) {
836  int ic2 = 2 * ic;
837  wt1r = mult_udagv_r(&u[ic2 + ig], vt1, m_Nc);
838  wt1i = mult_udagv_i(&u[ic2 + ig], vt1, m_Nc);
839  wt2r = mult_udagv_r(&u[ic2 + ig], vt2, m_Nc);
840  wt2i = mult_udagv_i(&u[ic2 + ig], vt2, m_Nc);
841 
842  w2[ic2 + id1 + iv] += wt1r;
843  w2[ic2 + 1 + id1 + iv] += wt1i;
844  w2[ic2 + id2 + iv] += wt2r;
845  w2[ic2 + 1 + id2 + iv] += wt2i;
846  w2[ic2 + id3 + iv] += wt2r;
847  w2[ic2 + 1 + id3 + iv] += wt2i;
848  w2[ic2 + id4 + iv] += -wt1r;
849  w2[ic2 + 1 + id4 + iv] += -wt1i;
850  }
851  }
852  }
853  }
854  }
855 
856 }
857 
858 //====================================================================
860  double *vcp1, double *v1, int ieo)
861 {
862  int Nvc2 = 2 * m_Nvc;
863  int Nvcd = m_Nvc * m_Nd;
864  int Nvcd2 = Nvcd / 2;
865 
866  int id1 = 0;
867  int id2 = m_Nvc;
868  int id3 = m_Nvc * 2;
869  int id4 = m_Nvc * 3;
870 
871  int isite = m_arg[itask].isite;
872  int isite_cp = m_arg[itask].isite_cpz;
873 
874  double *w2 = &vcp1[Nvcd2 * isite_cp];
875  double *w1 = &v1[Nvcd * isite];
876 
877  int idir = 2;
878  double bc2 = m_boundary2[idir];
879 
880  if (m_arg[itask].kz0 == 1) {
881  int Nxy = m_Nx2 * m_Ny;
882  int iz = 0;
883  for (int it = 0; it < m_Mt; ++it) {
884  for (int ixy = 0; ixy < Nxy; ++ixy) {
885  int is = ixy + Nxy * (iz + m_Nz * it);
886  int is2 = ixy + Nxy * it;
887 
888  int in = Nvcd * is;
889  int ix1 = Nvc2 * is2;
890  int ix2 = ix1 + m_Nvc;
891 
892  for (int ic = 0; ic < m_Nc; ++ic) {
893  w2[2 * ic + ix1] = bc2 * (w1[2 * ic + id1 + in] - w1[2 * ic + 1 + id3 + in]);
894  w2[2 * ic + 1 + ix1] = bc2 * (w1[2 * ic + 1 + id1 + in] + w1[2 * ic + id3 + in]);
895  w2[2 * ic + ix2] = bc2 * (w1[2 * ic + id2 + in] + w1[2 * ic + 1 + id4 + in]);
896  w2[2 * ic + 1 + ix2] = bc2 * (w1[2 * ic + 1 + id2 + in] - w1[2 * ic + id4 + in]);
897  }
898  }
899  }
900  }
901 
902 }
903 
904 //====================================================================
906  double *v2, double *vcp2, int ieo)
907 {
908  int Nvc2 = 2 * m_Nvc;
909  int Nvcd = m_Nvc * m_Nd;
910  int Nvcd2 = Nvcd / 2;
911 
912  int id1 = 0;
913  int id2 = m_Nvc;
914  int id3 = m_Nvc * 2;
915  int id4 = m_Nvc * 3;
916 
917  int idir = 2;
918 
919  double wt1r, wt1i, wt2r, wt2i;
920 
921  int isite = m_arg[itask].isite;
922  int isite_cp = m_arg[itask].isite_cpz;
923 
924  double *w2 = &v2[Nvcd * isite];
925  double *w1 = &vcp2[Nvcd2 * isite_cp];
926  double *u = const_cast<Field_G *>(m_U)->ptr(
927  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
928 
929  if (m_arg[itask].kz1 == 1) {
930  int Nxy = m_Nx2 * m_Ny;
931  int iz = m_Mz - 1;
932  for (int it = 0; it < m_Mt; ++it) {
933  for (int ixy = 0; ixy < Nxy; ++ixy) {
934  int is = ixy + Nxy * (iz + m_Nz * it);
935  int is2 = ixy + Nxy * it;
936  int iv = Nvcd * is;
937  int ig = m_Ndf * is;
938  int ix1 = Nvc2 * is2;
939  int ix2 = ix1 + m_Nvc;
940 
941  for (int ic = 0; ic < m_Nc; ++ic) {
942  int ic2 = ic * m_Nvc;
943 
944  wt1r = mult_uv_r(&u[ic2 + ig], &w1[ix1], m_Nc);
945  wt1i = mult_uv_i(&u[ic2 + ig], &w1[ix1], m_Nc);
946  wt2r = mult_uv_r(&u[ic2 + ig], &w1[ix2], m_Nc);
947  wt2i = mult_uv_i(&u[ic2 + ig], &w1[ix2], m_Nc);
948 
949  w2[2 * ic + id1 + iv] += wt1r;
950  w2[2 * ic + 1 + id1 + iv] += wt1i;
951  w2[2 * ic + id2 + iv] += wt2r;
952  w2[2 * ic + 1 + id2 + iv] += wt2i;
953  w2[2 * ic + id3 + iv] += wt1i;
954  w2[2 * ic + 1 + id3 + iv] += -wt1r;
955  w2[2 * ic + id4 + iv] += -wt2i;
956  w2[2 * ic + 1 + id4 + iv] += wt2r;
957  }
958  }
959  }
960  }
961 
962 }
963 
964 //====================================================================
966  double *v2, double *v1, int ieo)
967 {
968  int Nvcd = m_Nvc * m_Nd;
969 
970  int id1 = 0;
971  int id2 = m_Nvc;
972  int id3 = m_Nvc * 2;
973  int id4 = m_Nvc * 3;
974 
975  int idir = 2;
976 
977  double vt1[m_Nvc], vt2[m_Nvc];
978  double wt1r, wt1i, wt2r, wt2i;
979 
980  int isite = m_arg[itask].isite;
981 
982  double *w2 = &v2[Nvcd * isite];
983  double *w1 = &v1[Nvcd * isite];
984  double *u = const_cast<Field_G *>(m_U)->ptr(
985  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
986 
987  int kz1 = m_arg[itask].kz1;
988  int Nxy = m_Nx2 * m_Ny;
989 
990  for (int it = 0; it < m_Mt; ++it) {
991  for (int iz = 0; iz < m_Mz - kz1; ++iz) {
992  for (int ixy = 0; ixy < Nxy; ++ixy) {
993  int is = ixy + Nxy * (iz + m_Nz * it);
994  int iv = Nvcd * is;
995  int in = Nvcd * (is + Nxy);
996  int ig = m_Ndf * is;
997 
998  for (int ic = 0; ic < m_Nc; ++ic) {
999  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + 1 + id3 + in];
1000  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] + w1[2 * ic + id3 + in];
1001  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + 1 + id4 + in];
1002  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + id4 + in];
1003  }
1004 
1005  for (int ic = 0; ic < m_Nc; ++ic) {
1006  int ic2 = ic * m_Nvc;
1007 
1008  wt1r = mult_uv_r(&u[ic2 + ig], vt1, m_Nc);
1009  wt1i = mult_uv_i(&u[ic2 + ig], vt1, m_Nc);
1010  wt2r = mult_uv_r(&u[ic2 + ig], vt2, m_Nc);
1011  wt2i = mult_uv_i(&u[ic2 + ig], vt2, m_Nc);
1012 
1013  w2[2 * ic + id1 + iv] += wt1r;
1014  w2[2 * ic + 1 + id1 + iv] += wt1i;
1015  w2[2 * ic + id2 + iv] += wt2r;
1016  w2[2 * ic + 1 + id2 + iv] += wt2i;
1017  w2[2 * ic + id3 + iv] += wt1i;
1018  w2[2 * ic + 1 + id3 + iv] += -wt1r;
1019  w2[2 * ic + id4 + iv] += -wt2i;
1020  w2[2 * ic + 1 + id4 + iv] += wt2r;
1021  }
1022  }
1023  }
1024  }
1025 
1026 }
1027 
1028 //====================================================================
1030  double *vcp1, double *v1, int ieo)
1031 {
1032  int Nvc2 = 2 * m_Nvc;
1033  int Nvcd = m_Nvc * m_Nd;
1034  int Nvcd2 = Nvcd / 2;
1035 
1036  int id1 = 0;
1037  int id2 = m_Nvc;
1038  int id3 = m_Nvc * 2;
1039  int id4 = m_Nvc * 3;
1040 
1041  int idir = 2;
1042 
1043  int isite = m_arg[itask].isite;
1044  int isite_cp = m_arg[itask].isite_cpz;
1045 
1046  double *w2 = &vcp1[Nvcd2 * isite_cp];
1047  double *w1 = &v1[Nvcd * isite];
1048  double *u = const_cast<Field_G *>(m_U)->ptr(
1049  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1050 
1051  double vt1[m_Nvc], vt2[m_Nvc];
1052 
1053  if (m_arg[itask].kz1 == 1) {
1054  int Nxy = m_Nx2 * m_Ny;
1055  int iz = m_Mz - 1;
1056  for (int it = 0; it < m_Mt; ++it) {
1057  for (int ixy = 0; ixy < Nxy; ++ixy) {
1058  int is = ixy + Nxy * (iz + m_Nz * it);
1059  int is2 = ixy + Nxy * it;
1060  int in = Nvcd * is;
1061  int ig = m_Ndf * is;
1062  int ix1 = Nvc2 * is2;
1063  int ix2 = ix1 + m_Nvc;
1064 
1065  for (int ic = 0; ic < m_Nc; ++ic) {
1066  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + 1 + id3 + in];
1067  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + id3 + in];
1068  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + 1 + id4 + in];
1069  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + id4 + in];
1070  }
1071 
1072  for (int ic = 0; ic < m_Nc; ++ic) {
1073  int icr = 2 * ic;
1074  w2[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1, m_Nc);
1075  w2[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1, m_Nc);
1076  w2[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2, m_Nc);
1077  w2[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2, m_Nc);
1078  }
1079  }
1080  }
1081  }
1082 
1083 }
1084 
1085 //====================================================================
1087  double *v2, double *vcp2, int ieo)
1088 {
1089  int Nvc2 = 2 * m_Nvc;
1090  int Nvcd = m_Nvc * m_Nd;
1091  int Nvcd2 = Nvcd / 2;
1092 
1093  int id1 = 0;
1094  int id2 = m_Nvc;
1095  int id3 = m_Nvc * 2;
1096  int id4 = m_Nvc * 3;
1097 
1098  int idir = 2;
1099  double bc2 = m_boundary2[idir];
1100 
1101  // double wt1r, wt1i, wt2r, wt2i;
1102 
1103  int isite = m_arg[itask].isite;
1104  int isite_cp = m_arg[itask].isite_cpz;
1105 
1106  double *w2 = &v2[Nvcd * isite];
1107  double *w1 = &vcp2[Nvcd2 * isite_cp];
1108 
1109  if (m_arg[itask].kz0 == 1) {
1110  int Nxy = m_Nx2 * m_Ny;
1111 
1112  int iz = 0;
1113  for (int it = 0; it < m_Mt; ++it) {
1114  for (int ixy = 0; ixy < Nxy; ++ixy) {
1115  int is = ixy + Nxy * (iz + m_Nz * it);
1116  int is2 = ixy + Nxy * it;
1117  int iv = Nvcd * is;
1118  int ix1 = Nvc2 * is2;
1119  int ix2 = ix1 + m_Nvc;
1120 
1121  for (int ic = 0; ic < m_Nc; ++ic) {
1122  int icr = 2 * ic;
1123  int ici = 2 * ic + 1;
1124  w2[icr + id1 + iv] += bc2 * w1[icr + ix1];
1125  w2[ici + id1 + iv] += bc2 * w1[ici + ix1];
1126  w2[icr + id2 + iv] += bc2 * w1[icr + ix2];
1127  w2[ici + id2 + iv] += bc2 * w1[ici + ix2];
1128  w2[icr + id3 + iv] += -bc2 * w1[ici + ix1];
1129  w2[ici + id3 + iv] += bc2 * w1[icr + ix1];
1130  w2[icr + id4 + iv] += bc2 * w1[ici + ix2];
1131  w2[ici + id4 + iv] += -bc2 * w1[icr + ix2];
1132  }
1133  }
1134  }
1135  }
1136 
1137 }
1138 
1139 //====================================================================
1141  double *v2, double *v1, int ieo)
1142 {
1143  int Nvcd = m_Nvc * m_Nd;
1144 
1145  int id1 = 0;
1146  int id2 = m_Nvc;
1147  int id3 = m_Nvc * 2;
1148  int id4 = m_Nvc * 3;
1149 
1150  int idir = 2;
1151 
1152  double vt1[m_Nvc], vt2[m_Nvc];
1153  double wt1r, wt1i, wt2r, wt2i;
1154 
1155  int isite = m_arg[itask].isite;
1156 
1157  double *w2 = &v2[Nvcd * isite];
1158  double *w1 = &v1[Nvcd * isite];
1159  double *u = const_cast<Field_G *>(m_U)->ptr(
1160  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1161 
1162  int kz0 = m_arg[itask].kz0;
1163  int Nxy = m_Nx2 * m_Ny;
1164 
1165  for (int it = 0; it < m_Mt; ++it) {
1166  for (int iz = kz0; iz < m_Mz; ++iz) {
1167  for (int ixy = 0; ixy < Nxy; ++ixy) {
1168  int is = ixy + Nxy * (iz + m_Nz * it);
1169  int iv = Nvcd * is;
1170  int in = Nvcd * (is - Nxy);
1171  int ig = m_Ndf * (is - Nxy);
1172 
1173  for (int ic = 0; ic < m_Nc; ++ic) {
1174  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + 1 + id3 + in];
1175  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + id3 + in];
1176  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + 1 + id4 + in];
1177  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + id4 + in];
1178  }
1179 
1180  for (int ic = 0; ic < m_Nc; ++ic) {
1181  int ic2 = 2 * ic;
1182  wt1r = mult_udagv_r(&u[ic2 + ig], vt1, m_Nc);
1183  wt1i = mult_udagv_i(&u[ic2 + ig], vt1, m_Nc);
1184  wt2r = mult_udagv_r(&u[ic2 + ig], vt2, m_Nc);
1185  wt2i = mult_udagv_i(&u[ic2 + ig], vt2, m_Nc);
1186 
1187  w2[ic2 + id1 + iv] += wt1r;
1188  w2[ic2 + 1 + id1 + iv] += wt1i;
1189  w2[ic2 + id2 + iv] += wt2r;
1190  w2[ic2 + 1 + id2 + iv] += wt2i;
1191  w2[ic2 + id3 + iv] += -wt1i;
1192  w2[ic2 + 1 + id3 + iv] += wt1r;
1193  w2[ic2 + id4 + iv] += wt2i;
1194  w2[ic2 + 1 + id4 + iv] += -wt2r;
1195  }
1196  }
1197  }
1198  }
1199 
1200 }
1201 
1202 //====================================================================
1204  double *vcp1, double *v1, int ieo)
1205 {
1206  int Nvc2 = 2 * m_Nvc;
1207  int Nvcd = m_Nvc * m_Nd;
1208  int Nvcd2 = Nvcd / 2;
1209 
1210  // int id1 = 0;
1211  // int id2 = m_Nvc;
1212  int id3 = m_Nvc * 2;
1213  int id4 = m_Nvc * 3;
1214 
1215  int isite = m_arg[itask].isite;
1216  int isite_cp = m_arg[itask].isite_cpt;
1217 
1218  double *w2 = &vcp1[Nvcd2 * isite_cp];
1219  double *w1 = &v1[Nvcd * isite];
1220 
1221  int idir = 3;
1222  double bc2 = m_boundary2[idir];
1223 
1224  if (m_arg[itask].kt0 == 1) {
1225  int Nxy = m_Nx2 * m_Ny;
1226  int it = 0;
1227  for (int iz = 0; iz < m_Mz; ++iz) {
1228  for (int ixy = 0; ixy < Nxy; ++ixy) {
1229  int is = ixy + Nxy * (iz + m_Nz * it);
1230  int is2 = ixy + Nxy * iz;
1231 
1232  int in = Nvcd * is;
1233  int ix1 = Nvc2 * is2;
1234  int ix2 = ix1 + m_Nvc;
1235 
1236  for (int ic = 0; ic < m_Nc; ++ic) {
1237  w2[2 * ic + ix1] = 2.0 * bc2 * w1[2 * ic + id3 + in];
1238  w2[2 * ic + 1 + ix1] = 2.0 * bc2 * w1[2 * ic + 1 + id3 + in];
1239  w2[2 * ic + ix2] = 2.0 * bc2 * w1[2 * ic + id4 + in];
1240  w2[2 * ic + 1 + ix2] = 2.0 * bc2 * w1[2 * ic + 1 + id4 + in];
1241  }
1242  }
1243  }
1244  }
1245 
1246 }
1247 
1248 //====================================================================
1250  double *v2, double *vcp2, int ieo)
1251 {
1252  int Nvc2 = 2 * m_Nvc;
1253  int Nvcd = m_Nvc * m_Nd;
1254  int Nvcd2 = Nvcd / 2;
1255 
1256  // int id1 = 0;
1257  // int id2 = m_Nvc;
1258  int id3 = m_Nvc * 2;
1259  int id4 = m_Nvc * 3;
1260 
1261  int idir = 3;
1262 
1263  double wt1r, wt1i, wt2r, wt2i;
1264 
1265  int isite = m_arg[itask].isite;
1266  int isite_cp = m_arg[itask].isite_cpt;
1267 
1268  double *w2 = &v2[Nvcd * isite];
1269  double *w1 = &vcp2[Nvcd2 * isite_cp];
1270  double *u = const_cast<Field_G *>(m_U)->ptr(
1271  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
1272 
1273  if (m_arg[itask].kt1 == 1) {
1274  int Nxy = m_Nx2 * m_Ny;
1275  int it = m_Mt - 1;
1276  for (int iz = 0; iz < m_Mz; ++iz) {
1277  for (int ixy = 0; ixy < Nxy; ++ixy) {
1278  int is = ixy + Nxy * (iz + m_Nz * it);
1279  int is2 = ixy + Nxy * iz;
1280  int iv = Nvcd * is;
1281  int ig = m_Ndf * is;
1282  int ix1 = Nvc2 * is2;
1283  int ix2 = ix1 + m_Nvc;
1284 
1285  for (int ic = 0; ic < m_Nc; ++ic) {
1286  int ic2 = ic * m_Nvc;
1287 
1288  wt1r = mult_uv_r(&u[ic2 + ig], &w1[ix1], m_Nc);
1289  wt1i = mult_uv_i(&u[ic2 + ig], &w1[ix1], m_Nc);
1290  wt2r = mult_uv_r(&u[ic2 + ig], &w1[ix2], m_Nc);
1291  wt2i = mult_uv_i(&u[ic2 + ig], &w1[ix2], m_Nc);
1292 
1293  w2[2 * ic + id3 + iv] += wt1r;
1294  w2[2 * ic + 1 + id3 + iv] += wt1i;
1295  w2[2 * ic + id4 + iv] += wt2r;
1296  w2[2 * ic + 1 + id4 + iv] += wt2i;
1297  }
1298  }
1299  }
1300  }
1301 
1302 }
1303 
1304 //====================================================================
1306  double *v2, double *v1, int ieo)
1307 {
1308  int Nvcd = m_Nvc * m_Nd;
1309 
1310  // int id1 = 0;
1311  // int id2 = m_Nvc;
1312  int id3 = m_Nvc * 2;
1313  int id4 = m_Nvc * 3;
1314 
1315  int idir = 3;
1316 
1317  double vt1[m_Nvc], vt2[m_Nvc];
1318  double wt1r, wt1i, wt2r, wt2i;
1319 
1320  int isite = m_arg[itask].isite;
1321 
1322  double *w2 = &v2[Nvcd * isite];
1323  double *w1 = &v1[Nvcd * isite];
1324  double *u = const_cast<Field_G *>(m_U)->ptr(
1325  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
1326 
1327  int kt1 = m_arg[itask].kt1;
1328  int Nxy = m_Nx2 * m_Ny;
1329  int Nxyz = Nxy * m_Nz;
1330 
1331  for (int it = 0; it < m_Mt - kt1; ++it) {
1332  for (int iz = 0; iz < m_Mz; ++iz) {
1333  for (int ixy = 0; ixy < Nxy; ++ixy) {
1334  int is = ixy + Nxy * (iz + m_Nz * it);
1335  int iv = Nvcd * is;
1336  int in = Nvcd * (is + Nxyz);
1337  int ig = m_Ndf * is;
1338 
1339  for (int ic = 0; ic < m_Nc; ++ic) {
1340  vt1[2 * ic] = 2.0 * w1[2 * ic + id3 + in];
1341  vt1[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id3 + in];
1342  vt2[2 * ic] = 2.0 * w1[2 * ic + id4 + in];
1343  vt2[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id4 + in];
1344  }
1345 
1346  for (int ic = 0; ic < m_Nc; ++ic) {
1347  int ic2 = ic * m_Nvc;
1348 
1349  wt1r = mult_uv_r(&u[ic2 + ig], vt1, m_Nc);
1350  wt1i = mult_uv_i(&u[ic2 + ig], vt1, m_Nc);
1351  wt2r = mult_uv_r(&u[ic2 + ig], vt2, m_Nc);
1352  wt2i = mult_uv_i(&u[ic2 + ig], vt2, m_Nc);
1353 
1354  w2[2 * ic + id3 + iv] += wt1r;
1355  w2[2 * ic + 1 + id3 + iv] += wt1i;
1356  w2[2 * ic + id4 + iv] += wt2r;
1357  w2[2 * ic + 1 + id4 + iv] += wt2i;
1358  }
1359  }
1360  }
1361  }
1362 
1363 }
1364 
1365 //====================================================================
1367  double *vcp1, double *v1, int ieo)
1368 {
1369  int Nvc2 = 2 * m_Nvc;
1370  int Nvcd = m_Nvc * m_Nd;
1371  int Nvcd2 = Nvcd / 2;
1372 
1373  int id1 = 0;
1374  int id2 = m_Nvc;
1375  // int id3 = m_Nvc * 2;
1376  // int id4 = m_Nvc * 3;
1377 
1378  int idir = 3;
1379 
1380  int isite = m_arg[itask].isite;
1381  int isite_cp = m_arg[itask].isite_cpt;
1382 
1383  double *w2 = &vcp1[Nvcd2 * isite_cp];
1384  double *w1 = &v1[Nvcd * isite];
1385  double *u = const_cast<Field_G *>(m_U)->ptr(
1386  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1387 
1388  double vt1[m_Nvc], vt2[m_Nvc];
1389 
1390  if (m_arg[itask].kt1 == 1) {
1391  int Nxy = m_Nx2 * m_Ny;
1392  int it = m_Mt - 1;
1393  for (int iz = 0; iz < m_Mz; ++iz) {
1394  for (int ixy = 0; ixy < Nxy; ++ixy) {
1395  int is = ixy + Nxy * (iz + m_Nz * it);
1396  int is2 = ixy + Nxy * iz;
1397  int in = Nvcd * is;
1398  int ig = m_Ndf * is;
1399  int ix1 = Nvc2 * is2;
1400  int ix2 = ix1 + m_Nvc;
1401 
1402  for (int ic = 0; ic < m_Nc; ++ic) {
1403  vt1[2 * ic] = 2.0 * w1[2 * ic + id1 + in];
1404  vt1[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id1 + in];
1405  vt2[2 * ic] = 2.0 * w1[2 * ic + id2 + in];
1406  vt2[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id2 + in];
1407  }
1408 
1409  for (int ic = 0; ic < m_Nc; ++ic) {
1410  int icr = 2 * ic;
1411  w2[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1, m_Nc);
1412  w2[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1, m_Nc);
1413  w2[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2, m_Nc);
1414  w2[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2, m_Nc);
1415  }
1416  }
1417  }
1418  }
1419 
1420 }
1421 
1422 //====================================================================
1424  double *v2, double *vcp2, int ieo)
1425 {
1426  int Nvc2 = 2 * m_Nvc;
1427  int Nvcd = m_Nvc * m_Nd;
1428  int Nvcd2 = Nvcd / 2;
1429 
1430  int id1 = 0;
1431  int id2 = m_Nvc;
1432  // int id3 = m_Nvc * 2;
1433  // int id4 = m_Nvc * 3;
1434 
1435  int idir = 3;
1436  double bc2 = m_boundary2[idir];
1437 
1438  // double wt1r, wt1i, wt2r, wt2i;
1439 
1440  int isite = m_arg[itask].isite;
1441  int isite_cp = m_arg[itask].isite_cpt;
1442 
1443  double *w2 = &v2[Nvcd * isite];
1444  double *w1 = &vcp2[Nvcd2 * isite_cp];
1445 
1446  if (m_arg[itask].kt0 == 1) {
1447  int Nxy = m_Nx2 * m_Ny;
1448  int it = 0;
1449  for (int iz = 0; iz < m_Mz; ++iz) {
1450  for (int ixy = 0; ixy < Nxy; ++ixy) {
1451  int is = ixy + Nxy * (iz + m_Nz * it);
1452  int is2 = ixy + Nxy * iz;
1453  int iv = Nvcd * is;
1454  int ix1 = Nvc2 * is2;
1455  int ix2 = ix1 + m_Nvc;
1456 
1457  for (int ic = 0; ic < m_Nc; ++ic) {
1458  int icr = 2 * ic;
1459  int ici = 2 * ic + 1;
1460  w2[icr + id1 + iv] += bc2 * w1[icr + ix1];
1461  w2[ici + id1 + iv] += bc2 * w1[ici + ix1];
1462  w2[icr + id2 + iv] += bc2 * w1[icr + ix2];
1463  w2[ici + id2 + iv] += bc2 * w1[ici + ix2];
1464  }
1465  }
1466  }
1467  }
1468 
1469 }
1470 
1471 //====================================================================
1473  double *v2, double *v1, int ieo)
1474 {
1475  int Nvcd = m_Nvc * m_Nd;
1476 
1477  int id1 = 0;
1478  int id2 = m_Nvc;
1479  // int id3 = m_Nvc * 2;
1480  // int id4 = m_Nvc * 3;
1481 
1482  int idir = 3;
1483 
1484  double vt1[m_Nvc], vt2[m_Nvc];
1485  double wt1r, wt1i, wt2r, wt2i;
1486 
1487  int isite = m_arg[itask].isite;
1488 
1489  double *w2 = &v2[Nvcd * isite];
1490  double *w1 = &v1[Nvcd * isite];
1491  double *u = const_cast<Field_G *>(m_U)->ptr(
1492  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1493 
1494  int kt0 = m_arg[itask].kt0;
1495  int Nxy = m_Nx2 * m_Ny;
1496  int Nxyz = Nxy * m_Nz;
1497 
1498  for (int it = kt0; it < m_Mt; ++it) {
1499  for (int iz = 0; iz < m_Mz; ++iz) {
1500  for (int ixy = 0; ixy < Nxy; ++ixy) {
1501  int is = ixy + Nxy * (iz + m_Nz * it);
1502  int iv = Nvcd * is;
1503  int in = Nvcd * (is - Nxyz);
1504  int ig = m_Ndf * (is - Nxyz);
1505 
1506  for (int ic = 0; ic < m_Nc; ++ic) {
1507  vt1[2 * ic] = 2.0 * w1[2 * ic + id1 + in];
1508  vt1[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id1 + in];
1509  vt2[2 * ic] = 2.0 * w1[2 * ic + id2 + in];
1510  vt2[2 * ic + 1] = 2.0 * w1[2 * ic + 1 + id2 + in];
1511  }
1512 
1513  for (int ic = 0; ic < m_Nc; ++ic) {
1514  int ic2 = 2 * ic;
1515  wt1r = mult_udagv_r(&u[ic2 + ig], vt1, m_Nc);
1516  wt1i = mult_udagv_i(&u[ic2 + ig], vt1, m_Nc);
1517  wt2r = mult_udagv_r(&u[ic2 + ig], vt2, m_Nc);
1518  wt2i = mult_udagv_i(&u[ic2 + ig], vt2, m_Nc);
1519 
1520  w2[ic2 + id1 + iv] += wt1r;
1521  w2[ic2 + 1 + id1 + iv] += wt1i;
1522  w2[ic2 + id2 + iv] += wt2r;
1523  w2[ic2 + 1 + id2 + iv] += wt2i;
1524  }
1525  }
1526  }
1527  }
1528 
1529 }
1530 
1531 //====================================================================
1533  double *vcp1, double *v1, int ieo)
1534 {
1535  int Nvc2 = 2 * m_Nvc;
1536  int Nvcd = m_Nvc * m_Nd;
1537  int Nvcd2 = Nvcd / 2;
1538 
1539  int id1 = 0;
1540  int id2 = m_Nvc;
1541  int id3 = m_Nvc * 2;
1542  int id4 = m_Nvc * 3;
1543 
1544  int isite = m_arg[itask].isite;
1545  int isite_cp = m_arg[itask].isite_cpt;
1546 
1547  double *w2 = &vcp1[Nvcd2 * isite_cp];
1548  double *w1 = &v1[Nvcd * isite];
1549 
1550  int idir = 3;
1551  double bc2 = m_boundary2[idir];
1552 
1553  if (m_arg[itask].kt0 == 1) {
1554  int Nxy = m_Nx2 * m_Ny;
1555  int it = 0;
1556  for (int iz = 0; iz < m_Mz; ++iz) {
1557  for (int ixy = 0; ixy < Nxy; ++ixy) {
1558  int is = ixy + Nxy * (iz + m_Nz * it);
1559  int is2 = ixy + Nxy * iz;
1560 
1561  int in = Nvcd * is;
1562  int ix1 = Nvc2 * is2;
1563  int ix2 = ix1 + m_Nvc;
1564 
1565  for (int ic = 0; ic < m_Nc; ++ic) {
1566  w2[2 * ic + ix1] = bc2 * (w1[2 * ic + id1 + in] + w1[2 * ic + id3 + in]);
1567  w2[2 * ic + 1 + ix1] = bc2 * (w1[2 * ic + 1 + id1 + in] + w1[2 * ic + 1 + id3 + in]);
1568  w2[2 * ic + ix2] = bc2 * (w1[2 * ic + id2 + in] + w1[2 * ic + id4 + in]);
1569  w2[2 * ic + 1 + ix2] = bc2 * (w1[2 * ic + 1 + id2 + in] + w1[2 * ic + 1 + id4 + in]);
1570  }
1571  }
1572  }
1573  }
1574 
1575 }
1576 
1577 //====================================================================
1579  double *v2, double *vcp2, int ieo)
1580 {
1581  int Nvc2 = 2 * m_Nvc;
1582  int Nvcd = m_Nvc * m_Nd;
1583  int Nvcd2 = Nvcd / 2;
1584 
1585  int id1 = 0;
1586  int id2 = m_Nvc;
1587  int id3 = m_Nvc * 2;
1588  int id4 = m_Nvc * 3;
1589 
1590  int idir = 3;
1591 
1592  double wt1r, wt1i, wt2r, wt2i;
1593 
1594  int isite = m_arg[itask].isite;
1595  int isite_cp = m_arg[itask].isite_cpt;
1596 
1597  double *w2 = &v2[Nvcd * isite];
1598  double *w1 = &vcp2[Nvcd2 * isite_cp];
1599  double *u = const_cast<Field_G *>(m_U)->ptr(
1600  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
1601 
1602  if (m_arg[itask].kt1 == 1) {
1603  int Nxy = m_Nx2 * m_Ny;
1604  int it = m_Mt - 1;
1605  for (int iz = 0; iz < m_Mz; ++iz) {
1606  for (int ixy = 0; ixy < Nxy; ++ixy) {
1607  int is = ixy + Nxy * (iz + m_Nz * it);
1608  int is2 = ixy + Nxy * iz;
1609  int iv = Nvcd * is;
1610  int ig = m_Ndf * is;
1611  int ix1 = Nvc2 * is2;
1612  int ix2 = ix1 + m_Nvc;
1613 
1614  for (int ic = 0; ic < m_Nc; ++ic) {
1615  int ic2 = ic * m_Nvc;
1616 
1617  wt1r = mult_uv_r(&u[ic2 + ig], &w1[ix1], m_Nc);
1618  wt1i = mult_uv_i(&u[ic2 + ig], &w1[ix1], m_Nc);
1619  wt2r = mult_uv_r(&u[ic2 + ig], &w1[ix2], m_Nc);
1620  wt2i = mult_uv_i(&u[ic2 + ig], &w1[ix2], m_Nc);
1621 
1622  w2[2 * ic + id1 + iv] += wt1r;
1623  w2[2 * ic + 1 + id1 + iv] += wt1i;
1624  w2[2 * ic + id2 + iv] += wt2r;
1625  w2[2 * ic + 1 + id2 + iv] += wt2i;
1626  w2[2 * ic + id3 + iv] += wt1r;
1627  w2[2 * ic + 1 + id3 + iv] += wt1i;
1628  w2[2 * ic + id4 + iv] += wt2r;
1629  w2[2 * ic + 1 + id4 + iv] += wt2i;
1630  }
1631  }
1632  }
1633  }
1634 
1635 }
1636 
1637 //====================================================================
1639  double *v2, double *v1, int ieo)
1640 {
1641  int Nvcd = m_Nvc * m_Nd;
1642 
1643  int id1 = 0;
1644  int id2 = m_Nvc;
1645  int id3 = m_Nvc * 2;
1646  int id4 = m_Nvc * 3;
1647 
1648  int idir = 3;
1649 
1650  double vt1[m_Nvc], vt2[m_Nvc];
1651  double wt1r, wt1i, wt2r, wt2i;
1652 
1653  int isite = m_arg[itask].isite;
1654 
1655  double *w2 = &v2[Nvcd * isite];
1656  double *w1 = &v1[Nvcd * isite];
1657  double *u = const_cast<Field_G *>(m_U)->ptr(
1658  m_Ndf * (isite + ieo * m_Nvol / 2 + idir * m_Nvol));
1659 
1660  int kt1 = m_arg[itask].kt1;
1661  int Nxy = m_Nx2 * m_Ny;
1662  int Nxyz = Nxy * m_Nz;
1663 
1664  for (int it = 0; it < m_Mt - kt1; ++it) {
1665  for (int iz = 0; iz < m_Mz; ++iz) {
1666  for (int ixy = 0; ixy < Nxy; ++ixy) {
1667  int is = ixy + Nxy * (iz + m_Nz * it);
1668  int iv = Nvcd * is;
1669  int in = Nvcd * (is + Nxyz);
1670  int ig = m_Ndf * is;
1671 
1672  for (int ic = 0; ic < m_Nc; ++ic) {
1673  vt1[2 * ic] = w1[2 * ic + id1 + in] + w1[2 * ic + id3 + in];
1674  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] + w1[2 * ic + 1 + id3 + in];
1675  vt2[2 * ic] = w1[2 * ic + id2 + in] + w1[2 * ic + id4 + in];
1676  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] + w1[2 * ic + 1 + id4 + in];
1677  }
1678 
1679  for (int ic = 0; ic < m_Nc; ++ic) {
1680  int ic2 = ic * m_Nvc;
1681 
1682  wt1r = mult_uv_r(&u[ic2 + ig], vt1, m_Nc);
1683  wt1i = mult_uv_i(&u[ic2 + ig], vt1, m_Nc);
1684  wt2r = mult_uv_r(&u[ic2 + ig], vt2, m_Nc);
1685  wt2i = mult_uv_i(&u[ic2 + ig], vt2, m_Nc);
1686 
1687  w2[2 * ic + id1 + iv] += wt1r;
1688  w2[2 * ic + 1 + id1 + iv] += wt1i;
1689  w2[2 * ic + id2 + iv] += wt2r;
1690  w2[2 * ic + 1 + id2 + iv] += wt2i;
1691  w2[2 * ic + id3 + iv] += wt1r;
1692  w2[2 * ic + 1 + id3 + iv] += wt1i;
1693  w2[2 * ic + id4 + iv] += wt2r;
1694  w2[2 * ic + 1 + id4 + iv] += wt2i;
1695  }
1696  }
1697  }
1698  }
1699 
1700 }
1701 
1702 //====================================================================
1704  double *vcp1, double *v1, int ieo)
1705 {
1706  int Nvc2 = 2 * m_Nvc;
1707  int Nvcd = m_Nvc * m_Nd;
1708  int Nvcd2 = Nvcd / 2;
1709 
1710  int id1 = 0;
1711  int id2 = m_Nvc;
1712  int id3 = m_Nvc * 2;
1713  int id4 = m_Nvc * 3;
1714 
1715  int idir = 3;
1716 
1717  int isite = m_arg[itask].isite;
1718  int isite_cp = m_arg[itask].isite_cpt;
1719 
1720  double *w2 = &vcp1[Nvcd2 * isite_cp];
1721  double *w1 = &v1[Nvcd * isite];
1722  double *u = const_cast<Field_G *>(m_U)->ptr(
1723  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1724 
1725  double vt1[m_Nvc], vt2[m_Nvc];
1726 
1727  if (m_arg[itask].kt1 == 1) {
1728  int Nxy = m_Nx2 * m_Ny;
1729  int it = m_Mt - 1;
1730  for (int iz = 0; iz < m_Mz; ++iz) {
1731  for (int ixy = 0; ixy < Nxy; ++ixy) {
1732  int is = ixy + Nxy * (iz + m_Nz * it);
1733  int is2 = ixy + Nxy * iz;
1734  int in = Nvcd * is;
1735  int ig = m_Ndf * is;
1736  int ix1 = Nvc2 * is2;
1737  int ix2 = ix1 + m_Nvc;
1738 
1739  for (int ic = 0; ic < m_Nc; ++ic) {
1740  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + id3 + in];
1741  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + 1 + id3 + in];
1742  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + id4 + in];
1743  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + 1 + id4 + in];
1744  }
1745 
1746  for (int ic = 0; ic < m_Nc; ++ic) {
1747  int icr = 2 * ic;
1748  w2[icr + ix1] = mult_udagv_r(&u[icr + ig], vt1, m_Nc);
1749  w2[icr + 1 + ix1] = mult_udagv_i(&u[icr + ig], vt1, m_Nc);
1750  w2[icr + ix2] = mult_udagv_r(&u[icr + ig], vt2, m_Nc);
1751  w2[icr + 1 + ix2] = mult_udagv_i(&u[icr + ig], vt2, m_Nc);
1752  }
1753  }
1754  }
1755  }
1756 
1757 }
1758 
1759 //====================================================================
1761  double *v2, double *vcp2, int ieo)
1762 {
1763  int Nvc2 = 2 * m_Nvc;
1764  int Nvcd = m_Nvc * m_Nd;
1765  int Nvcd2 = Nvcd / 2;
1766 
1767  int id1 = 0;
1768  int id2 = m_Nvc;
1769  int id3 = m_Nvc * 2;
1770  int id4 = m_Nvc * 3;
1771 
1772  int idir = 3;
1773  double bc2 = m_boundary2[idir];
1774 
1775  // double wt1r, wt1i, wt2r, wt2i;
1776 
1777  int isite = m_arg[itask].isite;
1778  int isite_cp = m_arg[itask].isite_cpt;
1779 
1780  double *w2 = &v2[Nvcd * isite];
1781  double *w1 = &vcp2[Nvcd2 * isite_cp];
1782 
1783  if (m_arg[itask].kt0 == 1) {
1784  int Nxy = m_Nx2 * m_Ny;
1785  int it = 0;
1786  for (int iz = 0; iz < m_Mz; ++iz) {
1787  for (int ixy = 0; ixy < Nxy; ++ixy) {
1788  int is = ixy + Nxy * (iz + m_Nz * it);
1789  int is2 = ixy + Nxy * iz;
1790  int iv = Nvcd * is;
1791  int ix1 = Nvc2 * is2;
1792  int ix2 = ix1 + m_Nvc;
1793 
1794  for (int ic = 0; ic < m_Nc; ++ic) {
1795  int icr = 2 * ic;
1796  int ici = 2 * ic + 1;
1797  w2[icr + id1 + iv] += bc2 * w1[icr + ix1];
1798  w2[ici + id1 + iv] += bc2 * w1[ici + ix1];
1799  w2[icr + id2 + iv] += bc2 * w1[icr + ix2];
1800  w2[ici + id2 + iv] += bc2 * w1[ici + ix2];
1801  w2[icr + id3 + iv] -= bc2 * w1[icr + ix1];
1802  w2[ici + id3 + iv] -= bc2 * w1[ici + ix1];
1803  w2[icr + id4 + iv] -= bc2 * w1[icr + ix2];
1804  w2[ici + id4 + iv] -= bc2 * w1[ici + ix2];
1805  }
1806  }
1807  }
1808  }
1809 
1810 }
1811 
1812 //====================================================================
1814  double *v2, double *v1, int ieo)
1815 {
1816  int Nvcd = m_Nvc * m_Nd;
1817 
1818  int id1 = 0;
1819  int id2 = m_Nvc;
1820  int id3 = m_Nvc * 2;
1821  int id4 = m_Nvc * 3;
1822 
1823  int idir = 3;
1824 
1825  double vt1[m_Nvc], vt2[m_Nvc];
1826  double wt1r, wt1i, wt2r, wt2i;
1827 
1828  int isite = m_arg[itask].isite;
1829 
1830  double *w2 = &v2[Nvcd * isite];
1831  double *w1 = &v1[Nvcd * isite];
1832  double *u = const_cast<Field_G *>(m_U)->ptr(
1833  m_Ndf * (isite + (1 - ieo) * m_Nvol / 2 + idir * m_Nvol));
1834 
1835  int kt0 = m_arg[itask].kt0;
1836  int Nxy = m_Nx2 * m_Ny;
1837  int Nxyz = Nxy * m_Nz;
1838 
1839  for (int it = kt0; it < m_Mt; ++it) {
1840  for (int iz = 0; iz < m_Mz; ++iz) {
1841  for (int ixy = 0; ixy < Nxy; ++ixy) {
1842  int is = ixy + Nxy * (iz + m_Nz * it);
1843  int iv = Nvcd * is;
1844  int in = Nvcd * (is - Nxyz);
1845  int ig = m_Ndf * (is - Nxyz);
1846 
1847  for (int ic = 0; ic < m_Nc; ++ic) {
1848  vt1[2 * ic] = w1[2 * ic + id1 + in] - w1[2 * ic + id3 + in];
1849  vt1[2 * ic + 1] = w1[2 * ic + 1 + id1 + in] - w1[2 * ic + 1 + id3 + in];
1850  vt2[2 * ic] = w1[2 * ic + id2 + in] - w1[2 * ic + id4 + in];
1851  vt2[2 * ic + 1] = w1[2 * ic + 1 + id2 + in] - w1[2 * ic + 1 + id4 + in];
1852  }
1853 
1854  for (int ic = 0; ic < m_Nc; ++ic) {
1855  int ic2 = 2 * ic;
1856  wt1r = mult_udagv_r(&u[ic2 + ig], vt1, m_Nc);
1857  wt1i = mult_udagv_i(&u[ic2 + ig], vt1, m_Nc);
1858  wt2r = mult_udagv_r(&u[ic2 + ig], vt2, m_Nc);
1859  wt2i = mult_udagv_i(&u[ic2 + ig], vt2, m_Nc);
1860 
1861  w2[ic2 + id1 + iv] += wt1r;
1862  w2[ic2 + 1 + id1 + iv] += wt1i;
1863  w2[ic2 + id2 + iv] += wt2r;
1864  w2[ic2 + 1 + id2 + iv] += wt2i;
1865  w2[ic2 + id3 + iv] -= wt1r;
1866  w2[ic2 + 1 + id3 + iv] -= wt1i;
1867  w2[ic2 + id4 + iv] -= wt2r;
1868  w2[ic2 + 1 + id4 + iv] -= wt2i;
1869  }
1870  }
1871  }
1872  }
1873 
1874 }
1875 
1876 //====================================================================
1878  double *v2, double *v1)
1879 {
1880  int Nvcd = m_Nvc * m_Nd;
1881  int Nxy = m_Nx2 * m_Ny;
1882 
1883  int id1 = 0;
1884  int id2 = m_Nvc;
1885  int id3 = m_Nvc * 2;
1886  int id4 = m_Nvc * 3;
1887 
1888  int isite = m_arg[itask].isite;
1889  double *w2 = &v2[Nvcd * isite];
1890  double *w1 = &v1[Nvcd * isite];
1891 
1892  for (int it = 0; it < m_Mt; ++it) {
1893  for (int iz = 0; iz < m_Mz; ++iz) {
1894  for (int ixy = 0; ixy < Nxy; ++ixy) {
1895  int iv = Nvcd * (ixy + Nxy * (iz + m_Nz * it));
1896  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
1897  w2[ivc + id1 + iv] = w1[ivc + id3 + iv];
1898  w2[ivc + id2 + iv] = w1[ivc + id4 + iv];
1899  w2[ivc + id3 + iv] = w1[ivc + id1 + iv];
1900  w2[ivc + id4 + iv] = w1[ivc + id2 + iv];
1901  }
1902  }
1903  }
1904  }
1905 
1906 }
1907 
1908 //====================================================================
1910  double *v2, double *v1)
1911 {
1912  int Nvcd = m_Nvc * m_Nd;
1913  int Nxy = m_Nx2 * m_Ny;
1914 
1915  int id1 = 0;
1916  int id2 = m_Nvc;
1917  int id3 = m_Nvc * 2;
1918  int id4 = m_Nvc * 3;
1919 
1920  int isite = m_arg[itask].isite;
1921  double *w2 = &v2[Nvcd * isite];
1922  double *w1 = &v1[Nvcd * isite];
1923 
1924  for (int it = 0; it < m_Mt; ++it) {
1925  for (int iz = 0; iz < m_Mz; ++iz) {
1926  for (int ixy = 0; ixy < Nxy; ++ixy) {
1927  int iv = Nvcd * (ixy + Nxy * (iz + m_Nz * it));
1928  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
1929  w2[ivc + id1 + iv] = w1[ivc + id1 + iv];
1930  w2[ivc + id2 + iv] = w1[ivc + id2 + iv];
1931  w2[ivc + id3 + iv] = -w1[ivc + id3 + iv];
1932  w2[ivc + id4 + iv] = -w1[ivc + id4 + iv];
1933  }
1934  }
1935  }
1936  }
1937 
1938 }
1939 
1940 //====================================================================
1942  double *v1)
1943 {
1944  int Nvcd = m_Nvc * m_Nd;
1945  int Nxy = m_Nx2 * m_Ny;
1946 
1947  int id1 = 0;
1948  int id2 = m_Nvc;
1949  int id3 = m_Nvc * 2;
1950  int id4 = m_Nvc * 3;
1951 
1952  int isite = m_arg[itask].isite;
1953  double *w1 = &v1[Nvcd * isite];
1954 
1955  for (int it = 0; it < m_Mt; ++it) {
1956  for (int iz = 0; iz < m_Mz; ++iz) {
1957  for (int ixy = 0; ixy < Nxy; ++ixy) {
1958  int iv = Nvcd * (ixy + Nxy * (iz + m_Nz * it));
1959  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
1960  double wt1 = w1[ivc + id1 + iv];
1961  double wt2 = w1[ivc + id2 + iv];
1962  w1[ivc + id1 + iv] = w1[ivc + id3 + iv];
1963  w1[ivc + id2 + iv] = w1[ivc + id4 + iv];
1964  w1[ivc + id3 + iv] = wt1;
1965  w1[ivc + id4 + iv] = wt2;
1966  }
1967  }
1968  }
1969  }
1970 
1971 }
1972 
1973 //====================================================================
1975  double *v1)
1976 {
1977  int Nvcd = m_Nvc * m_Nd;
1978  int Nxy = m_Nx2 * m_Ny;
1979 
1980  // int id1 = 0;
1981  // int id2 = m_Nvc;
1982  int id3 = m_Nvc * 2;
1983  int id4 = m_Nvc * 3;
1984 
1985  int isite = m_arg[itask].isite;
1986  double *w1 = &v1[Nvcd * isite];
1987 
1988  for (int it = 0; it < m_Mt; ++it) {
1989  for (int iz = 0; iz < m_Mz; ++iz) {
1990  for (int ixy = 0; ixy < Nxy; ++ixy) {
1991  int iv = Nvcd * (ixy + Nxy * (iz + m_Nz * it));
1992  for (int ivc = 0; ivc < m_Nvc; ++ivc) {
1993  w1[ivc + id3 + iv] = -w1[ivc + id3 + iv];
1994  w1[ivc + id4 + iv] = -w1[ivc + id4 + iv];
1995  }
1996  }
1997  }
1998  }
1999 
2000 }
2001 
2002 //====================================================================
2003 //============================================================END=====
void mult_tm2_dirac_thread(int, double *, double *, int)
const Field_F Meo(const Field_F &, const int ieo)
void mult_tp1_dirac_thread(int, double *, double *, int)
BridgeIO vout
Definition: bridgeIO.cpp:207
void mult_ym1_thread(int, double *, double *, int)
void mult_xp1_thread(int, double *, double *, int)
void general(const char *format,...)
Definition: bridgeIO.cpp:38
void mult_tm2_chiral_thread(int, double *, double *, int)
void mult_ymb_thread(int, double *, double *, int)
void mult_tmb_dirac_thread(int, double *, double *, int)
void mult_zp1_thread(int, double *, double *, int)
void mult_xp2_thread(int, double *, double *, int)
void mult_ypb_thread(int, double *, double *, int)
SU(N) gauge field.
Definition: field_G.h:36
void mult_tm1_dirac_thread(int, double *, double *, int)
void mult_zm2_thread(int, double *, double *, int)
void mult_tmb_chiral_thread(int, double *, double *, int)
void mult_xm2_thread(int, double *, double *, int)
void mult_ym2_thread(int, double *, double *, int)
static int get_num_threads_available()
returns number of threads (works outside of parallel region).
void mult_xm1_thread(int, double *, double *, int)
void mult_tp1_chiral_thread(int, double *, double *, int)
void crucial(const char *format,...)
Definition: bridgeIO.cpp:26
void mult_tp2_dirac_thread(int, double *, double *, int)
void mult_tp2_chiral_thread(int, double *, double *, int)
void mult_tpb_chiral_thread(int, double *, double *, int)
void mult_zp2_thread(int, double *, double *, int)
void mult_xmb_thread(int, double *, double *, int)
void mult_zmb_thread(int, double *, double *, int)
void mult_yp1_thread(int, double *, double *, int)
void mult_zm1_thread(int, double *, double *, int)
void mult_zpb_thread(int, double *, double *, int)
void mult_yp2_thread(int, double *, double *, int)
void mult_xpb_thread(int, double *, double *, int)
void mult_tm1_chiral_thread(int, double *, double *, int)
void mult_tpb_dirac_thread(int, double *, double *, int)