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