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