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