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