Bridge++  Version 1.4.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_CloverTerm.cpp
Go to the documentation of this file.
1 
14 #include "Fopr/fopr_CloverTerm.h"
15 
17 
18 
19 using std::string;
20 
21 namespace Imp_BGQ {
22 #if defined USE_GROUP_SU3
23 #include "fopr_Wilson_impl_SU3.inc"
24 #elif defined USE_GROUP_SU2
25 #include "fopr_Wilson_impl_SU2.inc"
26 #elif defined USE_GROUP_SU_N
27 #include "fopr_Wilson_impl_SU_N.inc"
28 #endif
29 
30 //====================================================================
31 
32  const std::string Fopr_CloverTerm::class_name = "Imp_BGQ::Fopr_CloverTerm";
33 
34 //====================================================================
36  {
37  const string str_vlevel = params.get_string("verbose_level");
38 
39  m_vl = vout.set_verbose_level(str_vlevel);
40 
41  //- fetch and check input parameters
42  double kappa, cSW;
43  std::vector<int> bc;
44 
45  int err = 0;
46  err += params.fetch_double("hopping_parameter", kappa);
47  err += params.fetch_double("clover_coefficient", cSW);
48  err += params.fetch_int_vector("boundary_condition", bc);
49 
50  if (err) {
51  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
52  exit(EXIT_FAILURE);
53  }
54 
55 
56  set_parameters(kappa, cSW, bc);
57  }
58 
59 
60 //====================================================================
61  void Fopr_CloverTerm::set_parameters(double kappa, double cSW,
62  std::vector<int> bc)
63  {
64  //- print input parameters
65  vout.general(m_vl, "%s:\n", class_name.c_str());
66  vout.general(m_vl, " kappa = %12.8f\n", kappa);
67  vout.general(m_vl, " cSW = %12.8f\n", cSW);
68  for (int mu = 0; mu < m_Ndim; ++mu) {
69  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
70  }
71 
72  //- range check
73  // NB. kappa,cSW == 0 is allowed.
74  assert(bc.size() == m_Ndim);
75 
76  //- store values
77  m_kappa = kappa;
78  m_cSW = cSW;
79 
80  // m_boundary.resize(m_Ndim); // already resized in init.
81  for (int mu = 0; mu < m_Ndim; ++mu) {
82  m_boundary[mu] = bc[mu];
83  }
84  }
85 
86 
87 //====================================================================
89  {
90  m_U = (Field_G *)U;
91  set_csw();
92  }
93 
94 
95 //====================================================================
96  void Fopr_CloverTerm::init(string repr)
97  {
102  m_NinF = 2 * m_Nc * m_Nd;
103 
104  m_U = 0;
105 
106  m_repr = repr;
107 
108  m_boundary.resize(m_Ndim);
109  m_SG.resize(m_Ndim * m_Ndim);
110 
111  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);
112 
113  m_GM5 = gmset->get_GM(gmset->GAMMA5);
114 
115  m_SG[sg_index(0, 1)] = gmset->get_GM(gmset->SIGMA12);
116  m_SG[sg_index(1, 2)] = gmset->get_GM(gmset->SIGMA23);
117  m_SG[sg_index(2, 0)] = gmset->get_GM(gmset->SIGMA31);
118  m_SG[sg_index(3, 0)] = gmset->get_GM(gmset->SIGMA41);
119  m_SG[sg_index(3, 1)] = gmset->get_GM(gmset->SIGMA42);
120  m_SG[sg_index(3, 2)] = gmset->get_GM(gmset->SIGMA43);
121 
122  m_SG[sg_index(1, 0)] = m_SG[sg_index(0, 1)].mult(-1);
123  m_SG[sg_index(2, 1)] = m_SG[sg_index(1, 2)].mult(-1);
124  m_SG[sg_index(0, 2)] = m_SG[sg_index(2, 0)].mult(-1);
125  m_SG[sg_index(0, 3)] = m_SG[sg_index(3, 0)].mult(-1);
126  m_SG[sg_index(1, 3)] = m_SG[sg_index(3, 1)].mult(-1);
127  m_SG[sg_index(2, 3)] = m_SG[sg_index(3, 2)].mult(-1);
128 
129  m_SG[sg_index(0, 0)] = gmset->get_GM(gmset->UNITY);
130  m_SG[sg_index(1, 1)] = gmset->get_GM(gmset->UNITY);
131  m_SG[sg_index(2, 2)] = gmset->get_GM(gmset->UNITY);
132  m_SG[sg_index(3, 3)] = gmset->get_GM(gmset->UNITY);
133  // these 4 gamma matrices are actually not used.
134 
135  // m_fopr_w = new Fopr_Wilson(repr);
136  if (m_repr == "Dirac") {
139  } else if (m_repr == "Chiral") {
142  }
143 
144  delete gmset;
145 
146  // m_w2.reset(m_NinF,m_Nvol,1);
147  }
148 
149 
150 //====================================================================
152  {
153  // nothing to do.
154  }
155 
156 
157 //====================================================================
159  {
160  (this->*m_gm5)(v, f);
161  }
162 
163 
164 //====================================================================
166  {
167  int Nvc = 2 * CommonParameters::Nc();
168  int Nd = CommonParameters::Nd();
169 
170  const double *v1 = f.ptr(0);
171  double *v2 = w.ptr(0);
172 
173  int id1 = 0;
174  int id2 = Nvc;
175  int id3 = Nvc * 2;
176  int id4 = Nvc * 3;
177 
178  // threadding applied.
181 
182  int is = m_Nvol * ith / nth;
183  int ns = m_Nvol * (ith + 1) / nth - is;
184 
185  for (int site = is; site < is + ns; ++site) {
186  // for (int site = 0; site < m_Nvol; ++site) {
187  for (int icc = 0; icc < Nvc; icc++) {
188  int in = Nvc * Nd * site;
189  v2[icc + id1 + in] = v1[icc + id3 + in];
190  v2[icc + id2 + in] = v1[icc + id4 + in];
191  v2[icc + id3 + in] = v1[icc + id1 + in];
192  v2[icc + id4 + in] = v1[icc + id2 + in];
193  }
194  }
195  }
196 
197 
198 //====================================================================
200  {
201  int Nvc = 2 * CommonParameters::Nc();
202  int Nd = CommonParameters::Nd();
203 
204  const double *v1 = f.ptr(0);
205  double *v2 = w.ptr(0);
206 
207  int id1 = 0;
208  int id2 = Nvc;
209  int id3 = Nvc * 2;
210  int id4 = Nvc * 3;
211 
212  // threadding applied.
215 
216  int is = m_Nvol * ith / nth;
217  int ns = m_Nvol * (ith + 1) / nth - is;
218 
219  for (int site = is; site < is + ns; ++site) {
220  // for (int site = 0; site < m_Nvol; ++site) {
221  for (int icc = 0; icc < Nvc; icc++) {
222  int in = Nvc * Nd * site;
223  v2[icc + id1 + in] = v1[icc + id1 + in];
224  v2[icc + id2 + in] = v1[icc + id2 + in];
225  v2[icc + id3 + in] = -v1[icc + id3 + in];
226  v2[icc + id4 + in] = -v1[icc + id4 + in];
227  }
228  }
229  }
230 
231 
232 //====================================================================
234  const int mu, const int nu)
235  {
236  assert(mu != nu);
237  mult_iGM(v, m_SG[sg_index(mu, nu)], w);
238  }
239 
240 
241 //====================================================================
243  {
244  mult_csw(v, f);
245  }
246 
247 
248 //====================================================================
250  {
251  (this->*m_csw)(v, w);
252  }
253 
254 
255 //====================================================================
257  {
258  assert(w.nex() == 1);
259 
260  int Nc = CommonParameters::Nc();
261  int Nvc = 2 * Nc;
262  int Ndf = 2 * Nc * Nc;
263  int Nd = CommonParameters::Nd();
264  int Nvol = w.nvol();
265 
266  int id1 = 0;
267  int id2 = Nvc;
268  int id3 = Nvc * 2;
269  int id4 = Nvc * 3;
270 
271  double kappa_cSW = m_kappa * m_cSW;
272 
273  const double *w2 = w.ptr(0);
274  double *v2 = v.ptr(0);
275 
276  double *Bx = m_Bx.ptr(0);
277  double *By = m_By.ptr(0);
278  double *Bz = m_Bz.ptr(0);
279  double *Ex = m_Ex.ptr(0);
280  double *Ey = m_Ey.ptr(0);
281  double *Ez = m_Ez.ptr(0);
282 
283  // threadding applied.
286 
287  int is = m_Nvol * ith / nth;
288  int ns = m_Nvol * (ith + 1) / nth - is;
289 
290  for (int site = is; site < is + ns; ++site) {
291  int iv = Nvc * Nd * site;
292  int ig = Ndf * site;
293 
294  for (int ic = 0; ic < Nc; ++ic) {
295  int icr = 2 * ic;
296  int ici = icr + 1;
297  int icg = ic * Nvc + ig;
298 
299  v2[icr + id1 + iv] = 0.0;
300  v2[ici + id1 + iv] = 0.0;
301  v2[icr + id2 + iv] = 0.0;
302  v2[ici + id2 + iv] = 0.0;
303 
304  v2[icr + id3 + iv] = 0.0;
305  v2[ici + id3 + iv] = 0.0;
306  v2[icr + id4 + iv] = 0.0;
307  v2[ici + id4 + iv] = 0.0;
308 
309  // isigma_23 * Bx
310  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv], Nc);
311  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv], Nc);
312  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv], Nc);
313  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv], Nc);
314 
315  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv], Nc);
316  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv], Nc);
317  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv], Nc);
318  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv], Nc);
319 
320  // isigma_31 * By
321  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv], Nc);
322  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv], Nc);
323  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv], Nc);
324  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv], Nc);
325 
326  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv], Nc);
327  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv], Nc);
328  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv], Nc);
329  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv], Nc);
330 
331  // isigma_12 * Bz
332  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv], Nc);
333  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv], Nc);
334  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv], Nc);
335  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv], Nc);
336 
337  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv], Nc);
338  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv], Nc);
339  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv], Nc);
340  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv], Nc);
341 
342  // isigma_41 * Ex
343  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv], Nc);
344  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv], Nc);
345  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv], Nc);
346  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv], Nc);
347 
348  v2[icr + id3 + iv] -= mult_uv_i(&Ex[icg], &w2[id4 + iv], Nc);
349  v2[ici + id3 + iv] += mult_uv_r(&Ex[icg], &w2[id4 + iv], Nc);
350  v2[icr + id4 + iv] -= mult_uv_i(&Ex[icg], &w2[id3 + iv], Nc);
351  v2[ici + id4 + iv] += mult_uv_r(&Ex[icg], &w2[id3 + iv], Nc);
352 
353  // isigma_42 * Ey
354  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv], Nc);
355  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv], Nc);
356  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv], Nc);
357  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv], Nc);
358 
359  v2[icr + id3 + iv] += mult_uv_r(&Ey[icg], &w2[id4 + iv], Nc);
360  v2[ici + id3 + iv] += mult_uv_i(&Ey[icg], &w2[id4 + iv], Nc);
361  v2[icr + id4 + iv] -= mult_uv_r(&Ey[icg], &w2[id3 + iv], Nc);
362  v2[ici + id4 + iv] -= mult_uv_i(&Ey[icg], &w2[id3 + iv], Nc);
363 
364  // isigma_43 * Ez
365  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv], Nc);
366  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv], Nc);
367  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv], Nc);
368  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv], Nc);
369 
370  v2[icr + id3 + iv] -= mult_uv_i(&Ez[icg], &w2[id3 + iv], Nc);
371  v2[ici + id3 + iv] += mult_uv_r(&Ez[icg], &w2[id3 + iv], Nc);
372  v2[icr + id4 + iv] += mult_uv_i(&Ez[icg], &w2[id4 + iv], Nc);
373  v2[ici + id4 + iv] -= mult_uv_r(&Ez[icg], &w2[id4 + iv], Nc);
374 
375  // v *= m_kappa * m_cSW;
376  v2[icr + id1 + iv] *= kappa_cSW;
377  v2[ici + id1 + iv] *= kappa_cSW;
378  v2[icr + id2 + iv] *= kappa_cSW;
379  v2[ici + id2 + iv] *= kappa_cSW;
380 
381  v2[icr + id3 + iv] *= kappa_cSW;
382  v2[ici + id3 + iv] *= kappa_cSW;
383  v2[icr + id4 + iv] *= kappa_cSW;
384  v2[ici + id4 + iv] *= kappa_cSW;
385  }
386  }
387 #pragma omp barrier
388  }
389 
390 
391 //====================================================================
393  {
394  assert(w.nex() == 1);
395 
396  int Nc = CommonParameters::Nc();
397  int Nvc = 2 * Nc;
398  int Ndf = 2 * Nc * Nc;
399  int Nd = CommonParameters::Nd();
400  int Nvol = w.nvol();
401 
402  int id1 = 0;
403  int id2 = Nvc;
404  int id3 = Nvc * 2;
405  int id4 = Nvc * 3;
406 
407  double kappa_cSW = m_kappa * m_cSW;
408 
409  const double *w2 = w.ptr(0);
410  double *v2 = v.ptr(0);
411 
412  double *Bx = m_Bx.ptr(0);
413  double *By = m_By.ptr(0);
414  double *Bz = m_Bz.ptr(0);
415  double *Ex = m_Ex.ptr(0);
416  double *Ey = m_Ey.ptr(0);
417  double *Ez = m_Ez.ptr(0);
418 
419  // threadding applied.
422 
423  int is = m_Nvol * ith / nth;
424  int ns = m_Nvol * (ith + 1) / nth - is;
425 
426  for (int site = is; site < is + ns; ++site) {
427  int iv = Nvc * Nd * site;
428  int ig = Ndf * site;
429 
430  for (int ic = 0; ic < Nc; ++ic) {
431  int icr = 2 * ic;
432  int ici = icr + 1;
433  int icg = ic * Nvc + ig;
434 
435  v2[icr + id1 + iv] = 0.0;
436  v2[ici + id1 + iv] = 0.0;
437  v2[icr + id2 + iv] = 0.0;
438  v2[ici + id2 + iv] = 0.0;
439 
440  v2[icr + id3 + iv] = 0.0;
441  v2[ici + id3 + iv] = 0.0;
442  v2[icr + id4 + iv] = 0.0;
443  v2[ici + id4 + iv] = 0.0;
444 
445  // isigma_23 * Bx
446  v2[icr + id1 + iv] -= mult_uv_i(&Bx[icg], &w2[id2 + iv], Nc);
447  v2[ici + id1 + iv] += mult_uv_r(&Bx[icg], &w2[id2 + iv], Nc);
448  v2[icr + id2 + iv] -= mult_uv_i(&Bx[icg], &w2[id1 + iv], Nc);
449  v2[ici + id2 + iv] += mult_uv_r(&Bx[icg], &w2[id1 + iv], Nc);
450 
451  v2[icr + id3 + iv] -= mult_uv_i(&Bx[icg], &w2[id4 + iv], Nc);
452  v2[ici + id3 + iv] += mult_uv_r(&Bx[icg], &w2[id4 + iv], Nc);
453  v2[icr + id4 + iv] -= mult_uv_i(&Bx[icg], &w2[id3 + iv], Nc);
454  v2[ici + id4 + iv] += mult_uv_r(&Bx[icg], &w2[id3 + iv], Nc);
455 
456  // isigma_31 * By
457  v2[icr + id1 + iv] += mult_uv_r(&By[icg], &w2[id2 + iv], Nc);
458  v2[ici + id1 + iv] += mult_uv_i(&By[icg], &w2[id2 + iv], Nc);
459  v2[icr + id2 + iv] -= mult_uv_r(&By[icg], &w2[id1 + iv], Nc);
460  v2[ici + id2 + iv] -= mult_uv_i(&By[icg], &w2[id1 + iv], Nc);
461 
462  v2[icr + id3 + iv] += mult_uv_r(&By[icg], &w2[id4 + iv], Nc);
463  v2[ici + id3 + iv] += mult_uv_i(&By[icg], &w2[id4 + iv], Nc);
464  v2[icr + id4 + iv] -= mult_uv_r(&By[icg], &w2[id3 + iv], Nc);
465  v2[ici + id4 + iv] -= mult_uv_i(&By[icg], &w2[id3 + iv], Nc);
466 
467  // isigma_12 * Bz
468  v2[icr + id1 + iv] -= mult_uv_i(&Bz[icg], &w2[id1 + iv], Nc);
469  v2[ici + id1 + iv] += mult_uv_r(&Bz[icg], &w2[id1 + iv], Nc);
470  v2[icr + id2 + iv] += mult_uv_i(&Bz[icg], &w2[id2 + iv], Nc);
471  v2[ici + id2 + iv] -= mult_uv_r(&Bz[icg], &w2[id2 + iv], Nc);
472 
473  v2[icr + id3 + iv] -= mult_uv_i(&Bz[icg], &w2[id3 + iv], Nc);
474  v2[ici + id3 + iv] += mult_uv_r(&Bz[icg], &w2[id3 + iv], Nc);
475  v2[icr + id4 + iv] += mult_uv_i(&Bz[icg], &w2[id4 + iv], Nc);
476  v2[ici + id4 + iv] -= mult_uv_r(&Bz[icg], &w2[id4 + iv], Nc);
477 
478  // isigma_41 * Ex
479  v2[icr + id1 + iv] += mult_uv_i(&Ex[icg], &w2[id4 + iv], Nc);
480  v2[ici + id1 + iv] -= mult_uv_r(&Ex[icg], &w2[id4 + iv], Nc);
481  v2[icr + id2 + iv] += mult_uv_i(&Ex[icg], &w2[id3 + iv], Nc);
482  v2[ici + id2 + iv] -= mult_uv_r(&Ex[icg], &w2[id3 + iv], Nc);
483 
484  v2[icr + id3 + iv] += mult_uv_i(&Ex[icg], &w2[id2 + iv], Nc);
485  v2[ici + id3 + iv] -= mult_uv_r(&Ex[icg], &w2[id2 + iv], Nc);
486  v2[icr + id4 + iv] += mult_uv_i(&Ex[icg], &w2[id1 + iv], Nc);
487  v2[ici + id4 + iv] -= mult_uv_r(&Ex[icg], &w2[id1 + iv], Nc);
488 
489  // isigma_42 * Ey
490  v2[icr + id1 + iv] -= mult_uv_r(&Ey[icg], &w2[id4 + iv], Nc);
491  v2[ici + id1 + iv] -= mult_uv_i(&Ey[icg], &w2[id4 + iv], Nc);
492  v2[icr + id2 + iv] += mult_uv_r(&Ey[icg], &w2[id3 + iv], Nc);
493  v2[ici + id2 + iv] += mult_uv_i(&Ey[icg], &w2[id3 + iv], Nc);
494 
495  v2[icr + id3 + iv] -= mult_uv_r(&Ey[icg], &w2[id2 + iv], Nc);
496  v2[ici + id3 + iv] -= mult_uv_i(&Ey[icg], &w2[id2 + iv], Nc);
497  v2[icr + id4 + iv] += mult_uv_r(&Ey[icg], &w2[id1 + iv], Nc);
498  v2[ici + id4 + iv] += mult_uv_i(&Ey[icg], &w2[id1 + iv], Nc);
499 
500  // isigma_43 * Ez
501  v2[icr + id1 + iv] += mult_uv_i(&Ez[icg], &w2[id3 + iv], Nc);
502  v2[ici + id1 + iv] -= mult_uv_r(&Ez[icg], &w2[id3 + iv], Nc);
503  v2[icr + id2 + iv] -= mult_uv_i(&Ez[icg], &w2[id4 + iv], Nc);
504  v2[ici + id2 + iv] += mult_uv_r(&Ez[icg], &w2[id4 + iv], Nc);
505 
506  v2[icr + id3 + iv] += mult_uv_i(&Ez[icg], &w2[id1 + iv], Nc);
507  v2[ici + id3 + iv] -= mult_uv_r(&Ez[icg], &w2[id1 + iv], Nc);
508  v2[icr + id4 + iv] -= mult_uv_i(&Ez[icg], &w2[id2 + iv], Nc);
509  v2[ici + id4 + iv] += mult_uv_r(&Ez[icg], &w2[id2 + iv], Nc);
510 
511  // v *= m_kappa * m_cSW;
512  v2[icr + id1 + iv] *= kappa_cSW;
513  v2[ici + id1 + iv] *= kappa_cSW;
514  v2[icr + id2 + iv] *= kappa_cSW;
515  v2[ici + id2 + iv] *= kappa_cSW;
516 
517  v2[icr + id3 + iv] *= kappa_cSW;
518  v2[ici + id3 + iv] *= kappa_cSW;
519  v2[icr + id4 + iv] *= kappa_cSW;
520  v2[ici + id4 + iv] *= kappa_cSW;
521  }
522  }
523 #pragma omp barrier
524  }
525 
526 
527 //====================================================================
529  {
530  set_fieldstrength(m_Bx, 1, 2);
531  set_fieldstrength(m_By, 2, 0);
532  set_fieldstrength(m_Bz, 0, 1);
533  set_fieldstrength(m_Ex, 3, 0);
534  set_fieldstrength(m_Ey, 3, 1);
535  set_fieldstrength(m_Ez, 3, 2);
536  }
537 
538 
539 //====================================================================
541  const int mu, const int nu)
542  {
544 
545  assert(nth == 1);
546  // this function must be called in single thread region.
547 
548  m_staple.upper(m_Cup, *m_U, mu, nu); // these staple constructions
549  m_staple.lower(m_Cdn, *m_U, mu, nu); // are multi-threaded.
550 
551  //#pragma omp parallel
552  {
553  mult_Field_Gnd(Fst, 0, *m_U, mu, m_Cup, 0);
554  multadd_Field_Gnd(Fst, 0, *m_U, mu, m_Cdn, 0, -1.0);
555 
556  mult_Field_Gdn(m_v1, 0, m_Cup, 0, *m_U, mu);
557  multadd_Field_Gdn(m_v1, 0, m_Cdn, 0, *m_U, mu, -1.0);
558  }
559 
560  m_shift.forward(m_v2, m_v1, mu);
561 
562  //#pragma omp parallel
563  {
564  axpy(Fst, 1.0, m_v2);
565  ah_Field_G(Fst, 0);
566  scal(Fst, 0.25);
567  }
568  }
569 
570 
571 //====================================================================
573  {
574  // The following counting explicitly depends on the implementation
575  // and to be recalculated when the code is modified.
576  // Present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]
577 
578  int Lvol = CommonParameters::Lvol();
579 
580  double flop_site
581  = static_cast<double>(m_Nc * m_Nd * (2 + 48 * m_Nc));
582  double flop = flop_site * static_cast<double>(Lvol);
583 
584  return flop;
585  }
586 
587 
588 //====================================================================
589 }
590 //============================================================END=====
void(Fopr_CloverTerm::* m_gm5)(Field &, const Field &)
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:282
static const std::string class_name
BridgeIO vout
Definition: bridgeIO.cpp:495
void ah_Field_G(Field_G &W, const int ex)
static int get_num_threads()
returns available number of threads.
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:142
void mult_csw_chiral(Field &, const Field &)
Field_G m_v2
for calculation of field strength.
void general(const char *format,...)
Definition: bridgeIO.cpp:195
GammaMatrix get_GM(GMspecies spec)
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)
Container of Field-type object.
Definition: field.h:39
int sg_index(int mu, int nu)
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:211
void mult_csw(Field &, const Field &)
int nvol() const
Definition: field.h:116
void mult_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Class for parameters.
Definition: parameters.h:46
std::vector< GammaMatrix > m_SG
static int Lvol()
void lower(Field_G &, const Field_G &, const int mu, const int nu)
constructs lower staple in mu-nu plane.
Definition: staple_lex.cpp:180
const Field_G * m_U
pointer to gauge configuration.
void set_config(Field *U)
setting pointer to the gauge configuration.
static int get_thread_id()
returns thread id.
Wilson-type fermion field.
Definition: field_F.h:37
void init(std::string repr)
void set_fieldstrength(Field_G &, const int, const int)
SU(N) gauge field.
Definition: field_G.h:38
void mult_csw_dirac(Field &, const Field &)
void mult_iGM(Field_F &y, const GammaMatrix &gm, const Field_F &x)
gamma matrix multiplication (i is multiplied)
std::vector< int > m_boundary
Bridge::VerboseLevel m_vl
Definition: fopr.h:128
void set_parameters(const Parameters &params)
double flop_count()
this returns the number of floating point operations.
Set of Gamma Matrices: basis class.
int nex() const
Definition: field.h:117
void upper(Field_G &, const Field_G &, const int mu, const int nu)
constructs upper staple in mu-nu plane.
Definition: staple_lex.cpp:156
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:168
void(Fopr_CloverTerm::* m_csw)(Field &, const Field &)
void crucial(const char *format,...)
Definition: bridgeIO.cpp:178
void gm5_dirac(Field &, const Field &)
void mult_gm5(Field &v, const Field &w)
gamma_5 multiplication. [31 Mar 2017 H.Matsufuru]
void mult_isigma(Field_F &, const Field_F &, const int mu, const int nu)
void gm5_chiral(Field &, const Field &)
string get_string(const string &key) const
Definition: parameters.cpp:116
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:294
Field_G m_Ez
field strength.
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131
void mult_sigmaF(Field &, const Field &)
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)
void forward(Field &, const Field &, const int mu)
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)