Bridge++  Ver. 1.2.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
field_G_imp.cpp
Go to the documentation of this file.
1 
14 #include "field_G.h"
15 
16 #include "bridgeIO.h"
17 using Bridge::vout;
18 
19 //#include "threadManager_OpenMP.h"
20 
21 #if defined USE_GROUP_SU3
22 #include "field_G_imp_SU3.inc"
23 #elif defined USE_GROUP_SU2
24 #include "field_G_imp_SU2.inc"
25 #elif defined USE_GROUP_SU_N
26 #include "field_G_imp_SU_N.inc"
27 #endif
28 
29 //====================================================================
30 void Field_G::check()
31 {
32  assert(NC == CommonParameters::Nc());
33 
34  check_Nc();
35 }
36 
37 
38 //====================================================================
39 void Field_G::set_unit()
40 {
41  Mat_SU_N ut(nc());
42 
43  ut.unit();
44 
45  // int nth = ThreadManager_OpenMP::get_num_threads();
46  // int ith = ThreadManager_OpenMP::get_thread_id();
47  int nth = 1;
48  int ith = 0;
49 
50  int is = nvol() * ith / nth;
51  int ns = nvol() * (ith + 1) / nth - is;
52 
53  for (int mu = 0; mu < nex(); ++mu) {
54  for (int site = is; site < is + ns; ++site) {
55  set_mat(site, mu, ut);
56  }
57  }
58 }
59 
60 
61 //====================================================================
62 void Field_G::set_random(RandomNumbers *rand)
63 {
64  rand->gauss_lex_global(*this);
65 
66  Mat_SU_N ut(m_Nc);
67 
68  for (int mu = 0; mu < m_Nex; ++mu) {
69  for (int site = 0; site < m_Nvol; ++site) {
70  this->mat(ut, site, mu);
71  ut.reunit();
72  this->set_mat(site, mu, ut);
73  }
74  }
75 }
76 
77 
78 //====================================================================
79 void Field_G::reunit()
80 {
81  Mat_SU_N ut(nc());
82 
83  // int nth = ThreadManager_OpenMP::get_num_threads();
84  // int ith = ThreadManager_OpenMP::get_thread_id();
85  int nth = 1;
86  int ith = 0;
87 
88  int is = nvol() * ith / nth;
89  int ns = nvol() * (ith + 1) / nth - is;
90 
91  for (int mu = 0; mu < nex(); ++mu) {
92  for (int site = is; site < is + ns; ++site) {
93  mat(ut, site, mu);
94  ut.reunit();
95  set_mat(site, mu, ut);
96  }
97  }
98 }
99 
100 
101 //====================================================================
102 void mult_Field_Gnn(Field_G& w, const int ex,
103  const Field_G& u1, int ex1,
104  const Field_G& u2, int ex2)
105 {
106  assert(ex < w.nex());
107  assert(ex1 < u1.nex());
108  assert(ex2 < u2.nex());
109  assert(u1.nvol() == w.nvol());
110  assert(u2.nvol() == w.nvol());
111 
112  // int nth = ThreadManager_OpenMP::get_num_threads();
113  // int ith = ThreadManager_OpenMP::get_thread_id();
114  int nth = 1;
115  int ith = 0;
116 
117  int is = w.nvol() * ith / nth;
118  int ns = w.nvol() * (ith + 1) / nth - is;
119 
120  double *g = w.ptr(0, is, ex);
121  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
122  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
123 
124  const int Nc = w.nc();
125  const int Nc2 = 2 * Nc;
126  const int Ndf = Nc2 * Nc;
127 
128  for (int site = 0; site < ns; ++site) {
129  int ig = Ndf * site;
130 
131  for (int ic2 = 0; ic2 < Nc; ++ic2) {
132  for (int ic1 = 0; ic1 < Nc; ++ic1) {
133  int jg2 = ic2 * 2 + ig;
134  int jg1 = ic1 * Nc2 + ig;
135  g[2 * ic2 + Nc2 * ic1 + ig] = mult_Gnn_r(&g1[jg1], &g2[jg2], Nc);
136  g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gnn_i(&g1[jg1], &g2[jg2], Nc);
137  }
138  }
139  }
140 }
141 
142 
143 //====================================================================
144 void mult_Field_Gdn(Field_G& w, const int ex,
145  const Field_G& u1, int ex1,
146  const Field_G& u2, int ex2)
147 {
148  assert(ex < w.nex());
149  assert(ex1 < u1.nex());
150  assert(ex2 < u2.nex());
151  assert(u1.nvol() == w.nvol());
152  assert(u2.nvol() == w.nvol());
153 
154  // int nth = ThreadManager_OpenMP::get_num_threads();
155  // int ith = ThreadManager_OpenMP::get_thread_id();
156  int nth = 1;
157  int ith = 0;
158 
159  int is = w.nvol() * ith / nth;
160  int ns = w.nvol() * (ith + 1) / nth - is;
161 
162  double *g = w.ptr(0, is, ex);
163  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
164  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
165 
166  const int Nc = w.nc();
167  const int Nc2 = 2 * Nc;
168  const int Ndf = Nc2 * Nc;
169 
170  for (int site = 0; site < ns; ++site) {
171  int ig = Ndf * site;
172 
173  for (int ic2 = 0; ic2 < Nc; ++ic2) {
174  for (int ic1 = 0; ic1 < Nc; ++ic1) {
175  int jg2 = ic2 * 2 + ig;
176  int jg1 = ic1 * 2 + ig;
177  g[2 * ic2 + Nc2 * ic1 + ig] = mult_Gdn_r(&g1[jg1], &g2[jg2], Nc);
178  g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gdn_i(&g1[jg1], &g2[jg2], Nc);
179  }
180  }
181  }
182 }
183 
184 
185 //====================================================================
186 void mult_Field_Gnd(Field_G& w, const int ex,
187  const Field_G& u1, int ex1,
188  const Field_G& u2, int ex2)
189 {
190  assert(ex < w.nex());
191  assert(ex1 < u1.nex());
192  assert(ex2 < u2.nex());
193  assert(u1.nvol() == w.nvol());
194  assert(u2.nvol() == w.nvol());
195 
196  // int nth = ThreadManager_OpenMP::get_num_threads();
197  // int ith = ThreadManager_OpenMP::get_thread_id();
198  int nth = 1;
199  int ith = 0;
200 
201  int is = w.nvol() * ith / nth;
202  int ns = w.nvol() * (ith + 1) / nth - is;
203 
204  double *g = w.ptr(0, is, ex);
205  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
206  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
207 
208  const int Nc = w.nc();
209  const int Nc2 = 2 * Nc;
210  const int Ndf = Nc2 * Nc;
211 
212  for (int site = 0; site < ns; ++site) {
213  int ig = Ndf * site;
214 
215  for (int ic2 = 0; ic2 < Nc; ++ic2) {
216  for (int ic1 = 0; ic1 < Nc; ++ic1) {
217  int jg2 = ic2 * Nc2 + ig;
218  int jg1 = ic1 * Nc2 + ig;
219  g[2 * ic2 + Nc2 * ic1 + ig] = mult_Gnd_r(&g1[jg1], &g2[jg2], Nc);
220  g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gnd_i(&g1[jg1], &g2[jg2], Nc);
221  }
222  }
223  }
224 }
225 
226 
227 //====================================================================
228 void mult_Field_Gdd(Field_G& w, const int ex,
229  const Field_G& u1, int ex1,
230  const Field_G& u2, int ex2)
231 {
232  assert(ex < w.nex());
233  assert(ex1 < u1.nex());
234  assert(ex2 < u2.nex());
235  assert(u1.nvol() == w.nvol());
236  assert(u2.nvol() == w.nvol());
237 
238  // int nth = ThreadManager_OpenMP::get_num_threads();
239  // int ith = ThreadManager_OpenMP::get_thread_id();
240  int nth = 1;
241  int ith = 0;
242 
243  int is = w.nvol() * ith / nth;
244  int ns = w.nvol() * (ith + 1) / nth - is;
245 
246  double *g = w.ptr(0, is, ex);
247  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
248  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
249 
250  const int Nc = w.nc();
251  const int Nc2 = 2 * Nc;
252  const int Ndf = Nc2 * Nc;
253 
254  for (int site = 0; site < ns; ++site) {
255  int ig = Ndf * site;
256 
257  for (int ic2 = 0; ic2 < Nc; ++ic2) {
258  for (int ic1 = 0; ic1 < Nc; ++ic1) {
259  int jg2 = ic2 * Nc2 + ig;
260  int jg1 = ic1 * 2 + ig;
261  g[2 * ic2 + Nc2 * ic1 + ig] = mult_Gdd_r(&g1[jg1], &g2[jg2], Nc);
262  g[2 * ic2 + 1 + Nc2 * ic1 + ig] = mult_Gdd_i(&g1[jg1], &g2[jg2], Nc);
263  }
264  }
265  }
266 }
267 
268 
269 //====================================================================
270 void multadd_Field_Gnn(Field_G& w, const int ex,
271  const Field_G& u1, int ex1,
272  const Field_G& u2, int ex2,
273  const double ff)
274 {
275  assert(ex < w.nex());
276  assert(ex1 < u1.nex());
277  assert(ex2 < u2.nex());
278  assert(u1.nvol() == w.nvol());
279  assert(u2.nvol() == w.nvol());
280 
281  // int nth = ThreadManager_OpenMP::get_num_threads();
282  // int ith = ThreadManager_OpenMP::get_thread_id();
283  int nth = 1;
284  int ith = 0;
285 
286  int is = w.nvol() * ith / nth;
287  int ns = w.nvol() * (ith + 1) / nth - is;
288 
289  double *g = w.ptr(0, is, ex);
290  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
291  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
292 
293  const int Nc = w.nc();
294  const int Nc2 = 2 * Nc;
295  const int Ndf = Nc2 * Nc;
296 
297  for (int site = 0; site < ns; ++site) {
298  int ig = Ndf * site;
299 
300  for (int ic2 = 0; ic2 < Nc; ++ic2) {
301  for (int ic1 = 0; ic1 < Nc; ++ic1) {
302  int jg2 = ic2 * 2 + ig;
303  int jg1 = ic1 * Nc2 + ig;
304  g[2 * ic2 + Nc2 * ic1 + ig] +=
305  ff * mult_Gnn_r(&g1[jg1], &g2[jg2], Nc);
306  g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
307  ff * mult_Gnn_i(&g1[jg1], &g2[jg2], Nc);
308  }
309  }
310  }
311 }
312 
313 
314 //====================================================================
315 void multadd_Field_Gdn(Field_G& w, const int ex,
316  const Field_G& u1, int ex1,
317  const Field_G& u2, int ex2,
318  const double ff)
319 {
320  assert(ex < w.nex());
321  assert(ex1 < u1.nex());
322  assert(ex2 < u2.nex());
323  assert(u1.nvol() == w.nvol());
324  assert(u2.nvol() == w.nvol());
325 
326  // int nth = ThreadManager_OpenMP::get_num_threads();
327  // int ith = ThreadManager_OpenMP::get_thread_id();
328  int nth = 1;
329  int ith = 0;
330 
331  int is = w.nvol() * ith / nth;
332  int ns = w.nvol() * (ith + 1) / nth - is;
333 
334  double *g = w.ptr(0, is, ex);
335  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
336  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
337 
338  const int Nc = w.nc();
339  const int Nc2 = 2 * Nc;
340  const int Ndf = Nc2 * Nc;
341 
342  for (int site = 0; site < ns; ++site) {
343  int ig = Ndf * site;
344 
345  for (int ic2 = 0; ic2 < Nc; ++ic2) {
346  for (int ic1 = 0; ic1 < Nc; ++ic1) {
347  int jg2 = ic2 * 2 + ig;
348  int jg1 = ic1 * 2 + ig;
349  g[2 * ic2 + Nc2 * ic1 + ig] +=
350  ff * mult_Gdn_r(&g1[jg1], &g2[jg2], Nc);
351  g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
352  ff * mult_Gdn_i(&g1[jg1], &g2[jg2], Nc);
353  }
354  }
355  }
356 }
357 
358 
359 //====================================================================
360 void multadd_Field_Gnd(Field_G& w, const int ex,
361  const Field_G& u1, int ex1,
362  const Field_G& u2, int ex2,
363  const double ff)
364 {
365  assert(ex < w.nex());
366  assert(ex1 < u1.nex());
367  assert(ex2 < u2.nex());
368  assert(u1.nvol() == w.nvol());
369  assert(u2.nvol() == w.nvol());
370 
371  // int nth = ThreadManager_OpenMP::get_num_threads();
372  // int ith = ThreadManager_OpenMP::get_thread_id();
373  int nth = 1;
374  int ith = 0;
375 
376  int is = w.nvol() * ith / nth;
377  int ns = w.nvol() * (ith + 1) / nth - is;
378 
379  double *g = w.ptr(0, is, ex);
380  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
381  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
382 
383  const int Nc = w.nc();
384  const int Nc2 = 2 * Nc;
385  const int Ndf = Nc2 * Nc;
386 
387  for (int site = 0; site < ns; ++site) {
388  int ig = Ndf * site;
389 
390  for (int ic2 = 0; ic2 < Nc; ++ic2) {
391  for (int ic1 = 0; ic1 < Nc; ++ic1) {
392  int jg2 = ic2 * Nc2 + ig;
393  int jg1 = ic1 * Nc2 + ig;
394  g[2 * ic2 + Nc2 * ic1 + ig] +=
395  ff * mult_Gnd_r(&g1[jg1], &g2[jg2], Nc);
396  g[2 * ic2 + 1 + Nc2 * ic1 + ig]
397  += ff * mult_Gnd_i(&g1[jg1], &g2[jg2], Nc);
398  }
399  }
400  }
401 }
402 
403 
404 //====================================================================
405 void multadd_Field_Gdd(Field_G& w, const int ex,
406  const Field_G& u1, int ex1,
407  const Field_G& u2, int ex2,
408  const double ff)
409 {
410  assert(ex < w.nex());
411  assert(ex1 < u1.nex());
412  assert(ex2 < u2.nex());
413  assert(u1.nvol() == w.nvol());
414  assert(u2.nvol() == w.nvol());
415 
416  // int nth = ThreadManager_OpenMP::get_num_threads();
417  // int ith = ThreadManager_OpenMP::get_thread_id();
418  int nth = 1;
419  int ith = 0;
420 
421  int is = w.nvol() * ith / nth;
422  int ns = w.nvol() * (ith + 1) / nth - is;
423 
424  double *g = w.ptr(0, is, ex);
425  double *g1 = const_cast<Field_G *>(&u1)->ptr(0, is, ex1);
426  double *g2 = const_cast<Field_G *>(&u2)->ptr(0, is, ex2);
427 
428  const int Nc = w.nc();
429  const int Nc2 = 2 * Nc;
430  const int Ndf = Nc2 * Nc;
431 
432  for (int site = 0; site < ns; ++site) {
433  int ig = Ndf * site;
434 
435  for (int ic2 = 0; ic2 < Nc; ++ic2) {
436  for (int ic1 = 0; ic1 < Nc; ++ic1) {
437  int jg2 = ic2 * Nc2 + ig;
438  int jg1 = ic1 * 2 + ig;
439  g[2 * ic2 + Nc2 * ic1 + ig] +=
440  ff * mult_Gdd_r(&g1[jg1], &g2[jg2], Nc);
441  g[2 * ic2 + 1 + Nc2 * ic1 + ig] +=
442  ff * mult_Gdd_i(&g1[jg1], &g2[jg2], Nc);
443  }
444  }
445  }
446 }
447 
448 
449 //====================================================================
450 void at_Field_G(Field_G& w, const int ex)
451 {
452  assert(ex < w.nex());
453 
454  // int nth = ThreadManager_OpenMP::get_num_threads();
455  // int ith = ThreadManager_OpenMP::get_thread_id();
456  int nth = 1;
457  int ith = 0;
458 
459  int is = w.nvol() * ith / nth;
460  int ns = w.nvol() * (ith + 1) / nth - is;
461 
462  double *g = w.ptr(0, is, ex);
463 
464  const int Nc = w.nc();
465  const int Ndf = 2 * Nc * Nc;
466 
467  for (int site = 0; site < ns; ++site) {
468  int ig = Ndf * site;
469 
470  for (int a = 0; a < Nc; ++a) {
471  for (int b = a + 1; b < Nc; ++b) {
472  double re = g[2 * (Nc * a + b) + ig] - g[2 * (Nc * b + a) + ig];
473  double im = g[2 * (Nc * a + b) + 1 + ig] + g[2 * (Nc * b + a) + 1 + ig];
474 
475  g[2 * (Nc * a + b) + ig] = 0.5 * re;
476  g[2 * (Nc * a + b) + 1 + ig] = 0.5 * im;
477 
478  g[2 * (Nc * b + a) + ig] = -0.5 * re;
479  g[2 * (Nc * b + a) + 1 + ig] = 0.5 * im;
480  }
481  }
482  double tr = 0.0;
483  for (int cc = 0; cc < Nc; ++cc) {
484  tr += g[2 * (Nc * cc + cc) + 1 + ig];
485  }
486  tr = tr / Nc;
487  for (int cc = 0; cc < Nc; ++cc) {
488  g[2 * (Nc * cc + cc) + ig] = 0.0;
489  g[2 * (Nc * cc + cc) + 1 + ig] -= tr;
490  }
491  }
492 }
493 
494 
495 //====================================================================
496 void ah_Field_G(Field_G& w, const int ex)
497 {
498  assert(ex < w.nex());
499 
500  // int nth = ThreadManager_OpenMP::get_num_threads();
501  // int ith = ThreadManager_OpenMP::get_thread_id();
502  int nth = 1;
503  int ith = 0;
504 
505  int is = w.nvol() * ith / nth;
506  int ns = w.nvol() * (ith + 1) / nth - is;
507 
508  double *g = w.ptr(0, is, ex);
509 
510  const int Nc = w.nc();
511  const int Ndf = 2 * Nc * Nc;
512 
513  for (int site = 0; site < ns; ++site) {
514  int ig = Ndf * site;
515 
516  for (int a = 0; a < Nc; ++a) {
517  for (int b = a; b < Nc; ++b) {
518  double re = g[2 * (Nc * a + b) + ig] - g[2 * (Nc * b + a) + ig];
519  double im = g[2 * (Nc * a + b) + 1 + ig] + g[2 * (Nc * b + a) + 1 + ig];
520 
521  g[2 * (Nc * a + b) + ig] = 0.5 * re;
522  g[2 * (Nc * a + b) + 1 + ig] = 0.5 * im;
523 
524  g[2 * (Nc * b + a) + ig] = -0.5 * re;
525  g[2 * (Nc * b + a) + 1 + ig] = 0.5 * im;
526  }
527  }
528  }
529 }
530 
531 
532 //====================================================================
533 //============================================================END=====
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:207
Field::m_Nex
int m_Nex
Definition: field.h:45
mult_Field_Gdn
void mult_Field_Gdn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
Definition: field_G_imp.cpp:144
Field::m_Nvol
int m_Nvol
Definition: field.h:44
Field::ptr
double * ptr(const int jin, const int site, const int jex)
Definition: field.h:118
Field::nvol
int nvol() const
Definition: field.h:101
Field_G::check
void check()
check of assumptions for performance implementation.
Definition: field_G_imp.cpp:30
RandomNumbers::gauss_lex_global
virtual void gauss_lex_global(Field &)
gaussian random number defined on global lattice.
Definition: randomNumbers.cpp:52
Field_G::set_random
void set_random(RandomNumbers *rand)
Definition: field_G_imp.cpp:62
mult_Field_Gdd
void mult_Field_Gdd(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
Definition: field_G_imp.cpp:228
multadd_Field_Gdd
void multadd_Field_Gdd(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2, const double ff)
Definition: field_G_imp.cpp:405
SU_N::Mat_SU_N::reunit
Mat_SU_N & reunit()
Definition: mat_SU_N.h:71
Field_G::set_unit
void set_unit()
Definition: field_G_imp.cpp:39
ah_Field_G
void ah_Field_G(Field_G &w, const int ex)
Definition: field_G_imp.cpp:496
Field_G::nc
int nc() const
Definition: field_G.h:80
SU_N::Mat_SU_N
Definition: mat_SU_N.h:35
Field_G
SU(N) gauge field.
Definition: field_G.h:36
multadd_Field_Gdn
void multadd_Field_Gdn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2, const double ff)
Definition: field_G_imp.cpp:315
mult_Field_Gnd
void mult_Field_Gnd(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
Definition: field_G_imp.cpp:186
mult_Field_Gnn
void mult_Field_Gnn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2)
Definition: field_G_imp.cpp:102
field_G.h
Field::nex
int nex() const
Definition: field.h:102
Field_G::reunit
void reunit()
Definition: field_G_imp.cpp:79
multadd_Field_Gnd
void multadd_Field_Gnd(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2, const double ff)
Definition: field_G_imp.cpp:360
at_Field_G
void at_Field_G(Field_G &w, const int ex)
Definition: field_G_imp.cpp:450
RandomNumbers
Base class of random number generators.
Definition: randomNumbers.h:40
SU_N::Mat_SU_N::unit
Mat_SU_N & unit()
Definition: mat_SU_N.h:373
multadd_Field_Gnn
void multadd_Field_Gnn(Field_G &w, const int ex, const Field_G &u1, const int ex1, const Field_G &u2, const int ex2, const double ff)
Definition: field_G_imp.cpp:270
NC
#define NC
Definition: tensorProd.cpp:14
Field_G::set_mat
void set_mat(const int site, const int mn, const Mat_SU_N &U)
Definition: field_G.h:156
Field_G::mat
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:110
Field_G::m_Nc
int m_Nc
Definition: field_G.h:39