Bridge++  Version 1.6.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
contract_4spinor.cpp
Go to the documentation of this file.
1 
14 #include "Measurements/Fermion/contract_4spinor.h"
15 #include "Field/index_lex.h"
16 
17 #include <cassert>
18 
19 #define NDIM 4
20 
21 #if defined USE_GROUP_SU3
22 #define NC 3
23 #define NC2 6
24 #define ND 4
25 #define NCD2 24
26 
27 //- NB. definition of C1,C2,C3 differs from those in Imp/
28 #define C1 0
29 #define C2 1
30 #define C3 2
31 #elif defined USE_GROUP_SU2
32 #define NC 2
33 #define NC2 4
34 #define ND 4
35 #define NCD2 16
36 #endif
37 
38 //====================================================================
39 void contract_at_t(dcomplex& corr,
40  const GammaMatrix& gm_sink,
41  const Field_F& v1, const Field_F& v2,
42  const int time)
43 {
44  const int Nc = CommonParameters::Nc();
45  const int Nd = CommonParameters::Nd();
46  const int Nx = CommonParameters::Nx();
47  const int Ny = CommonParameters::Ny();
48  const int Nz = CommonParameters::Nz();
49  const int Nt = CommonParameters::Nt();
50  const int Nvol = CommonParameters::Nvol();
51 
52  assert(Nvol == v1.nvol());
53  assert(Nvol == v2.nvol());
54  assert(time < Nt);
55 
56  Index_lex index;
57  std::vector<int> gm_index(Nd);
58  std::vector<double> corr_r(Nd), corr_i(Nd);
59 
60  for (int i = 0; i < Nd; ++i) {
61  gm_index[i] = gm_sink.index(i);
62  corr_r[i] = 0.0;
63  corr_i[i] = 0.0;
64  }
65 
66  for (int z = 0; z < Nz; ++z) {
67  for (int y = 0; y < Ny; ++y) {
68  for (int x = 0; x < Nx; ++x) {
69  int site = index.site(x, y, z, time);
70 
71  for (int s0 = 0; s0 < Nd; ++s0) {
72  int s1 = gm_index[s0];
73 
74  for (int c1 = 0; c1 < Nc; ++c1) {
75  corr_r[s0] += v1.cmp_r(c1, s1, site) * v2.cmp_r(c1, s0, site)
76  + v1.cmp_i(c1, s1, site) * v2.cmp_i(c1, s0, site);
77 
78  corr_i[s0] += -v1.cmp_r(c1, s1, site) * v2.cmp_i(c1, s0, site)
79  + v1.cmp_i(c1, s1, site) * v2.cmp_r(c1, s0, site);
80  }
81  }
82  }
83  }
84  }
85 
86  corr = cmplx(0.0, 0.0);
87  for (int s0 = 0; s0 < Nd; ++s0) {
88  corr += gm_sink.value(s0) * cmplx(corr_r[s0], corr_i[s0]);
89  }
90 }
91 
92 
93 //====================================================================
94 void contract_at_t(std::vector<dcomplex>& corr_global,
95  const GammaMatrix& gm_sink,
96  const Field_F& v1, const Field_F& v2)
97 {
98  const int Lt = CommonParameters::Lt();
99  const int Nt = CommonParameters::Nt();
100 
101  assert(corr_global.size() == Lt);
102 
103  std::vector<dcomplex> corr_local(Nt, 0.0);
104  for (int t = 0; t < Nt; ++t) {
105  dcomplex corr_t;
106  contract_at_t(corr_t, gm_sink, v1, v2, t);
107  corr_local[t] += corr_t;
108  }
109  global_corr_t(corr_global, corr_local);
110 }
111 
112 
113 //====================================================================
114 void contract_at_t(dcomplex& corr,
115  const std::vector<int>& momentum_sink,
116  const GammaMatrix& gm_sink,
117  const std::vector<int>& source_position,
118  const Field_F& v1, const Field_F& v2,
119  const int time)
120 {
121  const int Nc = CommonParameters::Nc();
122  const int Nd = CommonParameters::Nd();
123  const int Nx = CommonParameters::Nx();
124  const int Ny = CommonParameters::Ny();
125  const int Nz = CommonParameters::Nz();
126  const int Nt = CommonParameters::Nt();
127  const int Nvol = CommonParameters::Nvol();
128  const int Ndim = CommonParameters::Ndim();
129  const int Lx = CommonParameters::Lx();
130  const int Ly = CommonParameters::Ly();
131  const int Lz = CommonParameters::Lz();
132 
133  assert(Nvol == v1.nvol());
134  assert(Nvol == v2.nvol());
135  assert(time < Nt);
136  assert(momentum_sink.size() == Ndim - 1);
137 
138  Index_lex index;
139  std::vector<int> gm_index(Nd);
140  std::vector<double> corr_r(Nd), corr_i(Nd);
141 
142  for (int i = 0; i < Nd; ++i) {
143  gm_index[i] = gm_sink.index(i);
144  corr_r[i] = 0.0;
145  corr_i[i] = 0.0;
146  }
147 
148  static const double PI = 4.0 * atan(1.0);
149  std::vector<double> p_unit(Ndim - 1);
150  p_unit[0] = (2.0 * PI / Lx) * momentum_sink[0];
151  p_unit[1] = (2.0 * PI / Ly) * momentum_sink[1];
152  p_unit[2] = (2.0 * PI / Lz) * momentum_sink[2];
153 
154  std::vector<int> ipe(Ndim - 1);
155  ipe[0] = Communicator::ipe(0);
156  ipe[1] = Communicator::ipe(1);
157  ipe[2] = Communicator::ipe(2);
158 
159  for (int z = 0; z < Nz; ++z) {
160  for (int y = 0; y < Ny; ++y) {
161  for (int x = 0; x < Nx; ++x) {
162  int site = index.site(x, y, z, time);
163 
164  int x_global = x + ipe[0] * Nx;
165  int y_global = y + ipe[1] * Ny;
166  int z_global = z + ipe[2] * Nz;
167 
168  double p_x = p_unit[0] * (x_global - source_position[0]);
169  double p_y = p_unit[1] * (y_global - source_position[1]);
170  double p_z = p_unit[2] * (z_global - source_position[2]);
171 
172  double cos_p_xyz = cos(p_x + p_y + p_z);
173  double sin_p_xyz = sin(p_x + p_y + p_z);
174 
175  for (int s0 = 0; s0 < Nd; ++s0) {
176  int s1 = gm_index[s0];
177 
178  double v1_v2_r = 0.0;
179  double v1_v2_i = 0.0;
180 
181  for (int c1 = 0; c1 < Nc; ++c1) {
182  v1_v2_r += v1.cmp_r(c1, s1, site) * v2.cmp_r(c1, s0, site)
183  + v1.cmp_i(c1, s1, site) * v2.cmp_i(c1, s0, site);
184 
185  v1_v2_i += -v1.cmp_r(c1, s1, site) * v2.cmp_i(c1, s0, site)
186  + v1.cmp_i(c1, s1, site) * v2.cmp_r(c1, s0, site);
187  }
188 
189  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
190  corr_r[s0] += v1_v2_r * cos_p_xyz - v1_v2_i * sin_p_xyz;
191  corr_i[s0] += v1_v2_r * sin_p_xyz + v1_v2_i * cos_p_xyz;
192  }
193  }
194  }
195  }
196 
197  corr = cmplx(0.0, 0.0);
198  for (int s0 = 0; s0 < Nd; ++s0) {
199  corr += gm_sink.value(s0) * cmplx(corr_r[s0], corr_i[s0]);
200  }
201 }
202 
203 
204 //====================================================================
205 void contract_at_t(std::vector<dcomplex>& corr_global,
206  const std::vector<int>& momentum_sink,
207  const GammaMatrix& gm_sink,
208  const std::vector<int>& source_position,
209  const Field_F& v1, const Field_F& v2)
210 {
211  const int Lt = CommonParameters::Lt();
212  const int Nt = CommonParameters::Nt();
213 
214  assert(corr_global.size() == Lt);
215 
216  std::vector<dcomplex> corr_local(Nt, 0.0);
217  for (int t = 0; t < Nt; ++t) {
218  dcomplex corr_t;
219  contract_at_t(corr_t, momentum_sink, gm_sink, source_position, v1, v2, t);
220  corr_local[t] += corr_t;
221  }
222  global_corr_t(corr_global, corr_local);
223 }
224 
225 
226 //====================================================================
227 void contract_at_t_cos(dcomplex& corr,
228  const std::vector<int>& momentum_sink,
229  const GammaMatrix& gm_sink,
230  const std::vector<int>& source_position,
231  const Field_F& v1, const Field_F& v2,
232  const int time)
233 {
234 #if defined USE_GROUP_SU_N
235  const int NC = CommonParameters::Nc();
236  const int ND = CommonParameters::Nd();
237  const int NC2 = 2 * NC;
238  const int NCD2 = NC2 * ND;
239 #endif
240 
241  const int Nvol = v1.nvol();
242  const int Nvol_s = Nvol / CommonParameters::Nt();
243  const int Nx = CommonParameters::Nx();
244  const int Ny = CommonParameters::Ny();
245  const int Nz = CommonParameters::Nz();
246  const int Lx = CommonParameters::Lx();
247  const int Ly = CommonParameters::Ly();
248  const int Lz = CommonParameters::Lz();
249 
250  assert(Nvol == v2.nvol());
251  assert(v1.nex() == 1);
252  assert(v2.nex() == 1);
253  assert(momentum_sink.size() == NDIM - 1);
254 
255  const double *w1 = v1.ptr(0);
256  const double *w2 = v2.ptr(0);
257 
258  int id1[ND];
259  int id2[ND];
260  for (int id = 0; id < ND; ++id) {
261  id1[id] = id * NC2;
262  id2[id] = gm_sink.index(id) * NC2;
263  }
264 
265  double c_r[ND];
266  double c_i[ND];
267  for (int id = 0; id < ND; ++id) {
268  c_r[id] = 0.0;
269  c_i[id] = 0.0;
270  }
271 
272  static const double PI = 4.0 * atan(1.0);
273  std::vector<double> p_unit(ND - 1);
274  p_unit[0] = (2.0 * PI / Lx) * momentum_sink[0];
275  p_unit[1] = (2.0 * PI / Ly) * momentum_sink[1];
276  p_unit[2] = (2.0 * PI / Lz) * momentum_sink[2];
277 
278  std::vector<int> ipe(ND - 1);
279  ipe[0] = Communicator::ipe(0);
280  ipe[1] = Communicator::ipe(1);
281  ipe[2] = Communicator::ipe(2);
282 
283 
284  for (int ss = 0; ss < Nvol_s; ++ss) {
285  int site = NCD2 * (ss + time * Nvol_s);
286 
287  int x = ss % Nx;
288  int y = ss % (Nx * Ny) / Nx;
289  int z = ss % (Nx * Ny * Nz) / (Nx * Ny);
290 
291  int x_global = x + ipe[0] * Nx;
292  int y_global = y + ipe[1] * Ny;
293  int z_global = z + ipe[2] * Nz;
294 
295  double p_x = p_unit[0] * (x_global - source_position[0]);
296  double p_y = p_unit[1] * (y_global - source_position[1]);
297  double p_z = p_unit[2] * (z_global - source_position[2]);
298 
299  double cos_p_xyz = cos(p_x + p_y + p_z);
300  //double sin_p_xyz = sin(p_x + p_y + p_z);
301  double sin_p_xyz = 0;
302 
303  for (int cc = 0; cc < NC; ++cc) {
304  for (int id = 0; id < ND; ++id) {
305  int ic1_r = 2 * cc + id1[id] + site;
306  int ic2_r = 2 * cc + id2[id] + site;
307 
308  int ic1_i = 2 * cc + 1 + id1[id] + site;
309  int ic2_i = 2 * cc + 1 + id2[id] + site;
310 
311  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
312  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
313  //double w1_w2_i = 0;
314 
315  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
316  c_r[id] += w1_w2_r * cos_p_xyz - w1_w2_i * sin_p_xyz;
317  c_i[id] += w1_w2_r * sin_p_xyz + w1_w2_i * cos_p_xyz;
318  }
319  }
320  }
321 
322  corr = cmplx(0.0, 0.0);
323  for (int id = 0; id < ND; ++id) {
324  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
325  }
326 }
327 
328 
329 //====================================================================
330 void contract_at_t_cos(std::vector<dcomplex>& corr_global,
331  const std::vector<int>& momentum_sink,
332  const GammaMatrix& gm_sink,
333  const std::vector<int>& source_position,
334  const Field_F& v1, const Field_F& v2)
335 {
336  const int Lt = CommonParameters::Lt();
337  const int Nt = CommonParameters::Nt();
338 
339  assert(corr_global.size() == Lt);
340 
341  std::vector<dcomplex> corr_local(Nt, 0.0);
342  for (int t = 0; t < Nt; ++t) {
343  dcomplex corr_t;
344  contract_at_t_cos(corr_t, momentum_sink, gm_sink, source_position, v1, v2, t);
345  corr_local[t] += corr_t;
346  }
347  global_corr_t(corr_global, corr_local);
348 }
349 
350 
351 //====================================================================
352 void contract_at_t(dcomplex& corr,
353  const GammaMatrix& gm_sink,
354  const int i_alpha,
355  const Field_F& v1, const Field_F& v2, const Field_F& v3,
356  const int time)
357 {
358 #if defined USE_GROUP_SU3
359  const int Nc = CommonParameters::Nc();
360  const int Nd = CommonParameters::Nd();
361  const int Nx = CommonParameters::Nx();
362  const int Ny = CommonParameters::Ny();
363  const int Nz = CommonParameters::Nz();
364  const int Nt = CommonParameters::Nt();
365  const int Nvol = CommonParameters::Nvol();
366 
367  assert(Nc == 3);
368  assert(Nvol == v1.nvol());
369  assert(Nvol == v2.nvol());
370  assert(Nvol == v3.nvol());
371  assert(time < Nt);
372 
373  Index_lex index;
374  std::vector<int> gm_index(Nd);
375  std::vector<double> c_r(Nd), c_i(Nd);
376 
377  for (int i = 0; i < Nd; ++i) {
378  gm_index[i] = gm_sink.index(i);
379  c_r[i] = 0.0;
380  c_i[i] = 0.0;
381  }
382 
383  for (int z = 0; z < Nz; ++z) {
384  for (int y = 0; y < Ny; ++y) {
385  for (int x = 0; x < Nx; ++x) {
386  int site = index.site(x, y, z, time);
387 
388  for (int d1 = 0; d1 < Nd; ++d1) {
389  int d2 = gm_index[d1];
390  int d3 = i_alpha;
391 
392  c_r[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
393  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C3, d3, site)
394  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
395  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C3, d3, site);
396  c_i[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
397  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C3, d3, site)
398  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
399  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C3, d3, site);
400 
401  c_r[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
402  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C1, d3, site)
403  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
404  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C1, d3, site);
405  c_i[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
406  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C1, d3, site)
407  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
408  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C1, d3, site);
409 
410  c_r[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
411  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C2, d3, site)
412  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
413  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C2, d3, site);
414  c_i[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
415  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C2, d3, site)
416  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
417  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C2, d3, site);
418 
419  c_r[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
420  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C1, d3, site)
421  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
422  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C1, d3, site);
423  c_i[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
424  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C1, d3, site)
425  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
426  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C1, d3, site);
427 
428  c_r[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
429  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C3, d3, site)
430  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
431  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C3, d3, site);
432  c_i[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
433  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C3, d3, site)
434  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
435  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C3, d3, site);
436 
437  c_r[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
438  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C2, d3, site)
439  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
440  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C2, d3, site);
441  c_i[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
442  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C2, d3, site)
443  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
444  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C2, d3, site);
445  }
446  }
447  }
448  }
449 
450  corr = cmplx(0.0, 0.0);
451  for (int s0 = 0; s0 < Nd; ++s0) {
452  corr += gm_sink.value(s0) * cmplx(c_r[s0], c_i[s0]);
453  }
454 #endif // (USE_GROUP_SU3)
455 }
456 
457 
458 //====================================================================
459 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta at x
460 void contract_at_x(dcomplex& corr, const GammaMatrix& gm_sink,
461  const Field_F& v1, const Field_F& v2,
462  int x)
463 {
464 #if defined USE_GROUP_SU_N
465  const int NC = CommonParameters::Nc();
466  const int ND = CommonParameters::Nd();
467  const int NC2 = 2 * NC;
468  const int NCD2 = NC2 * ND;
469 #endif
470 
471  const int Nvol = v1.nvol();
472  const int Nvol_s = Nvol / CommonParameters::Nt();
473  const int Nx = CommonParameters::Nx();
474  const int Ny = CommonParameters::Ny();
475  const int Nz = CommonParameters::Nz();
476  const int Nt = CommonParameters::Nt();
477 
478  assert(Nvol == v2.nvol());
479  assert(v1.nex() == 1);
480  assert(v2.nex() == 1);
481 
482  const double *w1 = v1.ptr(0);
483  const double *w2 = v2.ptr(0);
484 
485  int id1[ND];
486  int id2[ND];
487  for (int id = 0; id < ND; ++id) {
488  id1[id] = id * NC2;
489  id2[id] = gm_sink.index(id) * NC2;
490  }
491 
492  double c_r[ND];
493  double c_i[ND];
494  for (int id = 0; id < ND; ++id) {
495  c_r[id] = 0.0;
496  c_i[id] = 0.0;
497  }
498 
499  for (int t = 0; t < Nt; ++t) {
500  for (int z = 0; z < Nz; ++z) {
501  for (int y = 0; y < Ny; ++y) {
502  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
503 
504  for (int cc = 0; cc < NC; ++cc) {
505  for (int id = 0; id < ND; ++id) {
506  int ic1_r = 2 * cc + id1[id] + site;
507  int ic2_r = 2 * cc + id2[id] + site;
508 
509  int ic1_i = 2 * cc + 1 + id1[id] + site;
510  int ic2_i = 2 * cc + 1 + id2[id] + site;
511 
512  c_r[id] += w1[ic2_r] * w2[ic1_r]
513  + w1[ic2_i] * w2[ic1_i];
514 
515  c_i[id] += -w1[ic2_r] * w2[ic1_i]
516  + w1[ic2_i] * w2[ic1_r];
517  }
518  }
519  }
520  }
521  }
522 
523  corr = cmplx(0.0, 0.0);
524  for (int id = 0; id < ND; ++id) {
525  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
526  }
527 }
528 
529 
530 //====================================================================
531 void contract_at_x(std::vector<dcomplex>& corr_global,
532  const GammaMatrix& gm_sink,
533  const Field_F& v1, const Field_F& v2)
534 {
535  const int Lx = CommonParameters::Lx();
536  const int Nx = CommonParameters::Nx();
537 
538  assert(corr_global.size() == Lx);
539 
540  std::vector<dcomplex> corr_local(Nx, 0.0);
541  for (int x = 0; x < Nx; ++x) {
542  dcomplex corr_x;
543  contract_at_x(corr_x, gm_sink, v1, v2, x);
544  corr_local[x] += corr_x;
545  }
546  global_corr_x(corr_global, corr_local);
547 }
548 
549 
550 //====================================================================
551 void contract_at_x(dcomplex& corr,
552  const std::vector<int>& momentum_sink,
553  const GammaMatrix& gm_sink,
554  const std::vector<int>& source_position,
555  const Field_F& v1, const Field_F& v2,
556  const int x)
557 {
558 #if defined USE_GROUP_SU_N
559  int NC = CommonParameters::Nc();
560  int ND = CommonParameters::Nd();
561  int NC2 = 2 * NC;
562  int NCD2 = NC2 * ND;
563 #endif
564 
565  const int Nvol = v1.nvol();
566  const int Nvol_s = Nvol / CommonParameters::Nt();
567  const int Nx = CommonParameters::Nx();
568  const int Ny = CommonParameters::Ny();
569  const int Nz = CommonParameters::Nz();
570  const int Nt = CommonParameters::Nt();
571  const int Lx = CommonParameters::Lx();
572  const int Ly = CommonParameters::Ly();
573  const int Lz = CommonParameters::Lz();
574  const int Lt = CommonParameters::Lt();
575 
576  assert(Nvol == v2.nvol());
577  assert(v1.nex() == 1);
578  assert(v2.nex() == 1);
579  assert(momentum_sink.size() == NDIM - 1);
580  assert(source_position.size() == NDIM);
581 
582  const double *w1 = v1.ptr(0);
583  const double *w2 = v2.ptr(0);
584 
585  int id1[ND];
586  int id2[ND];
587  for (int id = 0; id < ND; ++id) {
588  id1[id] = id * NC2;
589  id2[id] = gm_sink.index(id) * NC2;
590  }
591 
592  double c_r[ND];
593  double c_i[ND];
594  for (int id = 0; id < ND; ++id) {
595  c_r[id] = 0.0;
596  c_i[id] = 0.0;
597  }
598 
599  static const double PI = 4.0 * atan(1.0);
600  std::vector<double> p_unit(ND - 1);
601  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
602  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
603  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
604 
605  std::vector<int> ipe(ND);
606  ipe[0] = Communicator::ipe(0);
607  ipe[1] = Communicator::ipe(1);
608  ipe[2] = Communicator::ipe(2);
609  ipe[3] = Communicator::ipe(3);
610 
611  for (int t = 0; t < Nt; ++t) {
612  for (int z = 0; z < Nz; ++z) {
613  for (int y = 0; y < Ny; ++y) {
614  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
615 
616  int y_global = y + ipe[1] * Ny;
617  int z_global = z + ipe[2] * Nz;
618  int t_global = t + ipe[3] * Nt;
619 
620  double p_y = p_unit[0] * (y_global - source_position[1]);
621  double p_z = p_unit[1] * (z_global - source_position[2]);
622  double p_t = p_unit[2] * (t_global - source_position[3]);
623 
624  double cos_p_yzt = cos(p_t + p_y + p_z);
625  double sin_p_yzt = sin(p_t + p_y + p_z);
626 
627  for (int cc = 0; cc < NC; ++cc) {
628  for (int id = 0; id < ND; ++id) {
629  int ic1_r = 2 * cc + id1[id] + site;
630  int ic2_r = 2 * cc + id2[id] + site;
631 
632  int ic1_i = 2 * cc + 1 + id1[id] + site;
633  int ic2_i = 2 * cc + 1 + id2[id] + site;
634 
635  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
636  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
637 
638  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
639  c_r[id] += w1_w2_r * cos_p_yzt - w1_w2_i * sin_p_yzt;
640  c_i[id] += w1_w2_r * sin_p_yzt + w1_w2_i * cos_p_yzt;
641  }
642  }
643  }
644  }
645  }
646 
647  corr = cmplx(0.0, 0.0);
648  for (int id = 0; id < ND; ++id) {
649  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
650  }
651 }
652 
653 
654 //====================================================================
655 void contract_at_x(std::vector<dcomplex>& corr_global,
656  const std::vector<int>& momentum_sink,
657  const GammaMatrix& gm_sink,
658  const std::vector<int>& source_position,
659  const Field_F& v1, const Field_F& v2)
660 {
661  const int Lx = CommonParameters::Lx();
662  const int Nx = CommonParameters::Nx();
663 
664  assert(corr_global.size() == Lx);
665 
666  std::vector<dcomplex> corr_local(Nx, 0.0);
667  for (int x = 0; x < Nx; ++x) {
668  dcomplex corr_x;
669  contract_at_x(corr_x, momentum_sink, gm_sink, source_position, v1, v2, x);
670  corr_local[x] += corr_x;
671  }
672  global_corr_x(corr_global, corr_local);
673 }
674 
675 
676 //====================================================================
677 void contract_at_x_cos(dcomplex& corr,
678  const std::vector<int>& momentum_sink,
679  const GammaMatrix& gm_sink,
680  const std::vector<int>& source_position,
681  const Field_F& v1, const Field_F& v2,
682  const int x)
683 {
684 #if defined USE_GROUP_SU_N
685  int NC = CommonParameters::Nc();
686  int ND = CommonParameters::Nd();
687  int NC2 = 2 * NC;
688  int NCD2 = NC2 * ND;
689 #endif
690 
691  const int Nvol = v1.nvol();
692  const int Nvol_s = Nvol / CommonParameters::Nt();
693  const int Nx = CommonParameters::Nx();
694  const int Ny = CommonParameters::Ny();
695  const int Nz = CommonParameters::Nz();
696  const int Nt = CommonParameters::Nt();
697  const int Lx = CommonParameters::Lx();
698  const int Ly = CommonParameters::Ly();
699  const int Lz = CommonParameters::Lz();
700  const int Lt = CommonParameters::Lt();
701 
702  assert(Nvol == v2.nvol());
703  assert(v1.nex() == 1);
704  assert(v2.nex() == 1);
705  assert(momentum_sink.size() == NDIM - 1);
706  assert(source_position.size() == NDIM);
707 
708  const double *w1 = v1.ptr(0);
709  const double *w2 = v2.ptr(0);
710 
711  int id1[ND];
712  int id2[ND];
713  for (int id = 0; id < ND; ++id) {
714  id1[id] = id * NC2;
715  id2[id] = gm_sink.index(id) * NC2;
716  }
717 
718  double c_r[ND];
719  double c_i[ND];
720  for (int id = 0; id < ND; ++id) {
721  c_r[id] = 0.0;
722  c_i[id] = 0.0;
723  }
724 
725  static const double PI = 4.0 * atan(1.0);
726  std::vector<double> p_unit(ND - 1);
727  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
728  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
729  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
730 
731  std::vector<int> ipe(ND);
732  ipe[0] = Communicator::ipe(0);
733  ipe[1] = Communicator::ipe(1);
734  ipe[2] = Communicator::ipe(2);
735  ipe[3] = Communicator::ipe(3);
736 
737 
738  for (int t = 0; t < Nt; ++t) {
739  for (int z = 0; z < Nz; ++z) {
740  for (int y = 0; y < Ny; ++y) {
741  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
742 
743  int y_global = y + ipe[1] * Ny;
744  int z_global = z + ipe[2] * Nz;
745  int t_global = t + ipe[3] * Nt;
746 
747  double p_y = p_unit[0] * (y_global - source_position[1]);
748  double p_z = p_unit[1] * (z_global - source_position[2]);
749  double p_t = p_unit[2] * (t_global - source_position[3]);
750 
751  double cos_p_yzt = cos(p_t + p_y + p_z);
752  //double sin_p_yzt = sin(p_t + p_y + p_z);
753  double sin_p_yzt = 0;
754 
755  for (int cc = 0; cc < NC; ++cc) {
756  for (int id = 0; id < ND; ++id) {
757  int ic1_r = 2 * cc + id1[id] + site;
758  int ic2_r = 2 * cc + id2[id] + site;
759 
760  int ic1_i = 2 * cc + 1 + id1[id] + site;
761  int ic2_i = 2 * cc + 1 + id2[id] + site;
762 
763  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
764  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
765  //double w1_w2_i = 0;
766 
767  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
768  c_r[id] += w1_w2_r * cos_p_yzt - w1_w2_i * sin_p_yzt;
769  c_i[id] += w1_w2_r * sin_p_yzt + w1_w2_i * cos_p_yzt;
770  }
771  }
772  }
773  }
774  }
775 
776  corr = cmplx(0.0, 0.0);
777  for (int id = 0; id < ND; ++id) {
778  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
779  }
780 }
781 
782 
783 //====================================================================
784 void contract_at_x_cos(std::vector<dcomplex>& corr_global,
785  const std::vector<int>& momentum_sink,
786  const GammaMatrix& gm_sink,
787  const std::vector<int>& source_position,
788  const Field_F& v1, const Field_F& v2)
789 {
790  const int Lx = CommonParameters::Lx();
791  const int Nx = CommonParameters::Nx();
792 
793  assert(corr_global.size() == Lx);
794 
795  std::vector<dcomplex> corr_local(Nx, 0.0);
796  for (int x = 0; x < Nx; ++x) {
797  dcomplex corr_x;
798  contract_at_x_cos(corr_x, momentum_sink, gm_sink, source_position, v1, v2, x);
799  corr_local[x] += corr_x;
800  }
801  global_corr_x(corr_global, corr_local);
802 }
803 
804 
805 //====================================================================
806 void contract_at_x(dcomplex& corr,
807  const GammaMatrix& gm_sink,
808  const int i_alpha,
809  const Field_F& v1, const Field_F& v2, const Field_F& v3,
810  const int x)
811 {
812 #if defined USE_GROUP_SU3
813  const int Nvol = v1.nvol();
814  const int Nvol_s = Nvol / CommonParameters::Nt();
815  const int Nx = CommonParameters::Nx();
816  const int Ny = CommonParameters::Ny();
817  const int Nz = CommonParameters::Nz();
818  const int Nt = CommonParameters::Nt();
819 
820  assert(Nvol == v2.nvol());
821  assert(Nvol == v3.nvol());
822  assert(v1.nex() == 1);
823  assert(v2.nex() == 1);
824  assert(v3.nex() == 1);
825 
826  const double *w1 = v1.ptr(0);
827  const double *w2 = v2.ptr(0);
828  const double *w3 = v3.ptr(0);
829 
830  int id1[ND];
831  int id2[ND];
832  for (int id = 0; id < ND; ++id) {
833  id1[id] = id * NC2;
834  id2[id] = gm_sink.index(id) * NC2;
835  }
836  int id3 = i_alpha * NC2;
837 
838  double c_r[ND];
839  double c_i[ND];
840  for (int id = 0; id < ND; ++id) {
841  c_r[id] = 0.0;
842  c_i[id] = 0.0;
843  }
844 
845 
846  for (int t = 0; t < Nt; ++t) {
847  for (int z = 0; z < Nz; ++z) {
848  for (int y = 0; y < Ny; ++y) {
849  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
850 
851  for (int id = 0; id < ND; ++id) {
852  int ic11_r = C1 + id1[id] + site;
853  int ic22_r = C2 + id2[id] + site;
854  int ic33_r = C3 + id3 + site;
855 
856  int ic11_i = C1 + 1 + id1[id] + site;
857  int ic22_i = C2 + 1 + id2[id] + site;
858  int ic33_i = C3 + 1 + id3 + site;
859 
860  int ic21_r = C2 + id1[id] + site;
861  int ic32_r = C3 + id2[id] + site;
862  int ic13_r = C1 + id3 + site;
863 
864  int ic21_i = C2 + 1 + id1[id] + site;
865  int ic32_i = C3 + 1 + id2[id] + site;
866  int ic13_i = C1 + 1 + id3 + site;
867 
868  int ic31_r = C3 + id1[id] + site;
869  int ic12_r = C1 + id2[id] + site;
870  int ic23_r = C2 + id3 + site;
871 
872  int ic31_i = C3 + 1 + id1[id] + site;
873  int ic12_i = C1 + 1 + id2[id] + site;
874  int ic23_i = C2 + 1 + id3 + site;
875 
876  c_r[id] += (w1[ic11_r] * w2[ic22_r] - w1[ic11_i] * w2[ic22_i]) * w3[ic33_r]
877  - (w1[ic11_r] * w2[ic22_i] + w1[ic11_i] * w2[ic22_r]) * w3[ic33_i];
878  c_i[id] += (w1[ic11_r] * w2[ic22_r] - w1[ic11_i] * w2[ic22_i]) * w3[ic33_i]
879  + (w1[ic11_r] * w2[ic22_i] + w1[ic11_i] * w2[ic22_r]) * w3[ic33_r];
880 
881  c_r[id] += (w1[ic21_r] * w2[ic32_r] - w1[ic21_i] * w2[ic32_i]) * w3[ic13_r]
882  - (w1[ic21_r] * w2[ic32_i] + w1[ic21_i] * w2[ic32_r]) * w3[ic13_i];
883  c_i[id] += (w1[ic21_r] * w2[ic32_r] - w1[ic21_i] * w2[ic32_i]) * w3[ic13_i]
884  + (w1[ic21_r] * w2[ic32_i] + w1[ic21_i] * w2[ic32_r]) * w3[ic13_r];
885 
886  c_r[id] += (w1[ic31_r] * w2[ic12_r] - w1[ic31_i] * w2[ic12_i]) * w3[ic23_r]
887  - (w1[ic31_r] * w2[ic12_i] + w1[ic31_i] * w2[ic12_r]) * w3[ic23_i];
888  c_i[id] += (w1[ic31_r] * w2[ic12_r] - w1[ic31_i] * w2[ic12_i]) * w3[ic23_i]
889  + (w1[ic31_r] * w2[ic12_i] + w1[ic31_i] * w2[ic12_r]) * w3[ic23_r];
890 
891  c_r[id] -= (w1[ic31_r] * w2[ic22_r] - w1[ic31_i] * w2[ic22_i]) * w3[ic13_r]
892  - (w1[ic31_r] * w2[ic22_i] + w1[ic31_i] * w2[ic22_r]) * w3[ic13_i];
893  c_i[id] -= (w1[ic31_r] * w2[ic22_r] - w1[ic31_i] * w2[ic22_i]) * w3[ic13_i]
894  + (w1[ic31_r] * w2[ic22_i] + w1[ic31_i] * w2[ic22_r]) * w3[ic13_r];
895 
896  c_r[id] -= (w1[ic21_r] * w2[ic12_r] - w1[ic21_i] * w2[ic12_i]) * w3[ic33_r]
897  - (w1[ic21_r] * w2[ic12_i] + w1[ic21_i] * w2[ic12_r]) * w3[ic33_i];
898  c_i[id] -= (w1[ic21_r] * w2[ic12_r] - w1[ic21_i] * w2[ic12_i]) * w3[ic33_i]
899  + (w1[ic21_r] * w2[ic12_i] + w1[ic21_i] * w2[ic12_r]) * w3[ic33_r];
900 
901  c_r[id] -= (w1[ic11_r] * w2[ic32_r] - w1[ic11_i] * w2[ic32_i]) * w3[ic23_r]
902  - (w1[ic11_r] * w2[ic32_i] + w1[ic11_i] * w2[ic32_r]) * w3[ic23_i];
903  c_i[id] -= (w1[ic11_r] * w2[ic32_r] - w1[ic11_i] * w2[ic32_i]) * w3[ic23_i]
904  + (w1[ic11_r] * w2[ic32_i] + w1[ic11_i] * w2[ic32_r]) * w3[ic23_r];
905  }
906  }
907  }
908  }
909 
910  corr = cmplx(0.0, 0.0);
911  for (int id = 0; id < ND; ++id) {
912  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
913  }
914 #endif // (USE_GROUP_SU3)
915 }
916 
917 
918 //====================================================================
919 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta
920 void contract_at_y(dcomplex& corr,
921  const GammaMatrix& gm_sink,
922  const Field_F& v1, const Field_F& v2,
923  const int y)
924 {
925 #if defined USE_GROUP_SU_N
926  const int NC = CommonParameters::Nc();
927  const int ND = CommonParameters::Nd();
928  const int NC2 = 2 * NC;
929  const int NCD2 = NC2 * ND;
930 #endif
931 
932  const int Nvol = v1.nvol();
933  int Nx = CommonParameters::Nx();
934  int Ny = CommonParameters::Ny();
935  int Nz = CommonParameters::Nz();
936  int Nt = CommonParameters::Nt();
937 
938  assert(Nvol == v1.nvol());
939  assert(Nvol == v2.nvol());
940  assert(v1.nex() == 1);
941  assert(v2.nex() == 1);
942  assert(y < Ny);
943 
944  const double *w1 = v1.ptr(0);
945  const double *w2 = v2.ptr(0);
946 
947  int id1[ND];
948  int id2[ND];
949  for (int id = 0; id < ND; ++id) {
950  id1[id] = id * NC2;
951  id2[id] = gm_sink.index(id) * NC2;
952  }
953 
954  double c_r[ND];
955  double c_i[ND];
956  for (int id = 0; id < ND; ++id) {
957  c_r[id] = 0.0;
958  c_i[id] = 0.0;
959  }
960 
961  for (int t = 0; t < Nt; ++t) {
962  for (int z = 0; z < Nz; ++z) {
963  for (int x = 0; x < Nx; ++x) {
964  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
965  for (int cc = 0; cc < NC; ++cc) {
966  for (int id = 0; id < ND; ++id) {
967  int ic1_r = 2 * cc + id1[id] + site;
968  int ic2_r = 2 * cc + id2[id] + site;
969 
970  int ic1_i = 2 * cc + 1 + id1[id] + site;
971  int ic2_i = 2 * cc + 1 + id2[id] + site;
972 
973  c_r[id] += w1[ic2_r] * w2[ic1_r]
974  + w1[ic2_i] * w2[ic1_i];
975 
976  c_i[id] += -w1[ic2_r] * w2[ic1_i]
977  + w1[ic2_i] * w2[ic1_r];
978  }
979  }
980  }
981  }
982  }
983 
984  corr = cmplx(0.0, 0.0);
985  for (int id = 0; id < ND; ++id) {
986  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
987  }
988 }
989 
990 
991 //====================================================================
992 void contract_at_y(std::vector<dcomplex>& corr_global,
993  const GammaMatrix& gm_sink,
994  const Field_F& v1, const Field_F& v2)
995 {
996  const int Ly = CommonParameters::Ly();
997  const int Ny = CommonParameters::Ny();
998 
999  assert(corr_global.size() == Ly);
1000 
1001  std::vector<dcomplex> corr_local(Ny, 0.0);
1002  for (int y = 0; y < Ny; ++y) {
1003  dcomplex corr_y;
1004  contract_at_y(corr_y, gm_sink, v1, v2, y);
1005  corr_local[y] += corr_y;
1006  }
1007  global_corr_y(corr_global, corr_local);
1008 }
1009 
1010 
1011 //====================================================================
1012 void contract_at_y(dcomplex& corr,
1013  const std::vector<int>& momentum_sink,
1014  const GammaMatrix& gm_sink,
1015  const std::vector<int>& source_position,
1016  const Field_F& v1, const Field_F& v2,
1017  const int y)
1018 {
1019 #if defined USE_GROUP_SU_N
1020  int NC = CommonParameters::Nc();
1021  int ND = CommonParameters::Nd();
1022  int NC2 = 2 * NC;
1023  int NCD2 = NC2 * ND;
1024 #endif
1025 
1026  const int Nvol = v1.nvol();
1027  const int Nvol_s = Nvol / CommonParameters::Nt();
1028  const int Nx = CommonParameters::Nx();
1029  const int Ny = CommonParameters::Ny();
1030  const int Nz = CommonParameters::Nz();
1031  const int Nt = CommonParameters::Nt();
1032  const int Lx = CommonParameters::Lx();
1033  const int Ly = CommonParameters::Ly();
1034  const int Lz = CommonParameters::Lz();
1035  const int Lt = CommonParameters::Lt();
1036 
1037  assert(Nvol == v2.nvol());
1038  assert(v1.nex() == 1);
1039  assert(v2.nex() == 1);
1040  assert(momentum_sink.size() == NDIM - 1);
1041  assert(source_position.size() == NDIM);
1042 
1043  const double *w1 = v1.ptr(0);
1044  const double *w2 = v2.ptr(0);
1045 
1046  int id1[ND];
1047  int id2[ND];
1048  for (int id = 0; id < ND; ++id) {
1049  id1[id] = id * NC2;
1050  id2[id] = gm_sink.index(id) * NC2;
1051  }
1052 
1053  double c_r[ND];
1054  double c_i[ND];
1055  for (int id = 0; id < ND; ++id) {
1056  c_r[id] = 0.0;
1057  c_i[id] = 0.0;
1058  }
1059 
1060  static const double PI = 4.0 * atan(1.0);
1061  std::vector<double> p_unit(ND - 1);
1062  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1063  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1064  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1065 
1066  std::vector<int> ipe(ND);
1067  ipe[0] = Communicator::ipe(0);
1068  ipe[1] = Communicator::ipe(1);
1069  ipe[2] = Communicator::ipe(2);
1070  ipe[3] = Communicator::ipe(3);
1071 
1072  for (int t = 0; t < Nt; ++t) {
1073  for (int z = 0; z < Nz; ++z) {
1074  for (int x = 0; x < Nx; ++x) {
1075  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1076 
1077  int x_global = x + ipe[0] * Nx;
1078  int z_global = z + ipe[2] * Nz;
1079  int t_global = t + ipe[3] * Nt;
1080 
1081  double p_x = p_unit[0] * (x_global - source_position[0]);
1082  double p_z = p_unit[1] * (z_global - source_position[2]);
1083  double p_t = p_unit[2] * (t_global - source_position[3]);
1084 
1085  double cos_p_xzt = cos(p_t + p_x + p_z);
1086  double sin_p_xzt = sin(p_t + p_x + p_z);
1087 
1088  for (int cc = 0; cc < NC; ++cc) {
1089  for (int id = 0; id < ND; ++id) {
1090  int ic1_r = 2 * cc + id1[id] + site;
1091  int ic2_r = 2 * cc + id2[id] + site;
1092 
1093  int ic1_i = 2 * cc + 1 + id1[id] + site;
1094  int ic2_i = 2 * cc + 1 + id2[id] + site;
1095 
1096  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1097  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1098 
1099  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1100  c_r[id] += w1_w2_r * cos_p_xzt - w1_w2_i * sin_p_xzt;
1101  c_i[id] += w1_w2_r * sin_p_xzt + w1_w2_i * cos_p_xzt;
1102  }
1103  }
1104  }
1105  }
1106  }
1107 
1108  corr = cmplx(0.0, 0.0);
1109  for (int id = 0; id < ND; ++id) {
1110  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1111  }
1112 }
1113 
1114 
1115 //====================================================================
1116 void contract_at_y(std::vector<dcomplex>& corr_global,
1117  const std::vector<int>& momentum_sink,
1118  const GammaMatrix& gm_sink,
1119  const std::vector<int>& source_position,
1120  const Field_F& v1, const Field_F& v2)
1121 {
1122  const int Ly = CommonParameters::Ly();
1123  const int Ny = CommonParameters::Ny();
1124 
1125  assert(corr_global.size() == Ly);
1126 
1127  std::vector<dcomplex> corr_local(Ny, 0.0);
1128  for (int y = 0; y < Ny; ++y) {
1129  dcomplex corr_y;
1130  contract_at_y(corr_y, momentum_sink, gm_sink, source_position, v1, v2, y);
1131  corr_local[y] += corr_y;
1132  }
1133  global_corr_y(corr_global, corr_local);
1134 }
1135 
1136 
1137 //====================================================================
1138 void contract_at_y_cos(dcomplex& corr,
1139  const std::vector<int>& momentum_sink,
1140  const GammaMatrix& gm_sink,
1141  const std::vector<int>& source_position,
1142  const Field_F& v1, const Field_F& v2,
1143  const int y)
1144 {
1145 #if defined USE_GROUP_SU_N
1146  int NC = CommonParameters::Nc();
1147  int ND = CommonParameters::Nd();
1148  int NC2 = 2 * NC;
1149  int NCD2 = NC2 * ND;
1150 #endif
1151 
1152  const int Nvol = v1.nvol();
1153  const int Nvol_s = Nvol / CommonParameters::Nt();
1154  const int Nx = CommonParameters::Nx();
1155  const int Ny = CommonParameters::Ny();
1156  const int Nz = CommonParameters::Nz();
1157  const int Nt = CommonParameters::Nt();
1158  const int Lx = CommonParameters::Lx();
1159  const int Ly = CommonParameters::Ly();
1160  const int Lz = CommonParameters::Lz();
1161  const int Lt = CommonParameters::Lt();
1162 
1163  assert(Nvol == v2.nvol());
1164  assert(v1.nex() == 1);
1165  assert(v2.nex() == 1);
1166  assert(momentum_sink.size() == NDIM - 1);
1167  assert(source_position.size() == NDIM);
1168 
1169  const double *w1 = v1.ptr(0);
1170  const double *w2 = v2.ptr(0);
1171 
1172  int id1[ND];
1173  int id2[ND];
1174  for (int id = 0; id < ND; ++id) {
1175  id1[id] = id * NC2;
1176  id2[id] = gm_sink.index(id) * NC2;
1177  }
1178 
1179  double c_r[ND];
1180  double c_i[ND];
1181  for (int id = 0; id < ND; ++id) {
1182  c_r[id] = 0.0;
1183  c_i[id] = 0.0;
1184  }
1185 
1186  static const double PI = 4.0 * atan(1.0);
1187  std::vector<double> p_unit(ND - 1);
1188  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1189  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1190  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1191 
1192  std::vector<int> ipe(ND);
1193  ipe[0] = Communicator::ipe(0);
1194  ipe[1] = Communicator::ipe(1);
1195  ipe[2] = Communicator::ipe(2);
1196  ipe[3] = Communicator::ipe(3);
1197 
1198 
1199  for (int t = 0; t < Nt; ++t) {
1200  for (int z = 0; z < Nz; ++z) {
1201  for (int x = 0; x < Nx; ++x) {
1202  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1203 
1204  int x_global = x + ipe[0] * Nx;
1205  int z_global = z + ipe[2] * Nz;
1206  int t_global = t + ipe[3] * Nt;
1207 
1208  double p_x = p_unit[0] * (x_global - source_position[0]);
1209  double p_z = p_unit[1] * (z_global - source_position[2]);
1210  double p_t = p_unit[2] * (t_global - source_position[3]);
1211 
1212  double cos_p_xzt = cos(p_t + p_x + p_z);
1213  //double sin_p_xzt = sin(p_t + p_x + p_z);
1214  double sin_p_xzt = 0;
1215 
1216  for (int cc = 0; cc < NC; ++cc) {
1217  for (int id = 0; id < ND; ++id) {
1218  int ic1_r = 2 * cc + id1[id] + site;
1219  int ic2_r = 2 * cc + id2[id] + site;
1220 
1221  int ic1_i = 2 * cc + 1 + id1[id] + site;
1222  int ic2_i = 2 * cc + 1 + id2[id] + site;
1223 
1224  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1225  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1226  //double w1_w2_i = 0;
1227 
1228  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1229  c_r[id] += w1_w2_r * cos_p_xzt - w1_w2_i * sin_p_xzt;
1230  c_i[id] += w1_w2_r * sin_p_xzt + w1_w2_i * cos_p_xzt;
1231  }
1232  }
1233  }
1234  }
1235  }
1236 
1237  corr = cmplx(0.0, 0.0);
1238  for (int id = 0; id < ND; ++id) {
1239  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1240  }
1241 }
1242 
1243 
1244 //====================================================================
1245 void contract_at_y_cos(std::vector<dcomplex>& corr_global,
1246  const std::vector<int>& momentum_sink,
1247  const GammaMatrix& gm_sink,
1248  const std::vector<int>& source_position,
1249  const Field_F& v1, const Field_F& v2)
1250 {
1251  const int Ly = CommonParameters::Ly();
1252  const int Ny = CommonParameters::Ny();
1253 
1254  assert(corr_global.size() == Ly);
1255 
1256  std::vector<dcomplex> corr_local(Ny, 0.0);
1257  for (int y = 0; y < Ny; ++y) {
1258  dcomplex corr_y;
1259  contract_at_y_cos(corr_y, momentum_sink, gm_sink, source_position, v1, v2, y);
1260  corr_local[y] += corr_y;
1261  }
1262  global_corr_y(corr_global, corr_local);
1263 }
1264 
1265 
1266 //====================================================================
1267 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta
1268 void contract_at_z(dcomplex& corr,
1269  const GammaMatrix& gm_sink,
1270  const Field_F& v1, const Field_F& v2,
1271  const int z)
1272 {
1273 #if defined USE_GROUP_SU_N
1274  const int NC = CommonParameters::Nc();
1275  const int ND = CommonParameters::Nd();
1276  const int NC2 = 2 * NC;
1277  const int NCD2 = NC2 * ND;
1278 #endif
1279 
1280  int Nx = CommonParameters::Nx();
1281  int Ny = CommonParameters::Ny();
1282  int Nz = CommonParameters::Nz();
1283  int Nt = CommonParameters::Nt();
1284  const int Nvol = v1.nvol();
1285 
1286  assert(Nvol == v1.nvol());
1287  assert(Nvol == v2.nvol());
1288  assert(v1.nex() == 1);
1289  assert(v2.nex() == 1);
1290  assert(z < Nz);
1291 
1292  const double *w1 = v1.ptr(0);
1293  const double *w2 = v2.ptr(0);
1294 
1295  int id1[ND];
1296  int id2[ND];
1297  for (int id = 0; id < ND; ++id) {
1298  id1[id] = id * NC2;
1299  id2[id] = gm_sink.index(id) * NC2;
1300  }
1301 
1302  double c_r[ND];
1303  double c_i[ND];
1304  for (int id = 0; id < ND; ++id) {
1305  c_r[id] = 0.0;
1306  c_i[id] = 0.0;
1307  }
1308 
1309  for (int t = 0; t < Nt; ++t) {
1310  for (int y = 0; y < Ny; ++y) {
1311  for (int x = 0; x < Nx; ++x) {
1312  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1313  for (int cc = 0; cc < NC; ++cc) {
1314  for (int id = 0; id < ND; ++id) {
1315  int ic1_r = 2 * cc + id1[id] + site;
1316  int ic2_r = 2 * cc + id2[id] + site;
1317 
1318  int ic1_i = 2 * cc + 1 + id1[id] + site;
1319  int ic2_i = 2 * cc + 1 + id2[id] + site;
1320 
1321  c_r[id] += w1[ic2_r] * w2[ic1_r]
1322  + w1[ic2_i] * w2[ic1_i];
1323 
1324  c_i[id] += -w1[ic2_r] * w2[ic1_i]
1325  + w1[ic2_i] * w2[ic1_r];
1326  }
1327  }
1328  }
1329  }
1330  }
1331 
1332  corr = cmplx(0.0, 0.0);
1333  for (int id = 0; id < ND; ++id) {
1334  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1335  }
1336 }
1337 
1338 
1339 //====================================================================
1340 void contract_at_z(std::vector<dcomplex>& corr_global,
1341  const GammaMatrix& gm_sink,
1342  const Field_F& v1, const Field_F& v2)
1343 {
1344  const int Lz = CommonParameters::Lz();
1345  const int Nz = CommonParameters::Nz();
1346 
1347  assert(corr_global.size() == Lz);
1348 
1349  std::vector<dcomplex> corr_local(Nz, 0.0);
1350  for (int z = 0; z < Nz; ++z) {
1351  dcomplex corr_z;
1352  contract_at_z(corr_z, gm_sink, v1, v2, z);
1353  corr_local[z] += corr_z;
1354  }
1355  global_corr_z(corr_global, corr_local);
1356 }
1357 
1358 
1359 //====================================================================
1360 void contract_at_z(dcomplex& corr,
1361  const std::vector<int>& momentum_sink,
1362  const GammaMatrix& gm_sink,
1363  const std::vector<int>& source_position,
1364  const Field_F& v1, const Field_F& v2,
1365  const int z)
1366 {
1367 #if defined USE_GROUP_SU_N
1368  int NC = CommonParameters::Nc();
1369  int ND = CommonParameters::Nd();
1370  int NC2 = 2 * NC;
1371  int NCD2 = NC2 * ND;
1372 #endif
1373 
1374  const int Nvol = v1.nvol();
1375  const int Nvol_s = Nvol / CommonParameters::Nt();
1376  const int Nx = CommonParameters::Nx();
1377  const int Ny = CommonParameters::Ny();
1378  const int Nz = CommonParameters::Nz();
1379  const int Nt = CommonParameters::Nt();
1380  const int Lx = CommonParameters::Lx();
1381  const int Ly = CommonParameters::Ly();
1382  const int Lz = CommonParameters::Lz();
1383  const int Lt = CommonParameters::Lt();
1384 
1385  assert(Nvol == v2.nvol());
1386  assert(v1.nex() == 1);
1387  assert(v2.nex() == 1);
1388  assert(momentum_sink.size() == NDIM - 1);
1389  assert(source_position.size() == NDIM);
1390 
1391  const double *w1 = v1.ptr(0);
1392  const double *w2 = v2.ptr(0);
1393 
1394  int id1[ND];
1395  int id2[ND];
1396  for (int id = 0; id < ND; ++id) {
1397  id1[id] = id * NC2;
1398  id2[id] = gm_sink.index(id) * NC2;
1399  }
1400 
1401  double c_r[ND];
1402  double c_i[ND];
1403  for (int id = 0; id < ND; ++id) {
1404  c_r[id] = 0.0;
1405  c_i[id] = 0.0;
1406  }
1407 
1408  static const double PI = 4.0 * atan(1.0);
1409  std::vector<double> p_unit(ND - 1);
1410  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1411  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1412  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1413 
1414  std::vector<int> ipe(ND);
1415  ipe[0] = Communicator::ipe(0);
1416  ipe[1] = Communicator::ipe(1);
1417  ipe[2] = Communicator::ipe(2);
1418  ipe[3] = Communicator::ipe(3);
1419 
1420  for (int t = 0; t < Nt; ++t) {
1421  for (int y = 0; y < Ny; ++y) {
1422  for (int x = 0; x < Nx; ++x) {
1423  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1424 
1425  int x_global = x + ipe[0] * Nx;
1426  int y_global = y + ipe[1] * Ny;
1427  int t_global = t + ipe[3] * Nt;
1428 
1429  double p_x = p_unit[0] * (x_global - source_position[0]);
1430  double p_y = p_unit[1] * (y_global - source_position[1]);
1431  double p_t = p_unit[2] * (t_global - source_position[3]);
1432 
1433  double cos_p_xyt = cos(p_t + p_x + p_y);
1434  double sin_p_xyt = sin(p_t + p_x + p_y);
1435 
1436  for (int cc = 0; cc < NC; ++cc) {
1437  for (int id = 0; id < ND; ++id) {
1438  int ic1_r = 2 * cc + id1[id] + site;
1439  int ic2_r = 2 * cc + id2[id] + site;
1440 
1441  int ic1_i = 2 * cc + 1 + id1[id] + site;
1442  int ic2_i = 2 * cc + 1 + id2[id] + site;
1443 
1444  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1445  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1446 
1447  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1448  c_r[id] += w1_w2_r * cos_p_xyt - w1_w2_i * sin_p_xyt;
1449  c_i[id] += w1_w2_r * sin_p_xyt + w1_w2_i * cos_p_xyt;
1450  }
1451  }
1452  }
1453  }
1454  }
1455 
1456  corr = cmplx(0.0, 0.0);
1457  for (int id = 0; id < ND; ++id) {
1458  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1459  }
1460 }
1461 
1462 
1463 //====================================================================
1464 void contract_at_z(std::vector<dcomplex>& corr_global,
1465  const std::vector<int>& momentum_sink,
1466  const GammaMatrix& gm_sink,
1467  const std::vector<int>& source_position,
1468  const Field_F& v1, const Field_F& v2)
1469 {
1470  const int Lz = CommonParameters::Lz();
1471  const int Nz = CommonParameters::Nz();
1472 
1473  assert(corr_global.size() == Lz);
1474 
1475  std::vector<dcomplex> corr_local(Nz, 0.0);
1476  for (int z = 0; z < Nz; ++z) {
1477  dcomplex corr_z;
1478  contract_at_z(corr_z, momentum_sink, gm_sink, source_position, v1, v2, z);
1479  corr_local[z] += corr_z;
1480  }
1481  global_corr_z(corr_global, corr_local);
1482 }
1483 
1484 
1485 //====================================================================
1486 void contract_at_z_cos(dcomplex& corr,
1487  const std::vector<int>& momentum_sink,
1488  const GammaMatrix& gm_sink,
1489  const std::vector<int>& source_position,
1490  const Field_F& v1, const Field_F& v2,
1491  const int z)
1492 {
1493 #if defined USE_GROUP_SU_N
1494  int NC = CommonParameters::Nc();
1495  int ND = CommonParameters::Nd();
1496  int NC2 = 2 * NC;
1497  int NCD2 = NC2 * ND;
1498 #endif
1499 
1500  const int Nvol = v1.nvol();
1501  const int Nvol_s = Nvol / CommonParameters::Nt();
1502  const int Nx = CommonParameters::Nx();
1503  const int Ny = CommonParameters::Ny();
1504  const int Nz = CommonParameters::Nz();
1505  const int Nt = CommonParameters::Nt();
1506  const int Lx = CommonParameters::Lx();
1507  const int Ly = CommonParameters::Ly();
1508  const int Lz = CommonParameters::Lz();
1509  const int Lt = CommonParameters::Lt();
1510 
1511  assert(Nvol == v2.nvol());
1512  assert(v1.nex() == 1);
1513  assert(v2.nex() == 1);
1514  assert(momentum_sink.size() == NDIM - 1);
1515  assert(source_position.size() == NDIM);
1516 
1517  const double *w1 = v1.ptr(0);
1518  const double *w2 = v2.ptr(0);
1519 
1520  int id1[ND];
1521  int id2[ND];
1522  for (int id = 0; id < ND; ++id) {
1523  id1[id] = id * NC2;
1524  id2[id] = gm_sink.index(id) * NC2;
1525  }
1526 
1527  double c_r[ND];
1528  double c_i[ND];
1529  for (int id = 0; id < ND; ++id) {
1530  c_r[id] = 0.0;
1531  c_i[id] = 0.0;
1532  }
1533 
1534  static const double PI = 4.0 * atan(1.0);
1535  std::vector<double> p_unit(ND - 1);
1536  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1537  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1538  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1539 
1540  std::vector<int> ipe(ND);
1541  ipe[0] = Communicator::ipe(0);
1542  ipe[1] = Communicator::ipe(1);
1543  ipe[2] = Communicator::ipe(2);
1544  ipe[3] = Communicator::ipe(3);
1545 
1546 
1547  for (int t = 0; t < Nt; ++t) {
1548  for (int y = 0; y < Ny; ++y) {
1549  for (int x = 0; x < Nx; ++x) {
1550  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1551 
1552  int x_global = x + ipe[0] * Nx;
1553  int y_global = y + ipe[1] * Ny;
1554  int t_global = t + ipe[3] * Nt;
1555 
1556  double p_x = p_unit[0] * (x_global - source_position[0]);
1557  double p_y = p_unit[1] * (y_global - source_position[1]);
1558  double p_t = p_unit[2] * (t_global - source_position[2]);
1559 
1560  double cos_p_xyt = cos(p_t + p_x + p_y);
1561  //double sin_p_xyt = sin(p_t + p_x + p_y);
1562  double sin_p_xyt = 0;
1563 
1564  for (int cc = 0; cc < NC; ++cc) {
1565  for (int id = 0; id < ND; ++id) {
1566  int ic1_r = 2 * cc + id1[id] + site;
1567  int ic2_r = 2 * cc + id2[id] + site;
1568 
1569  int ic1_i = 2 * cc + 1 + id1[id] + site;
1570  int ic2_i = 2 * cc + 1 + id2[id] + site;
1571 
1572  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1573  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1574  //double w1_w2_i = 0;
1575 
1576  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1577  c_r[id] += w1_w2_r * cos_p_xyt - w1_w2_i * sin_p_xyt;
1578  c_i[id] += w1_w2_r * sin_p_xyt + w1_w2_i * cos_p_xyt;
1579  }
1580  }
1581  }
1582  }
1583  }
1584 
1585  corr = cmplx(0.0, 0.0);
1586  for (int id = 0; id < ND; ++id) {
1587  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1588  }
1589 }
1590 
1591 
1592 //====================================================================
1593 void contract_at_z_cos(std::vector<dcomplex>& corr_global,
1594  const std::vector<int>& momentum_sink,
1595  const GammaMatrix& gm_sink,
1596  const std::vector<int>& source_position,
1597  const Field_F& v1, const Field_F& v2)
1598 {
1599  const int Lz = CommonParameters::Lz();
1600  const int Nz = CommonParameters::Nz();
1601 
1602  assert(corr_global.size() == Lz);
1603 
1604  std::vector<dcomplex> corr_local(Nz, 0.0);
1605  for (int z = 0; z < Nz; ++z) {
1606  dcomplex corr_z;
1607  contract_at_z_cos(corr_z, momentum_sink, gm_sink, source_position, v1, v2, z);
1608  corr_local[z] += corr_z;
1609  }
1610  global_corr_z(corr_global, corr_local);
1611 }
1612 
1613 
1614 //====================================================================
1615 void global_corr_x(std::vector<dcomplex>& corr_global,
1616  std::vector<dcomplex>& corr_local)
1617 {
1618  int Lx = CommonParameters::Lx();
1619  int Nx = CommonParameters::Nx();
1620 
1621  assert(corr_global.size() == Lx);
1622  assert(corr_local.size() == Nx);
1623 
1624  std::vector<dcomplex> corr_tmp(Lx, 0);
1625 
1626  int ipex = Communicator::ipe(0);
1627 
1628  for (int x = 0; x < Nx; ++x) {
1629  int x2 = x + ipex * Nx;
1630  corr_tmp[x2] = corr_local[x];
1631  }
1632 
1633  for (int x = 0; x < Lx; ++x) {
1634  double crr = Communicator::reduce_sum(real(corr_tmp[x]));
1635  double cri = Communicator::reduce_sum(imag(corr_tmp[x]));
1636  corr_global[x] = cmplx(crr, cri);
1637  }
1638 }
1639 
1640 
1641 //====================================================================
1642 void global_corr_y(std::vector<dcomplex>& corr_global,
1643  std::vector<dcomplex>& corr_local)
1644 {
1645  int Ly = CommonParameters::Ly();
1646  int Ny = CommonParameters::Ny();
1647 
1648  assert(corr_global.size() == Ly);
1649  assert(corr_local.size() == Ny);
1650 
1651  std::vector<dcomplex> corr_tmp(Ly, 0);
1652 
1653  int ipey = Communicator::ipe(1);
1654 
1655  for (int y = 0; y < Ny; ++y) {
1656  int y2 = y + ipey * Ny;
1657  corr_tmp[y2] = corr_local[y];
1658  }
1659 
1660  for (int y = 0; y < Ly; ++y) {
1661  double crr = Communicator::reduce_sum(real(corr_tmp[y]));
1662  double cri = Communicator::reduce_sum(imag(corr_tmp[y]));
1663  corr_global[y] = cmplx(crr, cri);
1664  }
1665 }
1666 
1667 
1668 //====================================================================
1669 void global_corr_z(std::vector<dcomplex>& corr_global,
1670  std::vector<dcomplex>& corr_local)
1671 {
1672  int Lz = CommonParameters::Lz();
1673  int Nz = CommonParameters::Nz();
1674 
1675  assert(corr_global.size() == Lz);
1676  assert(corr_local.size() == Nz);
1677 
1678  std::vector<dcomplex> corr_tmp(Lz);
1679 
1680  int ipez = Communicator::ipe(2);
1681 
1682  for (int z = 0; z < Nz; ++z) {
1683  int z2 = z + ipez * Nz;
1684  corr_tmp[z2] = corr_local[z];
1685  }
1686 
1687  for (int z = 0; z < Lz; ++z) {
1688  double crr = Communicator::reduce_sum(real(corr_tmp[z]));
1689  double cri = Communicator::reduce_sum(imag(corr_tmp[z]));
1690  corr_global[z] = cmplx(crr, cri);
1691  }
1692 }
1693 
1694 
1695 //====================================================================
1696 void global_corr_t(std::vector<dcomplex>& corr_global,
1697  std::vector<dcomplex>& corr_local)
1698 {
1699  int Lt = CommonParameters::Lt();
1700  int Nt = CommonParameters::Nt();
1701 
1702  assert(corr_global.size() == Lt);
1703  assert(corr_local.size() == Nt);
1704 
1705  std::vector<dcomplex> corr_tmp(Lt, 0);
1706 
1707  int ipe_t = Communicator::ipe(3);
1708 
1709  for (int t = 0; t < Nt; ++t) {
1710  int t_global = t + ipe_t * Nt;
1711  corr_tmp[t_global] = corr_local[t];
1712  }
1713 
1714  for (int t_global = 0; t_global < Lt; ++t_global) {
1715  double cr_r = Communicator::reduce_sum(real(corr_tmp[t_global]));
1716  double cr_i = Communicator::reduce_sum(imag(corr_tmp[t_global]));
1717  corr_global[t_global] = cmplx(cr_r, cr_i);
1718  }
1719 }
1720 
1721 
1722 //====================================================================
1723 //============================================================END=====
NC2
#define NC2
Definition: field_F_imp_SU2-inc.h:3
Field_F::cmp_i
double cmp_i(const int cc, const int s, const int site, const int e=0) const
Definition: field_F.h:100
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:154
global_corr_z
void global_corr_z(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in z to global.
Definition: contract_4spinor.cpp:1686
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:53
global_corr_x
void global_corr_x(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in x to global.
Definition: contract_4spinor.cpp:1632
contract_at_t
void contract_at_t(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, const int time)
Contraction of hadron for 4-spinor fermion.
Definition: contract_4spinor.cpp:37
ND
#define ND
Definition: field_F_imp_SU2-inc.h:5
Field::nvol
int nvol() const
Definition: field.h:128
contract_at_x_cos
void contract_at_x_cos(dcomplex &corr, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const std::vector< int > &source_position, const Field_F &v1, const Field_F &v2, const int x)
contraction for meson at a given x with Fourier transformation, where (p_y,p_z,p_t) is given by momen...
Definition: contract_4spinor.cpp:694
NDIM
#define NDIM
Definition: contract_4spinor.cpp:19
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Communicator::reduce_sum
static int reduce_sum(int count, dcomplex *recv_buf, dcomplex *send_buf, int pattern=0)
make a global sum of an array of dcomplex over the communicator. pattern specifies the dimensions to ...
Definition: communicator.cpp:263
contract_at_y_cos
void contract_at_y_cos(dcomplex &corr, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const std::vector< int > &source_position, const Field_F &v1, const Field_F &v2, const int y)
contraction for meson at a given y with Fourier transformation, where (p_x,p_z,p_t) is given by momen...
Definition: contract_4spinor.cpp:1155
contract_at_z_cos
void contract_at_z_cos(dcomplex &corr, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const std::vector< int > &source_position, const Field_F &v1, const Field_F &v2, const int z)
contraction for meson at a given z with Fourier transformation, where (p_x,p_y,p_t) is given by momen...
Definition: contract_4spinor.cpp:1503
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
contract_at_z
void contract_at_z(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, const int z)
contraction for meson at a given z.
Definition: contract_4spinor.cpp:1285
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
NC
#define NC
Definition: field_F_imp_SU2-inc.h:2
index_lex.h
contract_at_x
void contract_at_x(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, int x)
Definition: contract_4spinor.cpp:477
global_corr_t
void global_corr_t(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in t to global.
Definition: contract_4spinor.cpp:1713
Field::nex
int nex() const
Definition: field.h:129
Index_lex
Lexical site index.
Definition: index_lex.h:34
global_corr_y
void global_corr_y(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in y to global.
Definition: contract_4spinor.cpp:1659
NCD2
#define NCD2
Definition: field_F_imp_SU2-inc.h:7
contract_4spinor.h
contract_at_t_cos
void contract_at_t_cos(dcomplex &corr, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const std::vector< int > &source_position, const Field_F &v1, const Field_F &v2, const int time)
Definition: contract_4spinor.cpp:245
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
contract_at_y
void contract_at_y(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, const int y)
contraction for meson at a given y.
Definition: contract_4spinor.cpp:937
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83
Field_F::cmp_r
double cmp_r(const int cc, const int s, const int site, const int e=0) const
Definition: field_F.h:94