Bridge++  Ver. 2.0.2
fopr_CloverTerm_impl.cpp
Go to the documentation of this file.
1 
14 #include "fopr_CloverTerm_impl.h"
15 
16 #include "Fopr/fopr_thread-inc.h"
17 
18 namespace Imp {
19 #if defined USE_GROUP_SU3
20 #include "fopr_Wilson_impl_SU3-inc.h"
21 #elif defined USE_GROUP_SU2
22 #include "fopr_Wilson_impl_SU2-inc.h"
23 #elif defined USE_GROUP_SU_N
24 #include "fopr_Wilson_impl_SU_N-inc.h"
25 #endif
26 
27  const std::string Fopr_CloverTerm::class_name = "Imp::Fopr_CloverTerm";
28 
29 //====================================================================
30  void Fopr_CloverTerm::init(const Parameters& params)
31  {
32  ThreadManager::assert_single_thread(class_name);
33 
34  m_vl = CommonParameters::Vlevel();
35 
36  vout.general(m_vl, "%s: construction\n", class_name.c_str());
37  vout.increase_indent();
38 
39  m_Nvol = CommonParameters::Nvol();
40  m_Ndim = CommonParameters::Ndim();
41  m_Nc = CommonParameters::Nc();
42  m_Nd = CommonParameters::Nd();
43  m_NinF = 2 * m_Nc * m_Nd;
44 
45  m_boundary.resize(m_Ndim);
46  m_SG.resize(m_Ndim * m_Ndim);
47 
48  int Ndf = 2 * m_Nc * m_Nc;
49  m_shift = new ShiftField_lex(Ndf);
50 
51  std::string repr;
52  if (!params.fetch_string("gamma_matrix_type", repr)) {
53  m_repr = repr;
54  } else {
55  m_repr = "Dirac"; // default
56  vout.general(m_vl, "gamma_matrix_type is not given: defalt = %s\n",
57  m_repr.c_str());
58  }
59  if ((m_repr != "Dirac") && (m_repr != "Chiral")) {
60  vout.crucial("Error in %s: irrelevant mult mode = %s\n",
61  class_name.c_str(), m_repr.c_str());
62  exit(EXIT_FAILURE);
63  }
64 
65  set_parameters(params);
66 
67  setup_gamma_matrices();
68 
69  m_U = 0;
70 
71  vout.decrease_indent();
72  vout.general(m_vl, "%s: construction finished.\n",
73  class_name.c_str());
74  }
75 
76 
77 //====================================================================
78  void Fopr_CloverTerm::init(const std::string repr)
79  {
80  ThreadManager::assert_single_thread(class_name);
81 
82  m_vl = CommonParameters::Vlevel();
83 
84  vout.general(m_vl, "%s: construction (obsolete)\n",
85  class_name.c_str());
86 
87  m_Nvol = CommonParameters::Nvol();
88  m_Ndim = CommonParameters::Ndim();
89  m_Nc = CommonParameters::Nc();
90  m_Nd = CommonParameters::Nd();
91  m_NinF = 2 * m_Nc * m_Nd;
92 
93  m_boundary.resize(m_Ndim);
94  m_SG.resize(m_Ndim * m_Ndim);
95 
96  m_repr = repr;
97 
98  setup_gamma_matrices();
99 
100  int Ndf = 2 * m_Nc * m_Nc;
101  m_shift = new ShiftField_lex(Ndf);
102 
103  m_U = 0;
104 
105  vout.general(m_vl, "%s: construction finished.\n",
106  class_name.c_str());
107  }
108 
109 
110 //====================================================================
111  void Fopr_CloverTerm::tidyup()
112  {
113  delete m_shift;
114  }
115 
116 
117 //====================================================================
118  void Fopr_CloverTerm::set_parameters(const Parameters& params)
119  {
120 #pragma omp barrier
121 
122  int ith = ThreadManager::get_thread_id();
123 
124  std::string vlevel;
125  if (!params.fetch_string("verbose_level", vlevel)) {
126  if (ith == 0) m_vl = vout.set_verbose_level(vlevel);
127  }
128 #pragma omp barrier
129 
130  //- fetch and check input parameters
131  double kappa, cSW;
132  std::vector<int> bc;
133 
134  int err = 0;
135  err += params.fetch_double("hopping_parameter", kappa);
136  err += params.fetch_double("clover_coefficient", cSW);
137  err += params.fetch_int_vector("boundary_condition", bc);
138 
139  if (err) {
140  vout.crucial("Error at %s: input parameter not found.\n",
141  class_name.c_str());
142  exit(EXIT_FAILURE);
143  }
144 
145  set_parameters(kappa, cSW, bc);
146  }
147 
148 
149 //====================================================================
150  void Fopr_CloverTerm::set_parameters(const double kappa,
151  const double cSW,
152  const std::vector<int> bc)
153  {
154  assert(bc.size() == m_Ndim);
155 
156 #pragma omp barrier
157  int ith = ThreadManager::get_thread_id();
158  if (ith == 0) {
159  m_kappa = kappa;
160  m_cSW = cSW;
161  m_boundary = bc;
162  }
163 #pragma omp barrier
164 
165  vout.general(m_vl, "%s: input parameters\n", class_name.c_str());
166  vout.general(m_vl, " gamma matrix type = %s\n", m_repr.c_str());
167  vout.general(m_vl, " kappa = %12.8f\n", kappa);
168  vout.general(m_vl, " cSW = %12.8f\n", cSW);
169  for (int mu = 0; mu < m_Ndim; ++mu) {
170  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
171  }
172  }
173 
174 
175 //====================================================================
176  void Fopr_CloverTerm::get_parameters(Parameters& params) const
177  {
178  params.set_double("hopping_parameter", m_kappa);
179  params.set_double("clover_coefficient", m_cSW);
180  params.set_int_vector("boundary_condition", m_boundary);
181  params.set_string("gamma_matrix_type", m_repr);
182 
183  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
184  }
185 
186 
187 //====================================================================
188  void Fopr_CloverTerm::setup_gamma_matrices()
189  {
190  GammaMatrixSet *gmset = GammaMatrixSet::New(m_repr);
191 
192  m_GM5 = gmset->get_GM(gmset->GAMMA5);
193 
194  m_SG[sg_index(0, 1)] = gmset->get_GM(gmset->SIGMA12);
195  m_SG[sg_index(1, 2)] = gmset->get_GM(gmset->SIGMA23);
196  m_SG[sg_index(2, 0)] = gmset->get_GM(gmset->SIGMA31);
197  m_SG[sg_index(3, 0)] = gmset->get_GM(gmset->SIGMA41);
198  m_SG[sg_index(3, 1)] = gmset->get_GM(gmset->SIGMA42);
199  m_SG[sg_index(3, 2)] = gmset->get_GM(gmset->SIGMA43);
200 
201  m_SG[sg_index(1, 0)] = m_SG[sg_index(0, 1)].mult(-1);
202  m_SG[sg_index(2, 1)] = m_SG[sg_index(1, 2)].mult(-1);
203  m_SG[sg_index(0, 2)] = m_SG[sg_index(2, 0)].mult(-1);
204  m_SG[sg_index(0, 3)] = m_SG[sg_index(3, 0)].mult(-1);
205  m_SG[sg_index(1, 3)] = m_SG[sg_index(3, 1)].mult(-1);
206  m_SG[sg_index(2, 3)] = m_SG[sg_index(3, 2)].mult(-1);
207 
208  m_SG[sg_index(0, 0)] = gmset->get_GM(gmset->UNITY);
209  m_SG[sg_index(1, 1)] = gmset->get_GM(gmset->UNITY);
210  m_SG[sg_index(2, 2)] = gmset->get_GM(gmset->UNITY);
211  m_SG[sg_index(3, 3)] = gmset->get_GM(gmset->UNITY);
212  // these 4 gamma matrices are actually not used.
213 
214  delete gmset;
215  }
216 
217 
218 //====================================================================
219  void Fopr_CloverTerm::set_config(Field *U)
220  {
221  int nth = ThreadManager::get_num_threads();
222  vout.detailed(m_vl, "%s: set_config is called: num_threads = %d\n",
223  class_name.c_str(), nth);
224 
225  if (nth > 1) {
226  set_config_impl(U);
227  } else {
228  set_config_omp(U);
229  }
230 
231  vout.detailed(m_vl, "%s: set_config finished.\n", class_name.c_str());
232  }
233 
234 
235 //====================================================================
236  void Fopr_CloverTerm::set_config_omp(Field *U)
237  {
238 #pragma omp parallel
239  {
240  set_config_impl(U);
241  }
242  }
243 
244 
245 //====================================================================
246  void Fopr_CloverTerm::set_config_impl(Field *U)
247  {
248 #pragma omp barrier
249 
250  int ith = ThreadManager::get_thread_id();
251  if (ith == 0) m_U = (Field_G *)U;
252 
253 #pragma omp barrier
254 
255  set_csw();
256 #pragma omp barrier
257  }
258 
259 
260 //====================================================================
261  void Fopr_CloverTerm::set_mode(const std::string mode)
262  {
263 #pragma omp barrier
264  int ith = ThreadManager::get_thread_id();
265  if (ith == 0) m_mode = mode;
266 #pragma omp barrier
267  }
268 
269 
270 //====================================================================
271  void Fopr_CloverTerm::mult(Field& v, const Field& w)
272  {
273  // csw kappa sigma_{mu nu} F_{mu nu}
274  if ((m_mode == "D") || (m_mode == "F")) {
275  mult_sigmaF(v, w);
276  } else {
277  vout.crucial("Error at %s: undefined mode = %s\n",
278  class_name.c_str(), m_mode.c_str());
279  exit(EXIT_FAILURE);
280  }
281  }
282 
283 
284 //====================================================================
285  void Fopr_CloverTerm::mult_dag(Field& v, const Field& w)
286  {
287  mult(v, w);
288  }
289 
290 
291 //====================================================================
292  void Fopr_CloverTerm::mult_gm5(Field& v, const Field& w)
293  {
294  if (m_repr == "Dirac") {
295  gm5_dirac(v, w);
296  } else if (m_repr == "Chiral") {
297  gm5_chiral(v, w);
298  }
299  }
300 
301 
302 //====================================================================
303  void Fopr_CloverTerm::gm5_dirac(Field& w, const Field& f)
304  {
305  const int Nvc = 2 * CommonParameters::Nc();
306  const int Nd = CommonParameters::Nd();
307 
308  const double *v1 = f.ptr(0);
309  double *v2 = w.ptr(0);
310 
311  const int id1 = 0;
312  const int id2 = Nvc;
313  const int id3 = Nvc * 2;
314  const int id4 = Nvc * 3;
315 
316  // threadding applied.
317  const int Nthread = ThreadManager::get_num_threads();
318  const int i_thread = ThreadManager::get_thread_id();
319 
320  const int is = m_Nvol * i_thread / Nthread;
321  const int ns = m_Nvol * (i_thread + 1) / Nthread - is;
322 
323  for (int site = is; site < is + ns; ++site) {
324  // for (int site = 0; site < m_Nvol; ++site) {
325  for (int icc = 0; icc < Nvc; icc++) {
326  int in = Nvc * Nd * site;
327 
328  v2[icc + id1 + in] = v1[icc + id3 + in];
329  v2[icc + id2 + in] = v1[icc + id4 + in];
330  v2[icc + id3 + in] = v1[icc + id1 + in];
331  v2[icc + id4 + in] = v1[icc + id2 + in];
332  }
333  }
334  }
335 
336 
337 //====================================================================
338  void Fopr_CloverTerm::gm5_chiral(Field& w, const Field& f)
339  {
340  const int Nvc = 2 * CommonParameters::Nc();
341  const int Nd = CommonParameters::Nd();
342 
343  const double *v1 = f.ptr(0);
344  double *v2 = w.ptr(0);
345 
346  const int id1 = 0;
347  const int id2 = Nvc;
348  const int id3 = Nvc * 2;
349  const int id4 = Nvc * 3;
350 
351  // threadding applied.
352  const int Nthread = ThreadManager::get_num_threads();
353  const int i_thread = ThreadManager::get_thread_id();
354 
355  const int is = m_Nvol * i_thread / Nthread;
356  const int ns = m_Nvol * (i_thread + 1) / Nthread - is;
357 
358  for (int site = is; site < is + ns; ++site) {
359  // for (int site = 0; site < m_Nvol; ++site) {
360  for (int icc = 0; icc < Nvc; icc++) {
361  int in = Nvc * Nd * site;
362 
363  v2[icc + id1 + in] = v1[icc + id1 + in];
364  v2[icc + id2 + in] = v1[icc + id2 + in];
365  v2[icc + id3 + in] = -v1[icc + id3 + in];
366  v2[icc + id4 + in] = -v1[icc + id4 + in];
367  }
368  }
369  }
370 
371 
372 //====================================================================
373  void Fopr_CloverTerm::mult_isigma(Field_F& v, const Field_F& w,
374  const int mu, const int nu)
375  {
376  assert(mu != nu);
377  mult_iGM(v, m_SG[sg_index(mu, nu)], w);
378  }
379 
380 
381 //====================================================================
382  void Fopr_CloverTerm::mult_sigmaF(Field& v, const Field& f)
383  {
384  // multiplies csw kappa sigma_{mu nu} F_{mu nu}
385  // NOTE: this is NOT 1 - csw kappa sigma_{mu nu} F_{mu nu}
386 
387  mult_csw(v, f);
388  }
389 
390 
391 //====================================================================
392  void Fopr_CloverTerm::mult_csw(Field& v, const Field& w)
393  {
394  // multiplies csw kappa sigma_{mu nu} F_{mu nu}
395  // NOTE: this is NOT 1 - csw kappa sigma_{mu nu} F_{mu nu}
396 
397  if (m_repr == "Dirac") {
398  mult_csw_dirac(v, w);
399  } else if (m_repr == "Chiral") {
400  mult_csw_chiral(v, w);
401  }
402  }
403 
404 
405 //====================================================================
406  void Fopr_CloverTerm::mult_csw_dirac(Field& v, const Field& w)
407  {
408  // multiplies csw kappa sigma_{mu nu} F_{mu nu}
409  // NOTE: this is NOT 1 - csw kappa sigma_{mu nu} F_{mu nu}
410  assert(w.nex() == 1);
411 
412  const int Nc = CommonParameters::Nc();
413  const int Nvc = 2 * Nc;
414  const int Ndf = 2 * Nc * Nc;
415  const int Nd = CommonParameters::Nd();
416  const int Nvol = w.nvol();
417 
418  const int id1 = 0;
419  const int id2 = Nvc;
420  const int id3 = Nvc * 2;
421  const int id4 = Nvc * 3;
422 
423  const double kappa_cSW = m_kappa * m_cSW;
424 
425  const double *w2 = w.ptr(0);
426  double *v2 = v.ptr(0);
427 
428  double *Bx = m_Bx.ptr(0);
429  double *By = m_By.ptr(0);
430  double *Bz = m_Bz.ptr(0);
431  double *Ex = m_Ex.ptr(0);
432  double *Ey = m_Ey.ptr(0);
433  double *Ez = m_Ez.ptr(0);
434 
435  // threading applied.
436  const int Nthread = ThreadManager::get_num_threads();
437  const int i_thread = ThreadManager::get_thread_id();
438 
439  const int is = m_Nvol * i_thread / Nthread;
440  const int ns = m_Nvol * (i_thread + 1) / Nthread - is;
441 
442  for (int site = is; site < is + ns; ++site) {
443  int iv = Nvc * Nd * site;
444  int ig = Ndf * site;
445 
446  for (int ic = 0; ic < Nc; ++ic) {
447  int ic_r = 2 * ic;
448  int ic_i = ic_r + 1;
449  int ic_g = ic * Nvc + ig;
450 
451  v2[ic_r + id1 + iv] = 0.0;
452  v2[ic_i + id1 + iv] = 0.0;
453  v2[ic_r + id2 + iv] = 0.0;
454  v2[ic_i + id2 + iv] = 0.0;
455 
456  v2[ic_r + id3 + iv] = 0.0;
457  v2[ic_i + id3 + iv] = 0.0;
458  v2[ic_r + id4 + iv] = 0.0;
459  v2[ic_i + id4 + iv] = 0.0;
460 
461  // isigma_23 * Bx
462  v2[ic_r + id1 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id2 + iv], Nc);
463  v2[ic_i + id1 + iv] += mult_uv_r(&Bx[ic_g], &w2[id2 + iv], Nc);
464  v2[ic_r + id2 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id1 + iv], Nc);
465  v2[ic_i + id2 + iv] += mult_uv_r(&Bx[ic_g], &w2[id1 + iv], Nc);
466 
467  v2[ic_r + id3 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id4 + iv], Nc);
468  v2[ic_i + id3 + iv] += mult_uv_r(&Bx[ic_g], &w2[id4 + iv], Nc);
469  v2[ic_r + id4 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id3 + iv], Nc);
470  v2[ic_i + id4 + iv] += mult_uv_r(&Bx[ic_g], &w2[id3 + iv], Nc);
471 
472  // isigma_31 * By
473  v2[ic_r + id1 + iv] += mult_uv_r(&By[ic_g], &w2[id2 + iv], Nc);
474  v2[ic_i + id1 + iv] += mult_uv_i(&By[ic_g], &w2[id2 + iv], Nc);
475  v2[ic_r + id2 + iv] -= mult_uv_r(&By[ic_g], &w2[id1 + iv], Nc);
476  v2[ic_i + id2 + iv] -= mult_uv_i(&By[ic_g], &w2[id1 + iv], Nc);
477 
478  v2[ic_r + id3 + iv] += mult_uv_r(&By[ic_g], &w2[id4 + iv], Nc);
479  v2[ic_i + id3 + iv] += mult_uv_i(&By[ic_g], &w2[id4 + iv], Nc);
480  v2[ic_r + id4 + iv] -= mult_uv_r(&By[ic_g], &w2[id3 + iv], Nc);
481  v2[ic_i + id4 + iv] -= mult_uv_i(&By[ic_g], &w2[id3 + iv], Nc);
482 
483  // isigma_12 * Bz
484  v2[ic_r + id1 + iv] -= mult_uv_i(&Bz[ic_g], &w2[id1 + iv], Nc);
485  v2[ic_i + id1 + iv] += mult_uv_r(&Bz[ic_g], &w2[id1 + iv], Nc);
486  v2[ic_r + id2 + iv] += mult_uv_i(&Bz[ic_g], &w2[id2 + iv], Nc);
487  v2[ic_i + id2 + iv] -= mult_uv_r(&Bz[ic_g], &w2[id2 + iv], Nc);
488 
489  v2[ic_r + id3 + iv] -= mult_uv_i(&Bz[ic_g], &w2[id3 + iv], Nc);
490  v2[ic_i + id3 + iv] += mult_uv_r(&Bz[ic_g], &w2[id3 + iv], Nc);
491  v2[ic_r + id4 + iv] += mult_uv_i(&Bz[ic_g], &w2[id4 + iv], Nc);
492  v2[ic_i + id4 + iv] -= mult_uv_r(&Bz[ic_g], &w2[id4 + iv], Nc);
493 
494  // isigma_41 * Ex
495  v2[ic_r + id1 + iv] += mult_uv_i(&Ex[ic_g], &w2[id4 + iv], Nc);
496  v2[ic_i + id1 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id4 + iv], Nc);
497  v2[ic_r + id2 + iv] += mult_uv_i(&Ex[ic_g], &w2[id3 + iv], Nc);
498  v2[ic_i + id2 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id3 + iv], Nc);
499 
500  v2[ic_r + id3 + iv] += mult_uv_i(&Ex[ic_g], &w2[id2 + iv], Nc);
501  v2[ic_i + id3 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id2 + iv], Nc);
502  v2[ic_r + id4 + iv] += mult_uv_i(&Ex[ic_g], &w2[id1 + iv], Nc);
503  v2[ic_i + id4 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id1 + iv], Nc);
504 
505  // isigma_42 * Ey
506  v2[ic_r + id1 + iv] -= mult_uv_r(&Ey[ic_g], &w2[id4 + iv], Nc);
507  v2[ic_i + id1 + iv] -= mult_uv_i(&Ey[ic_g], &w2[id4 + iv], Nc);
508  v2[ic_r + id2 + iv] += mult_uv_r(&Ey[ic_g], &w2[id3 + iv], Nc);
509  v2[ic_i + id2 + iv] += mult_uv_i(&Ey[ic_g], &w2[id3 + iv], Nc);
510 
511  v2[ic_r + id3 + iv] -= mult_uv_r(&Ey[ic_g], &w2[id2 + iv], Nc);
512  v2[ic_i + id3 + iv] -= mult_uv_i(&Ey[ic_g], &w2[id2 + iv], Nc);
513  v2[ic_r + id4 + iv] += mult_uv_r(&Ey[ic_g], &w2[id1 + iv], Nc);
514  v2[ic_i + id4 + iv] += mult_uv_i(&Ey[ic_g], &w2[id1 + iv], Nc);
515 
516  // isigma_43 * Ez
517  v2[ic_r + id1 + iv] += mult_uv_i(&Ez[ic_g], &w2[id3 + iv], Nc);
518  v2[ic_i + id1 + iv] -= mult_uv_r(&Ez[ic_g], &w2[id3 + iv], Nc);
519  v2[ic_r + id2 + iv] -= mult_uv_i(&Ez[ic_g], &w2[id4 + iv], Nc);
520  v2[ic_i + id2 + iv] += mult_uv_r(&Ez[ic_g], &w2[id4 + iv], Nc);
521 
522  v2[ic_r + id3 + iv] += mult_uv_i(&Ez[ic_g], &w2[id1 + iv], Nc);
523  v2[ic_i + id3 + iv] -= mult_uv_r(&Ez[ic_g], &w2[id1 + iv], Nc);
524  v2[ic_r + id4 + iv] -= mult_uv_i(&Ez[ic_g], &w2[id2 + iv], Nc);
525  v2[ic_i + id4 + iv] += mult_uv_r(&Ez[ic_g], &w2[id2 + iv], Nc);
526 
527  // v *= m_kappa * m_cSW;
528  v2[ic_r + id1 + iv] *= kappa_cSW;
529  v2[ic_i + id1 + iv] *= kappa_cSW;
530  v2[ic_r + id2 + iv] *= kappa_cSW;
531  v2[ic_i + id2 + iv] *= kappa_cSW;
532 
533  v2[ic_r + id3 + iv] *= kappa_cSW;
534  v2[ic_i + id3 + iv] *= kappa_cSW;
535  v2[ic_r + id4 + iv] *= kappa_cSW;
536  v2[ic_i + id4 + iv] *= kappa_cSW;
537  }
538  }
539 #pragma omp barrier
540  }
541 
542 
543 //====================================================================
544  void Fopr_CloverTerm::mult_csw_chiral(Field& v, const Field& w)
545  {
546  // multiplies csw kappa sigma_{mu nu} F_{mu nu}
547  // NOTE: this is NOT 1 - csw kappa sigma_{mu nu} F_{mu nu}
548  assert(w.nex() == 1);
549 
550  const int Nc = CommonParameters::Nc();
551  const int Nvc = 2 * Nc;
552  const int Ndf = 2 * Nc * Nc;
553  const int Nd = CommonParameters::Nd();
554  const int Nvol = w.nvol();
555 
556  const int id1 = 0;
557  const int id2 = Nvc;
558  const int id3 = Nvc * 2;
559  const int id4 = Nvc * 3;
560 
561  const double kappa_cSW = m_kappa * m_cSW;
562 
563  const double *w2 = w.ptr(0);
564  double *v2 = v.ptr(0);
565 
566  double *Bx = m_Bx.ptr(0);
567  double *By = m_By.ptr(0);
568  double *Bz = m_Bz.ptr(0);
569  double *Ex = m_Ex.ptr(0);
570  double *Ey = m_Ey.ptr(0);
571  double *Ez = m_Ez.ptr(0);
572 
573  int ith, nth, is, ns;
574  set_threadtask(ith, nth, is, ns, m_Nvol);
575 
576 #pragma omp barrier
577 
578  for (int site = is; site < ns; ++site) {
579  int iv = Nvc * Nd * site;
580  int ig = Ndf * site;
581 
582  for (int ic = 0; ic < Nc; ++ic) {
583  int ic_r = 2 * ic;
584  int ic_i = ic_r + 1;
585  int ic_g = ic * Nvc + ig;
586 
587  v2[ic_r + id1 + iv] = 0.0;
588  v2[ic_i + id1 + iv] = 0.0;
589  v2[ic_r + id2 + iv] = 0.0;
590  v2[ic_i + id2 + iv] = 0.0;
591 
592  v2[ic_r + id3 + iv] = 0.0;
593  v2[ic_i + id3 + iv] = 0.0;
594  v2[ic_r + id4 + iv] = 0.0;
595  v2[ic_i + id4 + iv] = 0.0;
596 
597  // isigma_23 * Bx
598  v2[ic_r + id1 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id2 + iv], Nc);
599  v2[ic_i + id1 + iv] += mult_uv_r(&Bx[ic_g], &w2[id2 + iv], Nc);
600  v2[ic_r + id2 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id1 + iv], Nc);
601  v2[ic_i + id2 + iv] += mult_uv_r(&Bx[ic_g], &w2[id1 + iv], Nc);
602 
603  v2[ic_r + id3 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id4 + iv], Nc);
604  v2[ic_i + id3 + iv] += mult_uv_r(&Bx[ic_g], &w2[id4 + iv], Nc);
605  v2[ic_r + id4 + iv] -= mult_uv_i(&Bx[ic_g], &w2[id3 + iv], Nc);
606  v2[ic_i + id4 + iv] += mult_uv_r(&Bx[ic_g], &w2[id3 + iv], Nc);
607 
608  // isigma_31 * By
609  v2[ic_r + id1 + iv] += mult_uv_r(&By[ic_g], &w2[id2 + iv], Nc);
610  v2[ic_i + id1 + iv] += mult_uv_i(&By[ic_g], &w2[id2 + iv], Nc);
611  v2[ic_r + id2 + iv] -= mult_uv_r(&By[ic_g], &w2[id1 + iv], Nc);
612  v2[ic_i + id2 + iv] -= mult_uv_i(&By[ic_g], &w2[id1 + iv], Nc);
613 
614  v2[ic_r + id3 + iv] += mult_uv_r(&By[ic_g], &w2[id4 + iv], Nc);
615  v2[ic_i + id3 + iv] += mult_uv_i(&By[ic_g], &w2[id4 + iv], Nc);
616  v2[ic_r + id4 + iv] -= mult_uv_r(&By[ic_g], &w2[id3 + iv], Nc);
617  v2[ic_i + id4 + iv] -= mult_uv_i(&By[ic_g], &w2[id3 + iv], Nc);
618 
619  // isigma_12 * Bz
620  v2[ic_r + id1 + iv] -= mult_uv_i(&Bz[ic_g], &w2[id1 + iv], Nc);
621  v2[ic_i + id1 + iv] += mult_uv_r(&Bz[ic_g], &w2[id1 + iv], Nc);
622  v2[ic_r + id2 + iv] += mult_uv_i(&Bz[ic_g], &w2[id2 + iv], Nc);
623  v2[ic_i + id2 + iv] -= mult_uv_r(&Bz[ic_g], &w2[id2 + iv], Nc);
624 
625  v2[ic_r + id3 + iv] -= mult_uv_i(&Bz[ic_g], &w2[id3 + iv], Nc);
626  v2[ic_i + id3 + iv] += mult_uv_r(&Bz[ic_g], &w2[id3 + iv], Nc);
627  v2[ic_r + id4 + iv] += mult_uv_i(&Bz[ic_g], &w2[id4 + iv], Nc);
628  v2[ic_i + id4 + iv] -= mult_uv_r(&Bz[ic_g], &w2[id4 + iv], Nc);
629 
630  // isigma_41 * Ex
631  v2[ic_r + id1 + iv] += mult_uv_i(&Ex[ic_g], &w2[id2 + iv], Nc);
632  v2[ic_i + id1 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id2 + iv], Nc);
633  v2[ic_r + id2 + iv] += mult_uv_i(&Ex[ic_g], &w2[id1 + iv], Nc);
634  v2[ic_i + id2 + iv] -= mult_uv_r(&Ex[ic_g], &w2[id1 + iv], Nc);
635 
636  v2[ic_r + id3 + iv] -= mult_uv_i(&Ex[ic_g], &w2[id4 + iv], Nc);
637  v2[ic_i + id3 + iv] += mult_uv_r(&Ex[ic_g], &w2[id4 + iv], Nc);
638  v2[ic_r + id4 + iv] -= mult_uv_i(&Ex[ic_g], &w2[id3 + iv], Nc);
639  v2[ic_i + id4 + iv] += mult_uv_r(&Ex[ic_g], &w2[id3 + iv], Nc);
640 
641  // isigma_42 * Ey
642  v2[ic_r + id1 + iv] -= mult_uv_r(&Ey[ic_g], &w2[id2 + iv], Nc);
643  v2[ic_i + id1 + iv] -= mult_uv_i(&Ey[ic_g], &w2[id2 + iv], Nc);
644  v2[ic_r + id2 + iv] += mult_uv_r(&Ey[ic_g], &w2[id1 + iv], Nc);
645  v2[ic_i + id2 + iv] += mult_uv_i(&Ey[ic_g], &w2[id1 + iv], Nc);
646 
647  v2[ic_r + id3 + iv] += mult_uv_r(&Ey[ic_g], &w2[id4 + iv], Nc);
648  v2[ic_i + id3 + iv] += mult_uv_i(&Ey[ic_g], &w2[id4 + iv], Nc);
649  v2[ic_r + id4 + iv] -= mult_uv_r(&Ey[ic_g], &w2[id3 + iv], Nc);
650  v2[ic_i + id4 + iv] -= mult_uv_i(&Ey[ic_g], &w2[id3 + iv], Nc);
651 
652  // isigma_43 * Ez
653  v2[ic_r + id1 + iv] += mult_uv_i(&Ez[ic_g], &w2[id1 + iv], Nc);
654  v2[ic_i + id1 + iv] -= mult_uv_r(&Ez[ic_g], &w2[id1 + iv], Nc);
655  v2[ic_r + id2 + iv] -= mult_uv_i(&Ez[ic_g], &w2[id2 + iv], Nc);
656  v2[ic_i + id2 + iv] += mult_uv_r(&Ez[ic_g], &w2[id2 + iv], Nc);
657 
658  v2[ic_r + id3 + iv] -= mult_uv_i(&Ez[ic_g], &w2[id3 + iv], Nc);
659  v2[ic_i + id3 + iv] += mult_uv_r(&Ez[ic_g], &w2[id3 + iv], Nc);
660  v2[ic_r + id4 + iv] += mult_uv_i(&Ez[ic_g], &w2[id4 + iv], Nc);
661  v2[ic_i + id4 + iv] -= mult_uv_r(&Ez[ic_g], &w2[id4 + iv], Nc);
662 
663  // v *= m_kappa * m_cSW;
664  v2[ic_r + id1 + iv] *= kappa_cSW;
665  v2[ic_i + id1 + iv] *= kappa_cSW;
666  v2[ic_r + id2 + iv] *= kappa_cSW;
667  v2[ic_i + id2 + iv] *= kappa_cSW;
668 
669  v2[ic_r + id3 + iv] *= kappa_cSW;
670  v2[ic_i + id3 + iv] *= kappa_cSW;
671  v2[ic_r + id4 + iv] *= kappa_cSW;
672  v2[ic_i + id4 + iv] *= kappa_cSW;
673  }
674  }
675 
676 #pragma omp barrier
677  }
678 
679 
680 //====================================================================
681  void Fopr_CloverTerm::set_csw()
682  {
683  set_fieldstrength(m_Bx, 1, 2);
684  set_fieldstrength(m_By, 2, 0);
685  set_fieldstrength(m_Bz, 0, 1);
686  set_fieldstrength(m_Ex, 3, 0);
687  set_fieldstrength(m_Ey, 3, 1);
688  set_fieldstrength(m_Ez, 3, 2);
689  }
690 
691 
692 //====================================================================
693  void Fopr_CloverTerm::set_fieldstrength(Field_G& Fst,
694  const int mu, const int nu)
695  {
696 #pragma omp barrier
697 
698  m_staple.upper(m_Cup, *m_U, mu, nu); // these staple constructions
699  m_staple.lower(m_Cdn, *m_U, mu, nu); // are multi-threaded.
700 
701  mult_Field_Gnd(Fst, 0, *m_U, mu, m_Cup, 0);
702  multadd_Field_Gnd(Fst, 0, *m_U, mu, m_Cdn, 0, -1.0);
703 
704  mult_Field_Gdn(m_v1, 0, m_Cup, 0, *m_U, mu);
705  multadd_Field_Gdn(m_v1, 0, m_Cdn, 0, *m_U, mu, -1.0);
706 
707  m_shift->forward(m_v2, m_v1, mu);
708 
709  axpy(Fst, 1.0, m_v2);
710 #pragma omp barrier
711 
712  ah_Field_G(Fst, 0);
713 #pragma omp barrier
714 
715  scal(Fst, 0.25);
716 #pragma omp barrier
717  }
718 
719 
720 //====================================================================
721  double Fopr_CloverTerm::flop_count()
722  {
723  // Counting of floating point operations in giga unit.
724  // The following counting explicitly depends on the implementation
725  // and to be recalculated when the code is modified.
726  // Present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]
727 
728  const int Nvol = CommonParameters::Nvol();
729  const int NPE = CommonParameters::NPE();
730 
731  const int flop_site = m_Nc * m_Nd * (2 + 48 * m_Nc);
732 
733  const double gflop = flop_site * (Nvol * (NPE / 1.0e+9));
734 
735  return gflop;
736  }
737 
738 
739 //====================================================================
740 }
741 //============================================================END=====
GammaMatrixSet
Set of Gamma Matrices: basis class.
Definition: gammaMatrixSet.h:37
Imp::Fopr_CloverTerm::set_config_omp
void set_config_omp(Field *U)
Definition: fopr_CloverTerm_impl.cpp:236
Imp::Fopr_CloverTerm::mult_csw_dirac
void mult_csw_dirac(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:406
GammaMatrixSet::GAMMA5
@ GAMMA5
Definition: gammaMatrixSet.h:48
fopr_thread-inc.h
Imp::Fopr_CloverTerm::m_v2
Field_G m_v2
working vectors
Definition: fopr_CloverTerm_impl.h:87
mult_Field_Gdn
void mult_Field_Gdn(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:134
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
Imp::Fopr_CloverTerm::m_mode
std::string m_mode
Definition: fopr_CloverTerm_impl.h:68
Imp::Fopr_CloverTerm::set_csw
void set_csw()
Definition: fopr_CloverTerm_impl.cpp:681
fopr_Wilson_impl_SU_N-inc.h
fopr_Wilson_impl_SU3-inc.h
ShiftField_lex::forward
void forward(Field &, const Field &, const int mu)
Definition: shiftField_lex.cpp:79
ThreadManager::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager.cpp:246
Imp::Fopr_CloverTerm::set_mode
void set_mode(const std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: fopr_CloverTerm_impl.cpp:261
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Parameters
Class for parameters.
Definition: parameters.h:46
Imp::Fopr_CloverTerm::mult_sigmaF
void mult_sigmaF(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:382
Imp::Fopr_CloverTerm::m_By
Field_G m_By
Definition: fopr_CloverTerm_impl.h:81
GammaMatrixSet::UNITY
@ UNITY
Definition: gammaMatrixSet.h:48
Imp::Fopr_CloverTerm::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: fopr_CloverTerm_impl.cpp:219
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.h:86
Imp::Fopr_CloverTerm::m_staple
Staple_lex m_staple
Definition: fopr_CloverTerm_impl.h:79
Imp::Fopr_CloverTerm::class_name
static const std::string class_name
Definition: fopr_CloverTerm_impl.h:58
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:219
Field::nex
int nex() const
Definition: field.h:128
Imp::Fopr_CloverTerm::mult_csw
void mult_csw(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:392
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Staple_lex::lower
void lower(Field_G &, const Field_G &, const int mu, const int nu)
constructs lower staple in mu-nu plane.
Definition: staple_lex.cpp:286
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
Imp::Fopr_CloverTerm::mult_isigma
void mult_isigma(Field_F &, const Field_F &, const int mu, const int nu)
Definition: fopr_CloverTerm_impl.cpp:373
Imp::Fopr_CloverTerm::m_Cdn
Field_G m_Cdn
Definition: fopr_CloverTerm_impl.h:87
Imp::Fopr_CloverTerm::set_config_impl
void set_config_impl(Field *U)
Definition: fopr_CloverTerm_impl.cpp:246
Imp::Fopr_CloverTerm::m_shift
ShiftField_lex * m_shift
Definition: fopr_CloverTerm_impl.h:78
Imp::Fopr_CloverTerm::gm5_dirac
void gm5_dirac(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:303
GammaMatrixSet::SIGMA41
@ SIGMA41
Definition: gammaMatrixSet.h:52
Imp::Fopr_CloverTerm::flop_count
double flop_count()
this returns the number of floating point operations.
Definition: fopr_CloverTerm_impl.cpp:721
Imp::Fopr_CloverTerm::m_vl
Bridge::VerboseLevel m_vl
Definition: fopr_CloverTerm_impl.h:61
Imp::Fopr_CloverTerm::m_Ex
Field_G m_Ex
Definition: fopr_CloverTerm_impl.h:84
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Imp::Fopr_CloverTerm::m_Ndim
int m_Ndim
Definition: fopr_CloverTerm_impl.h:73
Imp::Fopr_CloverTerm::m_Nvol
int m_Nvol
Definition: fopr_CloverTerm_impl.h:74
Imp::Fopr_CloverTerm::sg_index
int sg_index(const int mu, const int nu)
Definition: fopr_CloverTerm_impl.h:164
ah_Field_G
void ah_Field_G(Field_G &W, const int ex)
Definition: field_G_imp.cpp:462
mult_iGM
void mult_iGM(Field_F &y, const GammaMatrix &gm, const Field_F &x)
gamma matrix multiplication (i is multiplied)
Definition: field_F_imp.cpp:250
Imp::Fopr_CloverTerm::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_CloverTerm_impl.cpp:176
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
Imp::Fopr_CloverTerm::m_cSW
double m_cSW
Definition: fopr_CloverTerm_impl.h:64
Imp::Fopr_CloverTerm::m_U
const Field_G * m_U
pointer to gauge configuration.
Definition: fopr_CloverTerm_impl.h:76
Imp::Fopr_CloverTerm::m_repr
std::string m_repr
Definition: fopr_CloverTerm_impl.h:66
Imp::Fopr_CloverTerm::m_Cup
Field_G m_Cup
Definition: fopr_CloverTerm_impl.h:87
Field::nvol
int nvol() const
Definition: field.h:127
Imp::Fopr_CloverTerm::m_NinF
int m_NinF
Definition: fopr_CloverTerm_impl.h:73
Imp::Fopr_CloverTerm::init
void init(const std::string repr)
Definition: fopr_CloverTerm_impl.cpp:78
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Imp::Fopr_CloverTerm::gm5_chiral
void gm5_chiral(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:338
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
GammaMatrixSet::SIGMA42
@ SIGMA42
Definition: gammaMatrixSet.h:52
Imp::Fopr_CloverTerm::setup_gamma_matrices
void setup_gamma_matrices()
Definition: fopr_CloverTerm_impl.cpp:188
ShiftField_lex
Methods to shift a field in the lexical site index.
Definition: shiftField_lex.h:39
multadd_Field_Gnd
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)
Definition: field_G_imp.cpp:335
Imp::Fopr_CloverTerm::m_kappa
double m_kappa
Definition: fopr_CloverTerm_impl.h:63
Imp::Fopr_CloverTerm::mult_gm5
void mult_gm5(Field &v, const Field &w)
multiplies gamma_5 matrix.
Definition: fopr_CloverTerm_impl.cpp:292
Imp::Fopr_CloverTerm::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: fopr_CloverTerm_impl.cpp:118
Imp::Fopr_CloverTerm::m_Bx
Field_G m_Bx
Definition: fopr_CloverTerm_impl.h:81
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
Staple_lex::upper
void upper(Field_G &, const Field_G &, const int mu, const int nu)
constructs upper staple in mu-nu plane.
Definition: staple_lex.cpp:263
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
fopr_Wilson_impl_SU2-inc.h
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
Imp::Fopr_CloverTerm::mult_dag
void mult_dag(Field &v, const Field &f)
hermitian conjugate of mult.
Definition: fopr_CloverTerm_impl.cpp:285
Imp::Fopr_CloverTerm::m_Nd
int m_Nd
Definition: fopr_CloverTerm_impl.h:73
Imp::Fopr_CloverTerm::m_boundary
std::vector< int > m_boundary
Definition: fopr_CloverTerm_impl.h:65
GammaMatrixSet::SIGMA23
@ SIGMA23
Definition: gammaMatrixSet.h:51
Imp::Fopr_CloverTerm::m_GM5
GammaMatrix m_GM5
Definition: fopr_CloverTerm_impl.h:90
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
Field_F
Wilson-type fermion field.
Definition: field_F.h:37
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
GammaMatrixSet::SIGMA12
@ SIGMA12
Definition: gammaMatrixSet.h:51
Imp::Fopr_CloverTerm::mult
void mult(Field &v, const Field &f)
multiplies fermion operator to a given field.
Definition: fopr_CloverTerm_impl.cpp:271
Imp::Fopr_CloverTerm::mult_csw_chiral
void mult_csw_chiral(Field &, const Field &)
Definition: fopr_CloverTerm_impl.cpp:544
GammaMatrixSet::get_GM
GammaMatrix get_GM(GMspecies spec)
Definition: gammaMatrixSet.h:76
Imp::Fopr_CloverTerm::m_v1
Field_G m_v1
Definition: fopr_CloverTerm_impl.h:87
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:180
Imp
Clover term operator.
Definition: fopr_CloverTerm_eo_impl.cpp:32
Imp::Fopr_CloverTerm::m_Ey
Field_G m_Ey
Definition: fopr_CloverTerm_impl.h:84
Field
Container of Field-type object.
Definition: field.h:46
GammaMatrixSet::SIGMA31
@ SIGMA31
Definition: gammaMatrixSet.h:51
Imp::Fopr_CloverTerm::m_Ez
Field_G m_Ez
field strength (electric components)
Definition: fopr_CloverTerm_impl.h:84
Imp::Fopr_CloverTerm::tidyup
void tidyup()
Definition: fopr_CloverTerm_impl.cpp:111
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
GammaMatrixSet::SIGMA43
@ SIGMA43
Definition: gammaMatrixSet.h:52
fopr_CloverTerm_impl.h
Imp::Fopr_CloverTerm::set_fieldstrength
void set_fieldstrength(Field_G &, const int, const int)
Definition: fopr_CloverTerm_impl.cpp:693
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Imp::Fopr_CloverTerm::m_Nc
int m_Nc
Definition: fopr_CloverTerm_impl.h:73
Imp::Fopr_CloverTerm::m_Bz
Field_G m_Bz
field strength (magnetic components)
Definition: fopr_CloverTerm_impl.h:81
Imp::Fopr_CloverTerm::m_SG
std::vector< GammaMatrix > m_SG
Definition: fopr_CloverTerm_impl.h:89
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
multadd_Field_Gdn
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)
Definition: field_G_imp.cpp:293
mult_Field_Gnd
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:173
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154