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