Bridge++  Ver. 1.3.x
field_F_imp.cpp
Go to the documentation of this file.
1 
14 #include "field_F.h"
15 
16 using std::valarray;
17 
18 //#include "threadManager_OpenMP.h"
19 
20 #if defined USE_GROUP_SU3
21 #include "field_F_imp_SU3.inc"
22 #elif defined USE_GROUP_SU2
23 #include "field_F_imp_SU2.inc"
24 #elif defined USE_GROUP_SU_N
25 #include "field_F_imp_SU_N.inc"
26 #endif
27 
28 //====================================================================
29 void Field_F::check()
30 {
31  // assert(NC == CommonParameters::Nc());
32  // assert(ND == CommonParameters::Nd());
33 
34  check_Nc();
35 }
36 
37 
38 //====================================================================
39 void mult_Field_Gn(Field_F& y, const int ex,
40  const Field_G& U, const int ex1,
41  const Field_F& x, const int ex2)
42 {
43  assert(ex < y.nex());
44  assert(ex1 < U.nex());
45  assert(ex2 < x.nex());
46  assert(U.nvol() == y.nvol());
47  assert(x.nvol() == y.nvol());
48 
49  const int Nd = x.nd();
50  const int Nc = x.nc();
51  const int Nc2 = x.nc2();
52  const int Ncd2 = Nc2 * Nd;
53  const int Ndf = Nc2 * Nc;
54 
55  // int Nthread = ThreadManager_OpenMP::get_num_threads();
56  // int i_thread = ThreadManager_OpenMP::get_thread_id();
57  int Nthread = 1;
58  int i_thread = 0;
59 
60  int is = y.nvol() * i_thread / Nthread;
61  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
62 
63  double *v = y.ptr(0, is, ex);
64  const double *g = U.ptr(0, is, ex1);
65  const double *w = x.ptr(0, is, ex2);
66 
67  for (int site = 0; site < ns; ++site) {
68  int ig = Ndf * site;
69  int iv = Ncd2 * site;
70  for (int s = 0; s < Nd; ++s) {
71  for (int ic = 0; ic < Nc; ++ic) {
72  int ig2 = ic * Nc2 + ig;
73  int iv2 = s * Nc2 + iv;
74  v[2 * ic + iv2] = mult_Gn_r(&g[ig2], &w[iv2], Nc);
75  v[2 * ic + 1 + iv2] = mult_Gn_i(&g[ig2], &w[iv2], Nc);
76  }
77  }
78  }
79 }
80 
81 
82 //====================================================================
83 void mult_Field_Gd(Field_F& y, const int ex,
84  const Field_G& U, const int ex1,
85  const Field_F& x, const int ex2)
86 {
87  assert(ex < y.nex());
88  assert(ex1 < U.nex());
89  assert(ex2 < x.nex());
90  assert(U.nvol() == y.nvol());
91  assert(x.nvol() == y.nvol());
92 
93  const int Nd = x.nd();
94  const int Nc = x.nc();
95  const int Nc2 = x.nc2();
96  const int Ncd2 = Nc2 * Nd;
97  const int Ndf = Nc2 * Nc;
98 
99  // int Nthread = ThreadManager_OpenMP::get_num_threads();
100  // int i_thread = ThreadManager_OpenMP::get_thread_id();
101  int Nthread = 1;
102  int i_thread = 0;
103 
104  int is = y.nvol() * i_thread / Nthread;
105  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
106 
107  double *v = y.ptr(0, is, ex);
108  const double *g = U.ptr(0, is, ex1);
109  const double *w = x.ptr(0, is, ex2);
110 
111  for (int site = 0; site < ns; ++site) {
112  int ig = Ndf * site;
113  int iv = Ncd2 * site;
114  for (int s = 0; s < Nd; ++s) {
115  for (int ic = 0; ic < Nc; ++ic) {
116  int ig2 = ic * 2 + ig;
117  int iv2 = s * Nc2 + iv;
118  v[2 * ic + iv2] = mult_Gd_r(&g[ig2], &w[iv2], Nc);
119  v[2 * ic + 1 + iv2] = mult_Gd_i(&g[ig2], &w[iv2], Nc);
120  }
121  }
122  }
123 }
124 
125 
126 //====================================================================
127 void multadd_Field_Gn(Field_F& y, const int ex,
128  const Field_G& U, const int ex1,
129  const Field_F& x, const int ex2,
130  const double a)
131 {
132  assert(ex < y.nex());
133  assert(ex1 < U.nex());
134  assert(ex2 < x.nex());
135  assert(U.nvol() == y.nvol());
136  assert(x.nvol() == y.nvol());
137 
138  const int Nd = x.nd();
139  const int Nc = x.nc();
140  const int Nc2 = x.nc2();
141  const int Ncd2 = Nc2 * Nd;
142  const int Ndf = Nc2 * Nc;
143 
144  // int Nthread = ThreadManager_OpenMP::get_num_threads();
145  // int i_thread = ThreadManager_OpenMP::get_thread_id();
146  int Nthread = 1;
147  int i_thread = 0;
148 
149  int is = y.nvol() * i_thread / Nthread;
150  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
151 
152  double *v = y.ptr(0, is, ex);
153  const double *g = U.ptr(0, is, ex1);
154  const double *w = x.ptr(0, is, ex2);
155 
156  for (int site = 0; site < ns; ++site) {
157  int ig = Ndf * site;
158  int iv = Ncd2 * site;
159  for (int s = 0; s < Nd; ++s) {
160  for (int ic = 0; ic < Nc; ++ic) {
161  int ig2 = ic * Nc2 + ig;
162  int iv2 = s * Nc2 + iv;
163  v[2 * ic + iv2] += a * mult_Gn_r(&g[ig2], &w[iv2], Nc);
164  v[2 * ic + 1 + iv2] += a * mult_Gn_i(&g[ig2], &w[iv2], Nc);
165  }
166  }
167  }
168 }
169 
170 
171 //====================================================================
172 void multadd_Field_Gd(Field_F& y, const int ex,
173  const Field_G& U, const int ex1,
174  const Field_F& x, const int ex2,
175  const double a)
176 {
177  assert(ex < y.nex());
178  assert(ex1 < U.nex());
179  assert(ex2 < x.nex());
180  assert(U.nvol() == y.nvol());
181  assert(x.nvol() == y.nvol());
182 
183  const int Nd = x.nd();
184  const int Nc = x.nc();
185  const int Nc2 = x.nc2();
186  const int Ncd2 = Nc2 * Nd;
187  const int Ndf = Nc2 * Nc;
188 
189  // int Nthread = ThreadManager_OpenMP::get_num_threads();
190  // int i_thread = ThreadManager_OpenMP::get_thread_id();
191  int Nthread = 1;
192  int i_thread = 0;
193 
194  int is = y.nvol() * i_thread / Nthread;
195  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
196 
197  double *v = y.ptr(0, is, ex);
198  const double *g = U.ptr(0, is, ex1);
199  const double *w = x.ptr(0, is, ex2);
200 
201  for (int site = 0; site < ns; ++site) {
202  int ig = Ndf * site;
203  int iv = Ncd2 * site;
204  for (int s = 0; s < Nd; ++s) {
205  for (int ic = 0; ic < Nc; ++ic) {
206  int ig2 = ic * 2 + ig;
207  int iv2 = s * Nc2 + iv;
208  v[2 * ic + iv2] += a * mult_Gd_r(&g[ig2], &w[iv2], Nc);
209  v[2 * ic + 1 + iv2] += a * mult_Gd_i(&g[ig2], &w[iv2], Nc);
210  }
211  }
212  }
213 }
214 
215 
216 //====================================================================
217 void mult_GM(Field_F& y, const GammaMatrix& gm, const Field_F& x)
218 {
219  assert(x.nex() == y.nex());
220  assert(x.nvol() == y.nvol());
221 
222  const int Nd = x.nd();
223  const int Nc = x.nc();
224  const int Nc2 = x.nc2();
225  const int Ncd2 = Nc2 * Nd;
226 
227  int id[Nd];
228  int idc_r[Nd];
229  int idc_i[Nd];
230  double gv_r[Nd];
231  double gv_i[Nd];
232 
233  for (int s = 0; s < Nd; ++s) {
234  id[s] = gm.index(s);
235  gv_r[s] = gm.value_r(s);
236  gv_i[s] = gm.value_i(s);
237  idc_r[s] = gm.index_c(s);
238  idc_i[s] = 1 - idc_r[s];
239  }
240 
241  // int Nthread = ThreadManager_OpenMP::get_num_threads();
242  // int i_thread = ThreadManager_OpenMP::get_thread_id();
243  int Nthread = 1;
244  int i_thread = 0;
245 
246  int is = y.nvol() * i_thread / Nthread;
247  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
248 
249  const int Nex = y.nex();
250 
251  for (int ex = 0; ex < Nex; ++ex) {
252  double *v = y.ptr(0, is, ex);
253  const double *w = x.ptr(0, is, ex);
254 
255  for (int site = 0; site < ns; ++site) {
256  int iv = Ncd2 * site;
257 
258  for (int s = 0; s < Nd; ++s) {
259  int iv2 = s * Nc2 + iv;
260  int iw2 = id[s] * Nc2 + iv;
261 
262  for (int ic = 0; ic < Nc; ++ic) {
263  v[2 * ic + iv2] = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
264  v[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
265  }
266  }
267  }
268  }
269 }
270 
271 
272 //====================================================================
273 void mult_iGM(Field_F& y, const GammaMatrix& gm, const Field_F& x)
274 {
275  assert(x.nex() == y.nex());
276  assert(x.nvol() == y.nvol());
277 
278  const int Nd = x.nd();
279  const int Nc = x.nc();
280  const int Nc2 = x.nc2();
281  const int Ncd2 = Nc2 * Nd;
282 
283  int id[Nd];
284  int idc_r[Nd];
285  int idc_i[Nd];
286  double gv_r[Nd];
287  double gv_i[Nd];
288 
289  for (int s = 0; s < Nd; ++s) {
290  id[s] = gm.index(s);
291  gv_r[s] = -gm.value_i(s);
292  gv_i[s] = gm.value_r(s);
293  idc_r[s] = 1 - gm.index_c(s);
294  idc_i[s] = 1 - idc_r[s];
295  }
296 
297  // int Nthread = ThreadManager_OpenMP::get_num_threads();
298  // int i_thread = ThreadManager_OpenMP::get_thread_id();
299  int Nthread = 1;
300  int i_thread = 0;
301 
302  int is = y.nvol() * i_thread / Nthread;
303  int ns = y.nvol() * (i_thread + 1) / Nthread - is;
304 
305  const int Nex = y.nex();
306 
307  for (int ex = 0; ex < Nex; ++ex) {
308  double *v = y.ptr(0, is, ex);
309  const double *w = x.ptr(0, is, ex);
310 
311  for (int site = 0; site < ns; ++site) {
312  int iv = Ncd2 * site;
313 
314  for (int s = 0; s < Nd; ++s) {
315  int iv2 = s * Nc2 + iv;
316  int iw2 = id[s] * Nc2 + iv;
317  for (int ic = 0; ic < Nc; ++ic) {
318  v[2 * ic + iv2] = gv_r[s] * w[2 * ic + idc_r[s] + iw2];
319  v[2 * ic + 1 + iv2] = gv_i[s] * w[2 * ic + idc_i[s] + iw2];
320  }
321  }
322  }
323  }
324 }
325 
326 
327 //====================================================================
328 void mult_GMproj(Field_F& y,
329  const int pm, const GammaMatrix& gm,
330  const Field_F& w)
331 {
332  assert(w.nvol() == y.nvol());
333  assert(w.nex() == y.nex());
334 
335  mult_GM(y, gm, w);
336 
337  if (pm == 1) {
338  axpy(y, 1.0, w); // y += w;
339  scal(y, 0.5); // y *= 0.5;
340  } else if (pm == -1) {
341  axpy(y, -1.0, w); // y -= w i.e. y = (gamma - 1) * w
342  scal(y, -0.5); // y *= -0.5 i.e. y = (1 - gamma)/2 * w
343  } else {
344  vout.crucial("%s: wrong pm.\n", __func__);
345  exit(EXIT_FAILURE);
346  }
347 }
348 
349 
350 //====================================================================
351 void mult_GMproj2(Field_F& y,
352  const int pm, const GammaMatrix& gm,
353  const Field_F& w)
354 {
355  assert(w.nvol() == y.nvol());
356  assert(w.nex() == y.nex());
357 
358  mult_GM(y, gm, w);
359 
360  if (pm == 1) {
361  axpy(y, 1.0, w); // y += w;
362  } else if (pm == -1) {
363  axpy(y, -1.0, w); // y -= w;
364  scal(y, -1.0); // y *= -1.0;
365  } else {
366  vout.crucial("%s: wrong pm.\n", __func__);
367  exit(EXIT_FAILURE);
368  }
369 }
370 
371 
372 //====================================================================
373 void mult_GMproj2(Field_F& y,
374  const double nu_s,
375  const int pm,
376  const double r_s,
377  const GammaMatrix& gm,
378  const Field_F& w)
379 {
380  assert(w.nvol() == y.nvol());
381  assert(w.nex() == y.nex());
382 
383  mult_GM(y, gm, w);
384  scal(y, nu_s); // y *= nu_s;
385 
386  if (pm == 1) {
387  axpy(y, r_s, w); // y += r_s * w;
388  } else if (pm == -1) {
389  axpy(y, -r_s, w); // y -= r_s * w i.e. y = (nu_s * gamma - r_s) * w
390  scal(y, -1.0); // y *= -1.0 i.e. y = (r_s - nu_s * gamma) * w
391  } else {
392  vout.crucial("%s: wrong pm = %d\n", __func__, pm);
393  exit(EXIT_FAILURE);
394  }
395 }
396 
397 
398 //====================================================================
399 //============================================================END=====
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:282
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:278
mult_iGM
void mult_iGM(Field_F &y, const GammaMatrix &gm, const Field_F &x)
gamma matrix multiplication (i is multiplied)
Definition: field_F_imp.cpp:273
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:133
Field_F::check
void check()
check several assumptions for performance implementation.
Definition: field_F_imp.cpp:29
Field::nvol
int nvol() const
Definition: field.h:116
mult_GM
void mult_GM(Field_F &y, const GammaMatrix &gm, const Field_F &x)
gamma matrix multiplication
Definition: field_F_imp.cpp:217
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
GammaMatrix::value_r
double value_r(int row) const
Definition: gammaMatrix.h:98
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Field::nex
int nex() const
Definition: field.h:117
field_F.h
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:168
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:48
multadd_Field_Gn
void multadd_Field_Gn(Field_F &y, const int ex, const Field_G &U, const int ex1, const Field_F &x, const int ex2, const double a)
Definition: field_F_imp.cpp:127
mult_Field_Gd
void mult_Field_Gd(Field_F &y, const int ex, const Field_G &U, const int ex1, const Field_F &x, const int ex2)
Definition: field_F_imp.cpp:83
mult_GMproj2
void mult_GMproj2(Field_F &y, const int pm, const GammaMatrix &gm, const Field_F &w)
projection with gamma matrix: (1 gamma)
Definition: field_F_imp.cpp:351
GammaMatrix::index_c
int index_c(int row) const
Definition: gammaMatrix.h:93
Field_F::nc
int nc() const
Definition: field_F.h:90
Field_F::nc2
int nc2() const
Definition: field_F.h:91
GammaMatrix::value_i
double value_i(int row) const
Definition: gammaMatrix.h:103
mult_GMproj
void mult_GMproj(Field_F &y, const int pm, const GammaMatrix &gm, const Field_F &w)
projection with gamma matrix: (1 gamma)/2
Definition: field_F_imp.cpp:328
Field_F::nd
int nd() const
Definition: field_F.h:92
mult_Field_Gn
void mult_Field_Gn(Field_F &y, const int ex, const Field_G &U, const int ex1, const Field_F &x, const int ex2)
Definition: field_F_imp.cpp:39
multadd_Field_Gd
void multadd_Field_Gd(Field_F &y, const int ex, const Field_G &U, const int ex1, const Field_F &x, const int ex2, const double a)
Definition: field_F_imp.cpp:172
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83