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