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