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