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