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