Bridge++  Ver. 2.0.2
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 Lx = CommonParameters::Lx();
127  const int Ly = CommonParameters::Ly();
128  const int Lz = CommonParameters::Lz();
129  const int Ndim = CommonParameters::Ndim();
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  const int Ndim = CommonParameters::Ndim();
248 
249  assert(Nvol == v2.nvol());
250  assert(v1.nex() == 1);
251  assert(v2.nex() == 1);
252  assert(momentum_sink.size() == Ndim - 1);
253 
254  const double *w1 = v1.ptr(0);
255  const double *w2 = v2.ptr(0);
256 
257  int id1[ND];
258  int id2[ND];
259  for (int id = 0; id < ND; ++id) {
260  id1[id] = id * NC2;
261  id2[id] = gm_sink.index(id) * NC2;
262  }
263 
264  double c_r[ND];
265  double c_i[ND];
266  for (int id = 0; id < ND; ++id) {
267  c_r[id] = 0.0;
268  c_i[id] = 0.0;
269  }
270 
271  static const double PI = 4.0 * atan(1.0);
272  std::vector<double> p_unit(ND - 1);
273  p_unit[0] = (2.0 * PI / Lx) * momentum_sink[0];
274  p_unit[1] = (2.0 * PI / Ly) * momentum_sink[1];
275  p_unit[2] = (2.0 * PI / Lz) * momentum_sink[2];
276 
277  std::vector<int> ipe(ND - 1);
278  ipe[0] = Communicator::ipe(0);
279  ipe[1] = Communicator::ipe(1);
280  ipe[2] = Communicator::ipe(2);
281 
282 
283  for (int ss = 0; ss < Nvol_s; ++ss) {
284  int site = NCD2 * (ss + time * Nvol_s);
285 
286  int x = ss % Nx;
287  int y = ss % (Nx * Ny) / Nx;
288  int z = ss % (Nx * Ny * Nz) / (Nx * Ny);
289 
290  int x_global = x + ipe[0] * Nx;
291  int y_global = y + ipe[1] * Ny;
292  int z_global = z + ipe[2] * Nz;
293 
294  double p_x = p_unit[0] * (x_global - source_position[0]);
295  double p_y = p_unit[1] * (y_global - source_position[1]);
296  double p_z = p_unit[2] * (z_global - source_position[2]);
297 
298  double cos_p_xyz = cos(p_x + p_y + p_z);
299  //double sin_p_xyz = sin(p_x + p_y + p_z);
300  double sin_p_xyz = 0;
301 
302  for (int cc = 0; cc < NC; ++cc) {
303  for (int id = 0; id < ND; ++id) {
304  int ic1_r = 2 * cc + id1[id] + site;
305  int ic2_r = 2 * cc + id2[id] + site;
306 
307  int ic1_i = 2 * cc + 1 + id1[id] + site;
308  int ic2_i = 2 * cc + 1 + id2[id] + site;
309 
310  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
311  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
312  //double w1_w2_i = 0;
313 
314  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
315  c_r[id] += w1_w2_r * cos_p_xyz - w1_w2_i * sin_p_xyz;
316  c_i[id] += w1_w2_r * sin_p_xyz + w1_w2_i * cos_p_xyz;
317  }
318  }
319  }
320 
321  corr = cmplx(0.0, 0.0);
322  for (int id = 0; id < ND; ++id) {
323  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
324  }
325 }
326 
327 
328 //====================================================================
329 void contract_at_t_cos(std::vector<dcomplex>& corr_global,
330  const std::vector<int>& momentum_sink,
331  const GammaMatrix& gm_sink,
332  const std::vector<int>& source_position,
333  const Field_F& v1, const Field_F& v2)
334 {
335  const int Lt = CommonParameters::Lt();
336  const int Nt = CommonParameters::Nt();
337 
338  assert(corr_global.size() == Lt);
339 
340  std::vector<dcomplex> corr_local(Nt, 0.0);
341  for (int t = 0; t < Nt; ++t) {
342  dcomplex corr_t;
343  contract_at_t_cos(corr_t, momentum_sink, gm_sink, source_position, v1, v2, t);
344  corr_local[t] += corr_t;
345  }
346  global_corr_t(corr_global, corr_local);
347 }
348 
349 
350 //====================================================================
351 void contract_at_t(dcomplex& corr,
352  const GammaMatrix& gm_sink,
353  const int i_alpha,
354  const Field_F& v1, const Field_F& v2, const Field_F& v3,
355  const int time)
356 {
357 #if defined USE_GROUP_SU3
358  const int Nc = CommonParameters::Nc();
359  const int Nd = CommonParameters::Nd();
360  const int Nx = CommonParameters::Nx();
361  const int Ny = CommonParameters::Ny();
362  const int Nz = CommonParameters::Nz();
363  const int Nt = CommonParameters::Nt();
364  const int Nvol = CommonParameters::Nvol();
365 
366  assert(Nc == 3);
367  assert(Nvol == v1.nvol());
368  assert(Nvol == v2.nvol());
369  assert(Nvol == v3.nvol());
370  assert(time < Nt);
371 
372  Index_lex index;
373  std::vector<int> gm_index(Nd);
374  std::vector<double> c_r(Nd), c_i(Nd);
375 
376  for (int i = 0; i < Nd; ++i) {
377  gm_index[i] = gm_sink.index(i);
378  c_r[i] = 0.0;
379  c_i[i] = 0.0;
380  }
381 
382  for (int z = 0; z < Nz; ++z) {
383  for (int y = 0; y < Ny; ++y) {
384  for (int x = 0; x < Nx; ++x) {
385  int site = index.site(x, y, z, time);
386 
387  for (int d1 = 0; d1 < Nd; ++d1) {
388  int d2 = gm_index[d1];
389  int d3 = i_alpha;
390 
391  c_r[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
392  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C3, d3, site)
393  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
394  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C3, d3, site);
395  c_i[d1] += (v1.cmp_r(C1, d1, site) * v2.cmp_r(C2, d2, site)
396  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C3, d3, site)
397  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C2, d2, site)
398  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C3, d3, site);
399 
400  c_r[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
401  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C1, d3, site)
402  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
403  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C1, d3, site);
404  c_i[d1] += (v1.cmp_r(C2, d1, site) * v2.cmp_r(C3, d2, site)
405  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C1, d3, site)
406  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C3, d2, site)
407  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C1, d3, site);
408 
409  c_r[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
410  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C2, d3, site)
411  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
412  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C2, d3, site);
413  c_i[d1] += (v1.cmp_r(C3, d1, site) * v2.cmp_r(C1, d2, site)
414  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C2, d3, site)
415  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C1, d2, site)
416  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C2, d3, site);
417 
418  c_r[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
419  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_r(C1, d3, site)
420  - (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
421  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_i(C1, d3, site);
422  c_i[d1] -= (v1.cmp_r(C3, d1, site) * v2.cmp_r(C2, d2, site)
423  - v1.cmp_i(C3, d1, site) * v2.cmp_i(C2, d2, site)) * v3.cmp_i(C1, d3, site)
424  + (v1.cmp_r(C3, d1, site) * v2.cmp_i(C2, d2, site)
425  + v1.cmp_i(C3, d1, site) * v2.cmp_r(C2, d2, site)) * v3.cmp_r(C1, d3, site);
426 
427  c_r[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
428  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_r(C3, d3, site)
429  - (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
430  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_i(C3, d3, site);
431  c_i[d1] -= (v1.cmp_r(C2, d1, site) * v2.cmp_r(C1, d2, site)
432  - v1.cmp_i(C2, d1, site) * v2.cmp_i(C1, d2, site)) * v3.cmp_i(C3, d3, site)
433  + (v1.cmp_r(C2, d1, site) * v2.cmp_i(C1, d2, site)
434  + v1.cmp_i(C2, d1, site) * v2.cmp_r(C1, d2, site)) * v3.cmp_r(C3, d3, site);
435 
436  c_r[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
437  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_r(C2, d3, site)
438  - (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
439  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_i(C2, d3, site);
440  c_i[d1] -= (v1.cmp_r(C1, d1, site) * v2.cmp_r(C3, d2, site)
441  - v1.cmp_i(C1, d1, site) * v2.cmp_i(C3, d2, site)) * v3.cmp_i(C2, d3, site)
442  + (v1.cmp_r(C1, d1, site) * v2.cmp_i(C3, d2, site)
443  + v1.cmp_i(C1, d1, site) * v2.cmp_r(C3, d2, site)) * v3.cmp_r(C2, d3, site);
444  }
445  }
446  }
447  }
448 
449  corr = cmplx(0.0, 0.0);
450  for (int s0 = 0; s0 < Nd; ++s0) {
451  corr += gm_sink.value(s0) * cmplx(c_r[s0], c_i[s0]);
452  }
453 #endif // (USE_GROUP_SU3)
454 }
455 
456 
457 //====================================================================
458 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta at x
459 void contract_at_x(dcomplex& corr, const GammaMatrix& gm_sink,
460  const Field_F& v1, const Field_F& v2,
461  int x)
462 {
463 #if defined USE_GROUP_SU_N
464  const int NC = CommonParameters::Nc();
465  const int ND = CommonParameters::Nd();
466  const int NC2 = 2 * NC;
467  const int NCD2 = NC2 * ND;
468 #endif
469 
470  const int Nvol = v1.nvol();
471  const int Nvol_s = Nvol / CommonParameters::Nt();
472  const int Nx = CommonParameters::Nx();
473  const int Ny = CommonParameters::Ny();
474  const int Nz = CommonParameters::Nz();
475  const int Nt = CommonParameters::Nt();
476 
477  assert(Nvol == v2.nvol());
478  assert(v1.nex() == 1);
479  assert(v2.nex() == 1);
480 
481  const double *w1 = v1.ptr(0);
482  const double *w2 = v2.ptr(0);
483 
484  int id1[ND];
485  int id2[ND];
486  for (int id = 0; id < ND; ++id) {
487  id1[id] = id * NC2;
488  id2[id] = gm_sink.index(id) * NC2;
489  }
490 
491  double c_r[ND];
492  double c_i[ND];
493  for (int id = 0; id < ND; ++id) {
494  c_r[id] = 0.0;
495  c_i[id] = 0.0;
496  }
497 
498  for (int t = 0; t < Nt; ++t) {
499  for (int z = 0; z < Nz; ++z) {
500  for (int y = 0; y < Ny; ++y) {
501  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
502 
503  for (int cc = 0; cc < NC; ++cc) {
504  for (int id = 0; id < ND; ++id) {
505  int ic1_r = 2 * cc + id1[id] + site;
506  int ic2_r = 2 * cc + id2[id] + site;
507 
508  int ic1_i = 2 * cc + 1 + id1[id] + site;
509  int ic2_i = 2 * cc + 1 + id2[id] + site;
510 
511  c_r[id] += w1[ic2_r] * w2[ic1_r]
512  + w1[ic2_i] * w2[ic1_i];
513 
514  c_i[id] += -w1[ic2_r] * w2[ic1_i]
515  + w1[ic2_i] * w2[ic1_r];
516  }
517  }
518  }
519  }
520  }
521 
522  corr = cmplx(0.0, 0.0);
523  for (int id = 0; id < ND; ++id) {
524  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
525  }
526 }
527 
528 
529 //====================================================================
530 void contract_at_x(std::vector<dcomplex>& corr_global,
531  const GammaMatrix& gm_sink,
532  const Field_F& v1, const Field_F& v2)
533 {
534  const int Lx = CommonParameters::Lx();
535  const int Nx = CommonParameters::Nx();
536 
537  assert(corr_global.size() == Lx);
538 
539  std::vector<dcomplex> corr_local(Nx, 0.0);
540  for (int x = 0; x < Nx; ++x) {
541  dcomplex corr_x;
542  contract_at_x(corr_x, gm_sink, v1, v2, x);
543  corr_local[x] += corr_x;
544  }
545  global_corr_x(corr_global, corr_local);
546 }
547 
548 
549 //====================================================================
550 void contract_at_x(dcomplex& corr,
551  const std::vector<int>& momentum_sink,
552  const GammaMatrix& gm_sink,
553  const std::vector<int>& source_position,
554  const Field_F& v1, const Field_F& v2,
555  const int x)
556 {
557 #if defined USE_GROUP_SU_N
558  int NC = CommonParameters::Nc();
559  int ND = CommonParameters::Nd();
560  int NC2 = 2 * NC;
561  int NCD2 = NC2 * ND;
562 #endif
563 
564  const int Nvol = v1.nvol();
565  const int Nvol_s = Nvol / CommonParameters::Nt();
566  const int Nx = CommonParameters::Nx();
567  const int Ny = CommonParameters::Ny();
568  const int Nz = CommonParameters::Nz();
569  const int Nt = CommonParameters::Nt();
570  const int Lx = CommonParameters::Lx();
571  const int Ly = CommonParameters::Ly();
572  const int Lz = CommonParameters::Lz();
573  const int Lt = CommonParameters::Lt();
574  const int Ndim = CommonParameters::Ndim();
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  const int Ndim = CommonParameters::Ndim();
702 
703  assert(Nvol == v2.nvol());
704  assert(v1.nex() == 1);
705  assert(v2.nex() == 1);
706  assert(momentum_sink.size() == Ndim - 1);
707  assert(source_position.size() == Ndim);
708 
709  const double *w1 = v1.ptr(0);
710  const double *w2 = v2.ptr(0);
711 
712  int id1[ND];
713  int id2[ND];
714  for (int id = 0; id < ND; ++id) {
715  id1[id] = id * NC2;
716  id2[id] = gm_sink.index(id) * NC2;
717  }
718 
719  double c_r[ND];
720  double c_i[ND];
721  for (int id = 0; id < ND; ++id) {
722  c_r[id] = 0.0;
723  c_i[id] = 0.0;
724  }
725 
726  static const double PI = 4.0 * atan(1.0);
727  std::vector<double> p_unit(ND - 1);
728  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
729  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
730  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
731 
732  std::vector<int> ipe(ND);
733  ipe[0] = Communicator::ipe(0);
734  ipe[1] = Communicator::ipe(1);
735  ipe[2] = Communicator::ipe(2);
736  ipe[3] = Communicator::ipe(3);
737 
738 
739  for (int t = 0; t < Nt; ++t) {
740  for (int z = 0; z < Nz; ++z) {
741  for (int y = 0; y < Ny; ++y) {
742  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
743 
744  int y_global = y + ipe[1] * Ny;
745  int z_global = z + ipe[2] * Nz;
746  int t_global = t + ipe[3] * Nt;
747 
748  double p_y = p_unit[0] * (y_global - source_position[1]);
749  double p_z = p_unit[1] * (z_global - source_position[2]);
750  double p_t = p_unit[2] * (t_global - source_position[3]);
751 
752  double cos_p_yzt = cos(p_t + p_y + p_z);
753  //double sin_p_yzt = sin(p_t + p_y + p_z);
754  double sin_p_yzt = 0;
755 
756  for (int cc = 0; cc < NC; ++cc) {
757  for (int id = 0; id < ND; ++id) {
758  int ic1_r = 2 * cc + id1[id] + site;
759  int ic2_r = 2 * cc + id2[id] + site;
760 
761  int ic1_i = 2 * cc + 1 + id1[id] + site;
762  int ic2_i = 2 * cc + 1 + id2[id] + site;
763 
764  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
765  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
766  //double w1_w2_i = 0;
767 
768  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
769  c_r[id] += w1_w2_r * cos_p_yzt - w1_w2_i * sin_p_yzt;
770  c_i[id] += w1_w2_r * sin_p_yzt + w1_w2_i * cos_p_yzt;
771  }
772  }
773  }
774  }
775  }
776 
777  corr = cmplx(0.0, 0.0);
778  for (int id = 0; id < ND; ++id) {
779  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
780  }
781 }
782 
783 
784 //====================================================================
785 void contract_at_x_cos(std::vector<dcomplex>& corr_global,
786  const std::vector<int>& momentum_sink,
787  const GammaMatrix& gm_sink,
788  const std::vector<int>& source_position,
789  const Field_F& v1, const Field_F& v2)
790 {
791  const int Lx = CommonParameters::Lx();
792  const int Nx = CommonParameters::Nx();
793 
794  assert(corr_global.size() == Lx);
795 
796  std::vector<dcomplex> corr_local(Nx, 0.0);
797  for (int x = 0; x < Nx; ++x) {
798  dcomplex corr_x;
799  contract_at_x_cos(corr_x, momentum_sink, gm_sink, source_position, v1, v2, x);
800  corr_local[x] += corr_x;
801  }
802  global_corr_x(corr_global, corr_local);
803 }
804 
805 
806 //====================================================================
807 void contract_at_x(dcomplex& corr,
808  const GammaMatrix& gm_sink,
809  const int i_alpha,
810  const Field_F& v1, const Field_F& v2, const Field_F& v3,
811  const int x)
812 {
813 #if defined USE_GROUP_SU3
814  const int Nvol = v1.nvol();
815  const int Nvol_s = Nvol / CommonParameters::Nt();
816  const int Nx = CommonParameters::Nx();
817  const int Ny = CommonParameters::Ny();
818  const int Nz = CommonParameters::Nz();
819  const int Nt = CommonParameters::Nt();
820 
821  assert(Nvol == v2.nvol());
822  assert(Nvol == v3.nvol());
823  assert(v1.nex() == 1);
824  assert(v2.nex() == 1);
825  assert(v3.nex() == 1);
826 
827  const double *w1 = v1.ptr(0);
828  const double *w2 = v2.ptr(0);
829  const double *w3 = v3.ptr(0);
830 
831  int id1[ND];
832  int id2[ND];
833  for (int id = 0; id < ND; ++id) {
834  id1[id] = id * NC2;
835  id2[id] = gm_sink.index(id) * NC2;
836  }
837  int id3 = i_alpha * NC2;
838 
839  double c_r[ND];
840  double c_i[ND];
841  for (int id = 0; id < ND; ++id) {
842  c_r[id] = 0.0;
843  c_i[id] = 0.0;
844  }
845 
846 
847  for (int t = 0; t < Nt; ++t) {
848  for (int z = 0; z < Nz; ++z) {
849  for (int y = 0; y < Ny; ++y) {
850  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
851 
852  for (int id = 0; id < ND; ++id) {
853  int ic11_r = C1 + id1[id] + site;
854  int ic22_r = C2 + id2[id] + site;
855  int ic33_r = C3 + id3 + site;
856 
857  int ic11_i = C1 + 1 + id1[id] + site;
858  int ic22_i = C2 + 1 + id2[id] + site;
859  int ic33_i = C3 + 1 + id3 + site;
860 
861  int ic21_r = C2 + id1[id] + site;
862  int ic32_r = C3 + id2[id] + site;
863  int ic13_r = C1 + id3 + site;
864 
865  int ic21_i = C2 + 1 + id1[id] + site;
866  int ic32_i = C3 + 1 + id2[id] + site;
867  int ic13_i = C1 + 1 + id3 + site;
868 
869  int ic31_r = C3 + id1[id] + site;
870  int ic12_r = C1 + id2[id] + site;
871  int ic23_r = C2 + id3 + site;
872 
873  int ic31_i = C3 + 1 + id1[id] + site;
874  int ic12_i = C1 + 1 + id2[id] + site;
875  int ic23_i = C2 + 1 + id3 + site;
876 
877  c_r[id] += (w1[ic11_r] * w2[ic22_r] - w1[ic11_i] * w2[ic22_i]) * w3[ic33_r]
878  - (w1[ic11_r] * w2[ic22_i] + w1[ic11_i] * w2[ic22_r]) * w3[ic33_i];
879  c_i[id] += (w1[ic11_r] * w2[ic22_r] - w1[ic11_i] * w2[ic22_i]) * w3[ic33_i]
880  + (w1[ic11_r] * w2[ic22_i] + w1[ic11_i] * w2[ic22_r]) * w3[ic33_r];
881 
882  c_r[id] += (w1[ic21_r] * w2[ic32_r] - w1[ic21_i] * w2[ic32_i]) * w3[ic13_r]
883  - (w1[ic21_r] * w2[ic32_i] + w1[ic21_i] * w2[ic32_r]) * w3[ic13_i];
884  c_i[id] += (w1[ic21_r] * w2[ic32_r] - w1[ic21_i] * w2[ic32_i]) * w3[ic13_i]
885  + (w1[ic21_r] * w2[ic32_i] + w1[ic21_i] * w2[ic32_r]) * w3[ic13_r];
886 
887  c_r[id] += (w1[ic31_r] * w2[ic12_r] - w1[ic31_i] * w2[ic12_i]) * w3[ic23_r]
888  - (w1[ic31_r] * w2[ic12_i] + w1[ic31_i] * w2[ic12_r]) * w3[ic23_i];
889  c_i[id] += (w1[ic31_r] * w2[ic12_r] - w1[ic31_i] * w2[ic12_i]) * w3[ic23_i]
890  + (w1[ic31_r] * w2[ic12_i] + w1[ic31_i] * w2[ic12_r]) * w3[ic23_r];
891 
892  c_r[id] -= (w1[ic31_r] * w2[ic22_r] - w1[ic31_i] * w2[ic22_i]) * w3[ic13_r]
893  - (w1[ic31_r] * w2[ic22_i] + w1[ic31_i] * w2[ic22_r]) * w3[ic13_i];
894  c_i[id] -= (w1[ic31_r] * w2[ic22_r] - w1[ic31_i] * w2[ic22_i]) * w3[ic13_i]
895  + (w1[ic31_r] * w2[ic22_i] + w1[ic31_i] * w2[ic22_r]) * w3[ic13_r];
896 
897  c_r[id] -= (w1[ic21_r] * w2[ic12_r] - w1[ic21_i] * w2[ic12_i]) * w3[ic33_r]
898  - (w1[ic21_r] * w2[ic12_i] + w1[ic21_i] * w2[ic12_r]) * w3[ic33_i];
899  c_i[id] -= (w1[ic21_r] * w2[ic12_r] - w1[ic21_i] * w2[ic12_i]) * w3[ic33_i]
900  + (w1[ic21_r] * w2[ic12_i] + w1[ic21_i] * w2[ic12_r]) * w3[ic33_r];
901 
902  c_r[id] -= (w1[ic11_r] * w2[ic32_r] - w1[ic11_i] * w2[ic32_i]) * w3[ic23_r]
903  - (w1[ic11_r] * w2[ic32_i] + w1[ic11_i] * w2[ic32_r]) * w3[ic23_i];
904  c_i[id] -= (w1[ic11_r] * w2[ic32_r] - w1[ic11_i] * w2[ic32_i]) * w3[ic23_i]
905  + (w1[ic11_r] * w2[ic32_i] + w1[ic11_i] * w2[ic32_r]) * w3[ic23_r];
906  }
907  }
908  }
909  }
910 
911  corr = cmplx(0.0, 0.0);
912  for (int id = 0; id < ND; ++id) {
913  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
914  }
915 #endif // (USE_GROUP_SU3)
916 }
917 
918 
919 //====================================================================
920 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta
921 void contract_at_y(dcomplex& corr,
922  const GammaMatrix& gm_sink,
923  const Field_F& v1, const Field_F& v2,
924  const int y)
925 {
926 #if defined USE_GROUP_SU_N
927  const int NC = CommonParameters::Nc();
928  const int ND = CommonParameters::Nd();
929  const int NC2 = 2 * NC;
930  const int NCD2 = NC2 * ND;
931 #endif
932 
933  const int Nvol = v1.nvol();
934  int Nx = CommonParameters::Nx();
935  int Ny = CommonParameters::Ny();
936  int Nz = CommonParameters::Nz();
937  int Nt = CommonParameters::Nt();
938 
939  assert(Nvol == v1.nvol());
940  assert(Nvol == v2.nvol());
941  assert(v1.nex() == 1);
942  assert(v2.nex() == 1);
943  assert(y < Ny);
944 
945  const double *w1 = v1.ptr(0);
946  const double *w2 = v2.ptr(0);
947 
948  int id1[ND];
949  int id2[ND];
950  for (int id = 0; id < ND; ++id) {
951  id1[id] = id * NC2;
952  id2[id] = gm_sink.index(id) * NC2;
953  }
954 
955  double c_r[ND];
956  double c_i[ND];
957  for (int id = 0; id < ND; ++id) {
958  c_r[id] = 0.0;
959  c_i[id] = 0.0;
960  }
961 
962  for (int t = 0; t < Nt; ++t) {
963  for (int z = 0; z < Nz; ++z) {
964  for (int x = 0; x < Nx; ++x) {
965  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
966  for (int cc = 0; cc < NC; ++cc) {
967  for (int id = 0; id < ND; ++id) {
968  int ic1_r = 2 * cc + id1[id] + site;
969  int ic2_r = 2 * cc + id2[id] + site;
970 
971  int ic1_i = 2 * cc + 1 + id1[id] + site;
972  int ic2_i = 2 * cc + 1 + id2[id] + site;
973 
974  c_r[id] += w1[ic2_r] * w2[ic1_r]
975  + w1[ic2_i] * w2[ic1_i];
976 
977  c_i[id] += -w1[ic2_r] * w2[ic1_i]
978  + w1[ic2_i] * w2[ic1_r];
979  }
980  }
981  }
982  }
983  }
984 
985  corr = cmplx(0.0, 0.0);
986  for (int id = 0; id < ND; ++id) {
987  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
988  }
989 }
990 
991 
992 //====================================================================
993 void contract_at_y(std::vector<dcomplex>& corr_global,
994  const GammaMatrix& gm_sink,
995  const Field_F& v1, const Field_F& v2)
996 {
997  const int Ly = CommonParameters::Ly();
998  const int Ny = CommonParameters::Ny();
999 
1000  assert(corr_global.size() == Ly);
1001 
1002  std::vector<dcomplex> corr_local(Ny, 0.0);
1003  for (int y = 0; y < Ny; ++y) {
1004  dcomplex corr_y;
1005  contract_at_y(corr_y, gm_sink, v1, v2, y);
1006  corr_local[y] += corr_y;
1007  }
1008  global_corr_y(corr_global, corr_local);
1009 }
1010 
1011 
1012 //====================================================================
1013 void contract_at_y(dcomplex& corr,
1014  const std::vector<int>& momentum_sink,
1015  const GammaMatrix& gm_sink,
1016  const std::vector<int>& source_position,
1017  const Field_F& v1, const Field_F& v2,
1018  const int y)
1019 {
1020 #if defined USE_GROUP_SU_N
1021  int NC = CommonParameters::Nc();
1022  int ND = CommonParameters::Nd();
1023  int NC2 = 2 * NC;
1024  int NCD2 = NC2 * ND;
1025 #endif
1026 
1027  const int Nvol = v1.nvol();
1028  const int Nvol_s = Nvol / CommonParameters::Nt();
1029  const int Nx = CommonParameters::Nx();
1030  const int Ny = CommonParameters::Ny();
1031  const int Nz = CommonParameters::Nz();
1032  const int Nt = CommonParameters::Nt();
1033  const int Lx = CommonParameters::Lx();
1034  const int Ly = CommonParameters::Ly();
1035  const int Lz = CommonParameters::Lz();
1036  const int Lt = CommonParameters::Lt();
1037  const int Ndim = CommonParameters::Ndim();
1038 
1039  assert(Nvol == v2.nvol());
1040  assert(v1.nex() == 1);
1041  assert(v2.nex() == 1);
1042  assert(momentum_sink.size() == Ndim - 1);
1043  assert(source_position.size() == Ndim);
1044 
1045  const double *w1 = v1.ptr(0);
1046  const double *w2 = v2.ptr(0);
1047 
1048  int id1[ND];
1049  int id2[ND];
1050  for (int id = 0; id < ND; ++id) {
1051  id1[id] = id * NC2;
1052  id2[id] = gm_sink.index(id) * NC2;
1053  }
1054 
1055  double c_r[ND];
1056  double c_i[ND];
1057  for (int id = 0; id < ND; ++id) {
1058  c_r[id] = 0.0;
1059  c_i[id] = 0.0;
1060  }
1061 
1062  static const double PI = 4.0 * atan(1.0);
1063  std::vector<double> p_unit(ND - 1);
1064  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1065  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1066  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1067 
1068  std::vector<int> ipe(ND);
1069  ipe[0] = Communicator::ipe(0);
1070  ipe[1] = Communicator::ipe(1);
1071  ipe[2] = Communicator::ipe(2);
1072  ipe[3] = Communicator::ipe(3);
1073 
1074  for (int t = 0; t < Nt; ++t) {
1075  for (int z = 0; z < Nz; ++z) {
1076  for (int x = 0; x < Nx; ++x) {
1077  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1078 
1079  int x_global = x + ipe[0] * Nx;
1080  int z_global = z + ipe[2] * Nz;
1081  int t_global = t + ipe[3] * Nt;
1082 
1083  double p_x = p_unit[0] * (x_global - source_position[0]);
1084  double p_z = p_unit[1] * (z_global - source_position[2]);
1085  double p_t = p_unit[2] * (t_global - source_position[3]);
1086 
1087  double cos_p_xzt = cos(p_t + p_x + p_z);
1088  double sin_p_xzt = sin(p_t + p_x + p_z);
1089 
1090  for (int cc = 0; cc < NC; ++cc) {
1091  for (int id = 0; id < ND; ++id) {
1092  int ic1_r = 2 * cc + id1[id] + site;
1093  int ic2_r = 2 * cc + id2[id] + site;
1094 
1095  int ic1_i = 2 * cc + 1 + id1[id] + site;
1096  int ic2_i = 2 * cc + 1 + id2[id] + site;
1097 
1098  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1099  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1100 
1101  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1102  c_r[id] += w1_w2_r * cos_p_xzt - w1_w2_i * sin_p_xzt;
1103  c_i[id] += w1_w2_r * sin_p_xzt + w1_w2_i * cos_p_xzt;
1104  }
1105  }
1106  }
1107  }
1108  }
1109 
1110  corr = cmplx(0.0, 0.0);
1111  for (int id = 0; id < ND; ++id) {
1112  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1113  }
1114 }
1115 
1116 
1117 //====================================================================
1118 void contract_at_y(std::vector<dcomplex>& corr_global,
1119  const std::vector<int>& momentum_sink,
1120  const GammaMatrix& gm_sink,
1121  const std::vector<int>& source_position,
1122  const Field_F& v1, const Field_F& v2)
1123 {
1124  const int Ly = CommonParameters::Ly();
1125  const int Ny = CommonParameters::Ny();
1126 
1127  assert(corr_global.size() == Ly);
1128 
1129  std::vector<dcomplex> corr_local(Ny, 0.0);
1130  for (int y = 0; y < Ny; ++y) {
1131  dcomplex corr_y;
1132  contract_at_y(corr_y, momentum_sink, gm_sink, source_position, v1, v2, y);
1133  corr_local[y] += corr_y;
1134  }
1135  global_corr_y(corr_global, corr_local);
1136 }
1137 
1138 
1139 //====================================================================
1140 void contract_at_y_cos(dcomplex& corr,
1141  const std::vector<int>& momentum_sink,
1142  const GammaMatrix& gm_sink,
1143  const std::vector<int>& source_position,
1144  const Field_F& v1, const Field_F& v2,
1145  const int y)
1146 {
1147 #if defined USE_GROUP_SU_N
1148  int NC = CommonParameters::Nc();
1149  int ND = CommonParameters::Nd();
1150  int NC2 = 2 * NC;
1151  int NCD2 = NC2 * ND;
1152 #endif
1153 
1154  const int Nvol = v1.nvol();
1155  const int Nvol_s = Nvol / CommonParameters::Nt();
1156  const int Nx = CommonParameters::Nx();
1157  const int Ny = CommonParameters::Ny();
1158  const int Nz = CommonParameters::Nz();
1159  const int Nt = CommonParameters::Nt();
1160  const int Lx = CommonParameters::Lx();
1161  const int Ly = CommonParameters::Ly();
1162  const int Lz = CommonParameters::Lz();
1163  const int Lt = CommonParameters::Lt();
1164  const int Ndim = CommonParameters::Ndim();
1165 
1166  assert(Nvol == v2.nvol());
1167  assert(v1.nex() == 1);
1168  assert(v2.nex() == 1);
1169  assert(momentum_sink.size() == Ndim - 1);
1170  assert(source_position.size() == Ndim);
1171 
1172  const double *w1 = v1.ptr(0);
1173  const double *w2 = v2.ptr(0);
1174 
1175  int id1[ND];
1176  int id2[ND];
1177  for (int id = 0; id < ND; ++id) {
1178  id1[id] = id * NC2;
1179  id2[id] = gm_sink.index(id) * NC2;
1180  }
1181 
1182  double c_r[ND];
1183  double c_i[ND];
1184  for (int id = 0; id < ND; ++id) {
1185  c_r[id] = 0.0;
1186  c_i[id] = 0.0;
1187  }
1188 
1189  static const double PI = 4.0 * atan(1.0);
1190  std::vector<double> p_unit(ND - 1);
1191  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1192  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1193  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1194 
1195  std::vector<int> ipe(ND);
1196  ipe[0] = Communicator::ipe(0);
1197  ipe[1] = Communicator::ipe(1);
1198  ipe[2] = Communicator::ipe(2);
1199  ipe[3] = Communicator::ipe(3);
1200 
1201 
1202  for (int t = 0; t < Nt; ++t) {
1203  for (int z = 0; z < Nz; ++z) {
1204  for (int x = 0; x < Nx; ++x) {
1205  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1206 
1207  int x_global = x + ipe[0] * Nx;
1208  int z_global = z + ipe[2] * Nz;
1209  int t_global = t + ipe[3] * Nt;
1210 
1211  double p_x = p_unit[0] * (x_global - source_position[0]);
1212  double p_z = p_unit[1] * (z_global - source_position[2]);
1213  double p_t = p_unit[2] * (t_global - source_position[3]);
1214 
1215  double cos_p_xzt = cos(p_t + p_x + p_z);
1216  //double sin_p_xzt = sin(p_t + p_x + p_z);
1217  double sin_p_xzt = 0;
1218 
1219  for (int cc = 0; cc < NC; ++cc) {
1220  for (int id = 0; id < ND; ++id) {
1221  int ic1_r = 2 * cc + id1[id] + site;
1222  int ic2_r = 2 * cc + id2[id] + site;
1223 
1224  int ic1_i = 2 * cc + 1 + id1[id] + site;
1225  int ic2_i = 2 * cc + 1 + id2[id] + site;
1226 
1227  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1228  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1229  //double w1_w2_i = 0;
1230 
1231  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1232  c_r[id] += w1_w2_r * cos_p_xzt - w1_w2_i * sin_p_xzt;
1233  c_i[id] += w1_w2_r * sin_p_xzt + w1_w2_i * cos_p_xzt;
1234  }
1235  }
1236  }
1237  }
1238  }
1239 
1240  corr = cmplx(0.0, 0.0);
1241  for (int id = 0; id < ND; ++id) {
1242  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1243  }
1244 }
1245 
1246 
1247 //====================================================================
1248 void contract_at_y_cos(std::vector<dcomplex>& corr_global,
1249  const std::vector<int>& momentum_sink,
1250  const GammaMatrix& gm_sink,
1251  const std::vector<int>& source_position,
1252  const Field_F& v1, const Field_F& v2)
1253 {
1254  const int Ly = CommonParameters::Ly();
1255  const int Ny = CommonParameters::Ny();
1256 
1257  assert(corr_global.size() == Ly);
1258 
1259  std::vector<dcomplex> corr_local(Ny, 0.0);
1260  for (int y = 0; y < Ny; ++y) {
1261  dcomplex corr_y;
1262  contract_at_y_cos(corr_y, momentum_sink, gm_sink, source_position, v1, v2, y);
1263  corr_local[y] += corr_y;
1264  }
1265  global_corr_y(corr_global, corr_local);
1266 }
1267 
1268 
1269 //====================================================================
1270 // corr=(v2^*)_alpha (gm)_{alpha,beta} (v1)_beta
1271 void contract_at_z(dcomplex& corr,
1272  const GammaMatrix& gm_sink,
1273  const Field_F& v1, const Field_F& v2,
1274  const int z)
1275 {
1276 #if defined USE_GROUP_SU_N
1277  const int NC = CommonParameters::Nc();
1278  const int ND = CommonParameters::Nd();
1279  const int NC2 = 2 * NC;
1280  const int NCD2 = NC2 * ND;
1281 #endif
1282 
1283  int Nx = CommonParameters::Nx();
1284  int Ny = CommonParameters::Ny();
1285  int Nz = CommonParameters::Nz();
1286  int Nt = CommonParameters::Nt();
1287  const int Nvol = v1.nvol();
1288 
1289  assert(Nvol == v1.nvol());
1290  assert(Nvol == v2.nvol());
1291  assert(v1.nex() == 1);
1292  assert(v2.nex() == 1);
1293  assert(z < Nz);
1294 
1295  const double *w1 = v1.ptr(0);
1296  const double *w2 = v2.ptr(0);
1297 
1298  int id1[ND];
1299  int id2[ND];
1300  for (int id = 0; id < ND; ++id) {
1301  id1[id] = id * NC2;
1302  id2[id] = gm_sink.index(id) * NC2;
1303  }
1304 
1305  double c_r[ND];
1306  double c_i[ND];
1307  for (int id = 0; id < ND; ++id) {
1308  c_r[id] = 0.0;
1309  c_i[id] = 0.0;
1310  }
1311 
1312  for (int t = 0; t < Nt; ++t) {
1313  for (int y = 0; y < Ny; ++y) {
1314  for (int x = 0; x < Nx; ++x) {
1315  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1316  for (int cc = 0; cc < NC; ++cc) {
1317  for (int id = 0; id < ND; ++id) {
1318  int ic1_r = 2 * cc + id1[id] + site;
1319  int ic2_r = 2 * cc + id2[id] + site;
1320 
1321  int ic1_i = 2 * cc + 1 + id1[id] + site;
1322  int ic2_i = 2 * cc + 1 + id2[id] + site;
1323 
1324  c_r[id] += w1[ic2_r] * w2[ic1_r]
1325  + w1[ic2_i] * w2[ic1_i];
1326 
1327  c_i[id] += -w1[ic2_r] * w2[ic1_i]
1328  + w1[ic2_i] * w2[ic1_r];
1329  }
1330  }
1331  }
1332  }
1333  }
1334 
1335  corr = cmplx(0.0, 0.0);
1336  for (int id = 0; id < ND; ++id) {
1337  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1338  }
1339 }
1340 
1341 
1342 //====================================================================
1343 void contract_at_z(std::vector<dcomplex>& corr_global,
1344  const GammaMatrix& gm_sink,
1345  const Field_F& v1, const Field_F& v2)
1346 {
1347  const int Lz = CommonParameters::Lz();
1348  const int Nz = CommonParameters::Nz();
1349 
1350  assert(corr_global.size() == Lz);
1351 
1352  std::vector<dcomplex> corr_local(Nz, 0.0);
1353  for (int z = 0; z < Nz; ++z) {
1354  dcomplex corr_z;
1355  contract_at_z(corr_z, gm_sink, v1, v2, z);
1356  corr_local[z] += corr_z;
1357  }
1358  global_corr_z(corr_global, corr_local);
1359 }
1360 
1361 
1362 //====================================================================
1363 void contract_at_z(dcomplex& corr,
1364  const std::vector<int>& momentum_sink,
1365  const GammaMatrix& gm_sink,
1366  const std::vector<int>& source_position,
1367  const Field_F& v1, const Field_F& v2,
1368  const int z)
1369 {
1370 #if defined USE_GROUP_SU_N
1371  int NC = CommonParameters::Nc();
1372  int ND = CommonParameters::Nd();
1373  int NC2 = 2 * NC;
1374  int NCD2 = NC2 * ND;
1375 #endif
1376 
1377  const int Nvol = v1.nvol();
1378  const int Nvol_s = Nvol / CommonParameters::Nt();
1379  const int Nx = CommonParameters::Nx();
1380  const int Ny = CommonParameters::Ny();
1381  const int Nz = CommonParameters::Nz();
1382  const int Nt = CommonParameters::Nt();
1383  const int Lx = CommonParameters::Lx();
1384  const int Ly = CommonParameters::Ly();
1385  const int Lz = CommonParameters::Lz();
1386  const int Lt = CommonParameters::Lt();
1387  const int Ndim = CommonParameters::Ndim();
1388 
1389  assert(Nvol == v2.nvol());
1390  assert(v1.nex() == 1);
1391  assert(v2.nex() == 1);
1392  assert(momentum_sink.size() == Ndim - 1);
1393  assert(source_position.size() == Ndim);
1394 
1395  const double *w1 = v1.ptr(0);
1396  const double *w2 = v2.ptr(0);
1397 
1398  int id1[ND];
1399  int id2[ND];
1400  for (int id = 0; id < ND; ++id) {
1401  id1[id] = id * NC2;
1402  id2[id] = gm_sink.index(id) * NC2;
1403  }
1404 
1405  double c_r[ND];
1406  double c_i[ND];
1407  for (int id = 0; id < ND; ++id) {
1408  c_r[id] = 0.0;
1409  c_i[id] = 0.0;
1410  }
1411 
1412  static const double PI = 4.0 * atan(1.0);
1413  std::vector<double> p_unit(ND - 1);
1414  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1415  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1416  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1417 
1418  std::vector<int> ipe(ND);
1419  ipe[0] = Communicator::ipe(0);
1420  ipe[1] = Communicator::ipe(1);
1421  ipe[2] = Communicator::ipe(2);
1422  ipe[3] = Communicator::ipe(3);
1423 
1424  for (int t = 0; t < Nt; ++t) {
1425  for (int y = 0; y < Ny; ++y) {
1426  for (int x = 0; x < Nx; ++x) {
1427  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1428 
1429  int x_global = x + ipe[0] * Nx;
1430  int y_global = y + ipe[1] * Ny;
1431  int t_global = t + ipe[3] * Nt;
1432 
1433  double p_x = p_unit[0] * (x_global - source_position[0]);
1434  double p_y = p_unit[1] * (y_global - source_position[1]);
1435  double p_t = p_unit[2] * (t_global - source_position[3]);
1436 
1437  double cos_p_xyt = cos(p_t + p_x + p_y);
1438  double sin_p_xyt = sin(p_t + p_x + p_y);
1439 
1440  for (int cc = 0; cc < NC; ++cc) {
1441  for (int id = 0; id < ND; ++id) {
1442  int ic1_r = 2 * cc + id1[id] + site;
1443  int ic2_r = 2 * cc + id2[id] + site;
1444 
1445  int ic1_i = 2 * cc + 1 + id1[id] + site;
1446  int ic2_i = 2 * cc + 1 + id2[id] + site;
1447 
1448  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1449  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1450 
1451  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1452  c_r[id] += w1_w2_r * cos_p_xyt - w1_w2_i * sin_p_xyt;
1453  c_i[id] += w1_w2_r * sin_p_xyt + w1_w2_i * cos_p_xyt;
1454  }
1455  }
1456  }
1457  }
1458  }
1459 
1460  corr = cmplx(0.0, 0.0);
1461  for (int id = 0; id < ND; ++id) {
1462  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1463  }
1464 }
1465 
1466 
1467 //====================================================================
1468 void contract_at_z(std::vector<dcomplex>& corr_global,
1469  const std::vector<int>& momentum_sink,
1470  const GammaMatrix& gm_sink,
1471  const std::vector<int>& source_position,
1472  const Field_F& v1, const Field_F& v2)
1473 {
1474  const int Lz = CommonParameters::Lz();
1475  const int Nz = CommonParameters::Nz();
1476 
1477  assert(corr_global.size() == Lz);
1478 
1479  std::vector<dcomplex> corr_local(Nz, 0.0);
1480  for (int z = 0; z < Nz; ++z) {
1481  dcomplex corr_z;
1482  contract_at_z(corr_z, momentum_sink, gm_sink, source_position, v1, v2, z);
1483  corr_local[z] += corr_z;
1484  }
1485  global_corr_z(corr_global, corr_local);
1486 }
1487 
1488 
1489 //====================================================================
1490 void contract_at_z_cos(dcomplex& corr,
1491  const std::vector<int>& momentum_sink,
1492  const GammaMatrix& gm_sink,
1493  const std::vector<int>& source_position,
1494  const Field_F& v1, const Field_F& v2,
1495  const int z)
1496 {
1497 #if defined USE_GROUP_SU_N
1498  int NC = CommonParameters::Nc();
1499  int ND = CommonParameters::Nd();
1500  int NC2 = 2 * NC;
1501  int NCD2 = NC2 * ND;
1502 #endif
1503 
1504  const int Nvol = v1.nvol();
1505  const int Nvol_s = Nvol / CommonParameters::Nt();
1506  const int Nx = CommonParameters::Nx();
1507  const int Ny = CommonParameters::Ny();
1508  const int Nz = CommonParameters::Nz();
1509  const int Nt = CommonParameters::Nt();
1510  const int Lx = CommonParameters::Lx();
1511  const int Ly = CommonParameters::Ly();
1512  const int Lz = CommonParameters::Lz();
1513  const int Lt = CommonParameters::Lt();
1514  const int Ndim = CommonParameters::Ndim();
1515 
1516  assert(Nvol == v2.nvol());
1517  assert(v1.nex() == 1);
1518  assert(v2.nex() == 1);
1519  assert(momentum_sink.size() == Ndim - 1);
1520  assert(source_position.size() == Ndim);
1521 
1522  const double *w1 = v1.ptr(0);
1523  const double *w2 = v2.ptr(0);
1524 
1525  int id1[ND];
1526  int id2[ND];
1527  for (int id = 0; id < ND; ++id) {
1528  id1[id] = id * NC2;
1529  id2[id] = gm_sink.index(id) * NC2;
1530  }
1531 
1532  double c_r[ND];
1533  double c_i[ND];
1534  for (int id = 0; id < ND; ++id) {
1535  c_r[id] = 0.0;
1536  c_i[id] = 0.0;
1537  }
1538 
1539  static const double PI = 4.0 * atan(1.0);
1540  std::vector<double> p_unit(ND - 1);
1541  p_unit[0] = (2.0 * PI / Ly) * momentum_sink[0];
1542  p_unit[1] = (2.0 * PI / Lz) * momentum_sink[1];
1543  p_unit[2] = (2.0 * PI / Lt) * momentum_sink[2];
1544 
1545  std::vector<int> ipe(ND);
1546  ipe[0] = Communicator::ipe(0);
1547  ipe[1] = Communicator::ipe(1);
1548  ipe[2] = Communicator::ipe(2);
1549  ipe[3] = Communicator::ipe(3);
1550 
1551 
1552  for (int t = 0; t < Nt; ++t) {
1553  for (int y = 0; y < Ny; ++y) {
1554  for (int x = 0; x < Nx; ++x) {
1555  int site = NCD2 * (x + Nx * (y + Ny * (z + Nz * t)));
1556 
1557  int x_global = x + ipe[0] * Nx;
1558  int y_global = y + ipe[1] * Ny;
1559  int t_global = t + ipe[3] * Nt;
1560 
1561  double p_x = p_unit[0] * (x_global - source_position[0]);
1562  double p_y = p_unit[1] * (y_global - source_position[1]);
1563  double p_t = p_unit[2] * (t_global - source_position[2]);
1564 
1565  double cos_p_xyt = cos(p_t + p_x + p_y);
1566  //double sin_p_xyt = sin(p_t + p_x + p_y);
1567  double sin_p_xyt = 0;
1568 
1569  for (int cc = 0; cc < NC; ++cc) {
1570  for (int id = 0; id < ND; ++id) {
1571  int ic1_r = 2 * cc + id1[id] + site;
1572  int ic2_r = 2 * cc + id2[id] + site;
1573 
1574  int ic1_i = 2 * cc + 1 + id1[id] + site;
1575  int ic2_i = 2 * cc + 1 + id2[id] + site;
1576 
1577  double w1_w2_r = w1[ic2_r] * w2[ic1_r] + w1[ic2_i] * w2[ic1_i];
1578  double w1_w2_i = -w1[ic2_r] * w2[ic1_i] + w1[ic2_i] * w2[ic1_r];
1579  //double w1_w2_i = 0;
1580 
1581  //- corr[s0] += v2^dagger * v1 * exp(i * p_i * x_i)
1582  c_r[id] += w1_w2_r * cos_p_xyt - w1_w2_i * sin_p_xyt;
1583  c_i[id] += w1_w2_r * sin_p_xyt + w1_w2_i * cos_p_xyt;
1584  }
1585  }
1586  }
1587  }
1588  }
1589 
1590  corr = cmplx(0.0, 0.0);
1591  for (int id = 0; id < ND; ++id) {
1592  corr += gm_sink.value(id) * cmplx(c_r[id], c_i[id]);
1593  }
1594 }
1595 
1596 
1597 //====================================================================
1598 void contract_at_z_cos(std::vector<dcomplex>& corr_global,
1599  const std::vector<int>& momentum_sink,
1600  const GammaMatrix& gm_sink,
1601  const std::vector<int>& source_position,
1602  const Field_F& v1, const Field_F& v2)
1603 {
1604  const int Lz = CommonParameters::Lz();
1605  const int Nz = CommonParameters::Nz();
1606 
1607  assert(corr_global.size() == Lz);
1608 
1609  std::vector<dcomplex> corr_local(Nz, 0.0);
1610  for (int z = 0; z < Nz; ++z) {
1611  dcomplex corr_z;
1612  contract_at_z_cos(corr_z, momentum_sink, gm_sink, source_position, v1, v2, z);
1613  corr_local[z] += corr_z;
1614  }
1615  global_corr_z(corr_global, corr_local);
1616 }
1617 
1618 
1619 //====================================================================
1620 void global_corr_x(std::vector<dcomplex>& corr_global,
1621  std::vector<dcomplex>& corr_local)
1622 {
1623  int Lx = CommonParameters::Lx();
1624  int Nx = CommonParameters::Nx();
1625 
1626  assert(corr_global.size() == Lx);
1627  assert(corr_local.size() == Nx);
1628 
1629  std::vector<dcomplex> corr_tmp(Lx, 0);
1630 
1631  int ipex = Communicator::ipe(0);
1632 
1633  for (int x = 0; x < Nx; ++x) {
1634  int x2 = x + ipex * Nx;
1635  corr_tmp[x2] = corr_local[x];
1636  }
1637 
1638  for (int x = 0; x < Lx; ++x) {
1639  double crr = Communicator::reduce_sum(real(corr_tmp[x]));
1640  double cri = Communicator::reduce_sum(imag(corr_tmp[x]));
1641  corr_global[x] = cmplx(crr, cri);
1642  }
1643 }
1644 
1645 
1646 //====================================================================
1647 void global_corr_y(std::vector<dcomplex>& corr_global,
1648  std::vector<dcomplex>& corr_local)
1649 {
1650  int Ly = CommonParameters::Ly();
1651  int Ny = CommonParameters::Ny();
1652 
1653  assert(corr_global.size() == Ly);
1654  assert(corr_local.size() == Ny);
1655 
1656  std::vector<dcomplex> corr_tmp(Ly, 0);
1657 
1658  int ipey = Communicator::ipe(1);
1659 
1660  for (int y = 0; y < Ny; ++y) {
1661  int y2 = y + ipey * Ny;
1662  corr_tmp[y2] = corr_local[y];
1663  }
1664 
1665  for (int y = 0; y < Ly; ++y) {
1666  double crr = Communicator::reduce_sum(real(corr_tmp[y]));
1667  double cri = Communicator::reduce_sum(imag(corr_tmp[y]));
1668  corr_global[y] = cmplx(crr, cri);
1669  }
1670 }
1671 
1672 
1673 //====================================================================
1674 void global_corr_z(std::vector<dcomplex>& corr_global,
1675  std::vector<dcomplex>& corr_local)
1676 {
1677  int Lz = CommonParameters::Lz();
1678  int Nz = CommonParameters::Nz();
1679 
1680  assert(corr_global.size() == Lz);
1681  assert(corr_local.size() == Nz);
1682 
1683  std::vector<dcomplex> corr_tmp(Lz);
1684 
1685  int ipez = Communicator::ipe(2);
1686 
1687  for (int z = 0; z < Nz; ++z) {
1688  int z2 = z + ipez * Nz;
1689  corr_tmp[z2] = corr_local[z];
1690  }
1691 
1692  for (int z = 0; z < Lz; ++z) {
1693  double crr = Communicator::reduce_sum(real(corr_tmp[z]));
1694  double cri = Communicator::reduce_sum(imag(corr_tmp[z]));
1695  corr_global[z] = cmplx(crr, cri);
1696  }
1697 }
1698 
1699 
1700 //====================================================================
1701 void global_corr_t(std::vector<dcomplex>& corr_global,
1702  std::vector<dcomplex>& corr_local)
1703 {
1704  int Lt = CommonParameters::Lt();
1705  int Nt = CommonParameters::Nt();
1706 
1707  assert(corr_global.size() == Lt);
1708  assert(corr_local.size() == Nt);
1709 
1710  std::vector<dcomplex> corr_tmp(Lt, 0);
1711 
1712  int ipe_t = Communicator::ipe(3);
1713 
1714  for (int t = 0; t < Nt; ++t) {
1715  int t_global = t + ipe_t * Nt;
1716  corr_tmp[t_global] = corr_local[t];
1717  }
1718 
1719  for (int t_global = 0; t_global < Lt; ++t_global) {
1720  double cr_r = Communicator::reduce_sum(real(corr_tmp[t_global]));
1721  double cr_i = Communicator::reduce_sum(imag(corr_tmp[t_global]));
1722  corr_global[t_global] = cmplx(cr_r, cr_i);
1723  }
1724 }
1725 
1726 
1727 //====================================================================
1728 //============================================================END=====
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
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
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
NC2
#define NC2
Definition: field_F_imp_SU2-inc.h:3
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
Index_lex
Lexical site index.
Definition: index_lex.h:34
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83
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
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
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
Field::nex
int nex() const
Definition: field.h:128
CommonParameters::Ly
static int Ly()
Definition: commonParameters.h:92
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
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
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
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
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
NC
#define NC
Definition: field_F_imp_SU2-inc.h:2
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
CommonParameters::Lz
static int Lz()
Definition: commonParameters.h:93
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
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
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
NCD2
#define NCD2
Definition: field_F_imp_SU2-inc.h:7
ND
#define ND
Definition: field_F_imp_SU2-inc.h:5
Field::nvol
int nvol() const
Definition: field.h:127
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
index_lex.h
contract_4spinor.h
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:55
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
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
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
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
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
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
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
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