Bridge++  Ver.2.1.3
afopr_CloverTerm-tmpl.h
Go to the documentation of this file.
1 
10 template<typename AFIELD>
11 const std::string AFopr_CloverTerm<AFIELD>::class_name
12  = "AFopr_CloverTerm<AFIELD>";
13 
14 //====================================================================
15 template<typename AFIELD>
16 void AFopr_CloverTerm<AFIELD>::init(const Parameters& params)
17 {
18  ThreadManager::assert_single_thread(class_name);
19 
20  std::string vlevel;
21  if (!params.fetch_string("verbose_level", vlevel)) {
22  m_vl = vout.set_verbose_level(vlevel);
23  } else {
24  m_vl = CommonParameters::Vlevel();
25  }
26 
27  vout.general(m_vl, "%s: construction\n", class_name.c_str());
28  vout.increase_indent();
29 
30  m_Nc = CommonParameters::Nc();
31  m_Nd = CommonParameters::Nd();
32  m_Ndim = CommonParameters::Ndim();
33  m_Nvc = m_Nc * 2;
34  m_Ndf = 2 * m_Nc * m_Nc;
35  m_Ndm2 = m_Nd * m_Nd / 2,
36 
37  m_Nx = CommonParameters::Nx();
38  m_Ny = CommonParameters::Ny();
39  m_Nz = CommonParameters::Nz();
40  m_Nt = CommonParameters::Nt();
41  m_Nst = CommonParameters::Nvol();
42 
43  m_Nsize[0] = m_Nx;
44  m_Nsize[1] = m_Ny;
45  m_Nsize[2] = m_Nz;
46  m_Nsize[3] = m_Nt;
47 
48  set_parameters(params);
49 
50  m_staple = new AStaple_lex<AFIELD>;
51 
52  // gauge configuration.
53  m_U.reset(m_Ndf, m_Nst, m_Ndim);
54 
55  // clover term.
56  m_T.reset(m_Ndf, m_Nst, m_Ndm2);
57  m_Tinv.reset(m_Ndf, m_Nst, m_Ndm2);
58 
59  // working vectors.
60  int NinF = 2 * m_Nc * m_Nd;
61  m_v1.reset(NinF, m_Nst, 1);
62  m_v2.reset(NinF, m_Nst, 1);
63 
64  // working gauge field.
65  m_ut1.reset(m_Ndf, m_Nst, 1), m_ut2.reset(m_Ndf, m_Nst, 1);
66 
67  m_U2.reset(m_Ndf, m_Nst, 4);
68  m_F2.reset(m_Ndf, m_Nst, 1);
69 
70  // setup solver.
71  double ecrit = 1.0e-30;
72  if(sizeof(real_t) == 4) ecrit = 1.0e-16;
73 
74  Parameters params_solver;
75  params_solver.set_string("solver_type", "CG");
76  params_solver.set_int("maximum_number_of_iteration", 100);
77  params_solver.set_int("maximum_number_of_restart", 10);
78  params_solver.set_double("convergence_criterion_squared", ecrit);
79  //- NB. set VerboseLevel to CRUCIAL to suppress frequent messages.
80  // params_solver.set_string("verbose_level", "Detailed");
81  params_solver.set_string("verbose_level", vlevel);
82 
83  m_solver = new ASolver_CG<AFIELD>(this);
84  m_solver->set_parameters(params_solver);
85 
86  vout.decrease_indent();
87  vout.detailed(m_vl, "%s: construction finished.\n",
88  class_name.c_str());
89 
90 }
91 
92 //====================================================================
93 template<typename AFIELD>
94 void AFopr_CloverTerm<AFIELD>::tidyup()
95 {
96  ThreadManager::assert_single_thread(class_name);
97 
98  delete m_staple;
99  delete m_solver;
100 
101 }
102 
103 //====================================================================
104 template<typename AFIELD>
105 void AFopr_CloverTerm<AFIELD>::set_parameters(const Parameters& params)
106 {
107  std::string vlevel;
108  if (!params.fetch_string("verbose_level", vlevel)) {
109  m_vl = vout.set_verbose_level(vlevel);
110  }
111 
112  //- fetch and check input parameters
113  double kappa, cSW;
114  std::vector<int> bc;
115 
116  int err = 0;
117  err += params.fetch_double("hopping_parameter", kappa);
118  err += params.fetch_double("clover_coefficient", cSW);
119 
120  err += params.fetch_int_vector("boundary_condition", bc);
121  if (err) {
122  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
123  class_name.c_str());
124  exit(EXIT_FAILURE);
125  }
126 
127  //- setting gamma matrix representation
128  std::string repr;
129  err = params.fetch_string("gamma_matrix_type", repr);
130  if(err){
131  vout.general(m_vl, " gamma_matrix_type is not given - set to Dirac\n");
132  m_repr = DIRAC;
133  }else if(repr == "Dirac"){
134  m_repr = DIRAC;
135  }else if(repr == "Chiral"){
136  m_repr = CHIRAL;
137  }else{
138  vout.crucial(m_vl, "Error in %s: irrelevant gamma_matrix_type: %s\n",
139  class_name.c_str(), repr.c_str());
140  exit(EXIT_FAILURE);
141  }
142 
143  set_parameters(real_t(kappa), real_t(cSW), bc);
144 
145 }
146 
147 //====================================================================
148 template<typename AFIELD>
149 void AFopr_CloverTerm<AFIELD>::set_parameters(const real_t CKs,
150  const real_t cSW,
151  const std::vector<int> bc)
152 {
153  assert(bc.size() == m_Ndim);
154 
155 #pragma omp barrier
156 
157  int ith = ThreadManager::get_thread_id();
158 
159  if (ith == 0) {
160  m_CKs = CKs;
161  m_cSW = cSW;
162  m_boundary.resize(m_Ndim);
163  for (int mu = 0; mu < m_Ndim; ++mu) {
164  m_boundary[mu] = bc[mu];
165  }
166  }
167 
168  //- print input parameters
169  vout.general(m_vl, "Parameters of %s:\n", class_name.c_str());
170  if(m_repr == DIRAC){
171  vout.general(m_vl, " gamma-matrix type = Dirac\n");
172  }else{
173  vout.general(m_vl, " gamma-matrix type = Chiral\n");
174  }
175  vout.general(m_vl, " kappa = %8.4f\n", m_CKs);
176  vout.general(m_vl, " cSW = %8.4f\n", m_cSW);
177  for (int mu = 0; mu < m_Ndim; ++mu) {
178  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
179  }
180 
181 #pragma omp barrier
182 }
183 
184 //====================================================================
185 template<typename AFIELD>
186 void AFopr_CloverTerm<AFIELD>::set_config(Field* u)
187 {
188  int nth = ThreadManager::get_num_threads();
189  if (nth > 1) {
190  set_config_impl(u);
191  } else {
192  set_config_omp(u);
193  }
194 
195 }
196 
197 //====================================================================
198 template<typename AFIELD>
199 void AFopr_CloverTerm<AFIELD>::set_config_omp(Field* u)
200 {
201 #pragma omp parallel
202  {
203  set_config_impl(u);
204  }
205 }
206 
207 //====================================================================
208 template<typename AFIELD>
209 void AFopr_CloverTerm<AFIELD>::set_config_impl(Field* u)
210 {
211 #pragma omp barrier
212 
213  vout.detailed(m_vl, "%s: set_config started.\n", class_name.c_str());
214 
215  Timer timer;
216  double elapsed_time;
217  timer.start();
218 
219  int ith = ThreadManager::get_thread_id();
220  if (ith == 0) m_conf = u;
221 
222  AIndex_lex<real_t,AFIELD::IMPL> index_lex;
223 
224  convert_gauge(index_lex, m_U, *u);
225 
226  timer.stop();
227  elapsed_time = timer.elapsed_sec();
228  vout.detailed(m_vl, " convert: %11.6f [sec]\n", elapsed_time);
229 
230  timer.reset();
231  timer.start();
232 
233  set_csw(*u);
234 
235  timer.stop();
236  elapsed_time = timer.elapsed_sec();
237  vout.detailed(m_vl, " set_csw: %11.6f [sec]\n", elapsed_time);
238 
239  timer.reset();
240  timer.start();
241 
242  solve_csw_inv();
243 
244  timer.stop();
245  elapsed_time = timer.elapsed_sec();
246  vout.detailed(m_vl, " set_csw_inv: %11.6f [sec]\n", elapsed_time);
247 
248  vout.detailed(m_vl, "%s: set_config finished.\n", class_name.c_str());
249 
250 #pragma omp barrier
251 }
252 
253 //====================================================================
254 template<typename AFIELD>
255 void AFopr_CloverTerm<AFIELD>::set_csw(Field& U)
256 {
257  if(m_repr == DIRAC){
258  set_csw_dirac(U);
259  }else if(m_repr == CHIRAL){
260  set_csw_chiral(U);
261  }else{
262  vout.crucial(m_vl, "%s: unsupported representation.\n",
263  class_name.c_str());
264  exit(EXIT_FAILURE);
265  }
266 
267 }
268 
269 //====================================================================
270 template<typename AFIELD>
271 void AFopr_CloverTerm<AFIELD>::solve_csw_inv()
272 {
273  vout.paranoiac(m_vl, " %s: solving inverse of clover term started.\n",
274  class_name.c_str());
275 
276  if(m_repr == DIRAC){
277  solve_csw_inv_dirac();
278  }else{
279  solve_csw_inv_chiral();
280  }
281 
282  vout.paranoiac(m_vl, " %s: solving inverse of clover term finished.\n",
283  class_name.c_str());
284 
285  // check
286  /*
287  vout.general(m_vl, " %s: check of clover term inverse\n",
288  class_name.c_str());
289 
290  int NinF = 2 * m_Nc * m_Nd;
291  AFIELD w(NinF, m_Nst, 1), w2(NinF, m_Nst, 1);
292  AFIELD w3(NinF, m_Nst, 1);
293 
294  AIndex_lex<real_t> index;
295 
296  for(int id = 0; id < m_Nd; ++id){
297  for(int ic = 0; ic < m_Nc; ++ic){
298 
299  w.set_host(0.0);
300  for(int site = 0; site < m_Nst; ++site){
301  w.set_host(index.idx_SPr(ic, id, site, 0), real_t(1.0));
302  }
303  w.update_device();
304 
305  mult_csw(w3, w);
306  mult_csw_inv(w2, w3);
307  axpy(w2, real_t(-1.0), w);
308  real_t ww = w2.norm2();
309 
310  vout.general(m_vl, " ic = %d id = %d diff2 = %f\n",
311  ic, id, ww);
312 
313  }
314  }
315  */
316 
317 }
318 
319 //====================================================================
320 template<typename AFIELD>
321 void AFopr_CloverTerm<AFIELD>::solve_csw_inv_dirac()
322 {
323 #pragma omp barrier
324 
325  set_mode("D");
326 
327  int Nconv;
328  real_t diff;
329  AIndex_lex<real_t,AFIELD::IMPL> index;
330 
331  m_Tinv.set_host(0.0);
332 
333  int ith, nth, is, ns;
334  set_threadtask(ith, nth, is, ns, m_Nst);
335 
336  int Nd2 = m_Nd/2;
337  for(int id = 0; id < Nd2; ++id){
338  for(int ic = 0; ic < m_Nc; ++ic){
339 
340  m_v1.set_host(0.0);
341  for(int site = is; site < ns; ++site){
342  m_v1.set_host(index.idx_SPr(ic, id, site, 0), real_t(1.0));
343  }
344 #pragma omp barrier
345 
346  update_device(m_v1);
347  copy(m_v2, m_v1);
348 
349  m_solver->solve(m_v2, m_v1, Nconv, diff);
350  vout.paranoiac(m_vl, " ic = %d id = %d Nconv = %d diff = %12.4e\n",
351  ic, id, Nconv, diff);
352 
353  update_host(m_v2);
354 
355  for(int site = is; site < ns; ++site){
356  for(int ic2 = 0; ic2 < m_Nc; ++ic2){
357  for(int id2 = 0; id2 < m_Nd; ++id2){
358  real_t re = m_v2.cmp_host(index.idx_SPr(ic2, id2, site, 0));
359  real_t im = m_v2.cmp_host(index.idx_SPi(ic2, id2, site, 0));
360  int iT = id2 + m_Nd * id;
361  m_Tinv.set_host(index.idx_Gr(ic2, ic, site, iT), re);
362  m_Tinv.set_host(index.idx_Gi(ic2, ic, site, iT), -im);
363  }
364  }
365  }
366 #pragma omp barrier
367 
368  }
369  }
370 
371  update_device(m_Tinv);
372 
373 #pragma omp barrier
374 
375 }
376 
377 //====================================================================
378 template<typename AFIELD>
379 void AFopr_CloverTerm<AFIELD>::solve_csw_inv_chiral()
380 {
381 #pragma omp barrier
382 
383  set_mode("D");
384 
385  int Nconv;
386  real_t diff;
387  AIndex_lex<real_t,AFIELD::IMPL> index;
388 
389  m_Tinv.set_host(0.0);
390 
391  int ith, nth, is, ns;
392  set_threadtask(ith, nth, is, ns, m_Nst);
393 
394  int Nd2 = m_Nd/2;
395  for(int id = 0; id < Nd2; ++id){
396  for(int ic = 0; ic < m_Nc; ++ic){
397 
398  m_v1.set_host(0.0);
399  for(int site = is; site < ns; ++site){
400  m_v1.set_host(index.idx_SPr(ic, id, site, 0), real_t(1.0));
401  m_v1.set_host(index.idx_SPr(ic, id+Nd2, site, 0), real_t(1.0));
402  }
403 #pragma omp barrier
404 
405  update_device(m_v1);
406  copy(m_v2, m_v1);
407 
408  m_solver->solve(m_v2, m_v1, Nconv, diff);
409  vout.paranoiac(m_vl, " ic = %d id = %d Nconv = %d diff = %12.4e\n",
410  ic, id, Nconv, diff);
411 
412  update_host(m_v2);
413 
414  for(int site = is; site < ns; ++site){
415 
416  for(int ic2 = 0; ic2 < m_Nc; ++ic2){
417  for(int id2 = 0; id2 < Nd2; ++id2){
418  real_t re = m_v2.cmp_host(index.idx_SPr(ic2, id2, site, 0));
419  real_t im = m_v2.cmp_host(index.idx_SPi(ic2, id2, site, 0));
420  int iT = id2 + Nd2 * id;
421  m_Tinv.set_host(index.idx_Gr(ic2, ic, site, iT), re);
422  m_Tinv.set_host(index.idx_Gi(ic2, ic, site, iT), -im);
423  }
424  }
425 
426  for(int ic2 = 0; ic2 < m_Nc; ++ic2){
427  for(int id2 = 0; id2 < Nd2; ++id2){
428  int jd2 = id2 + Nd2;
429  real_t re = m_v2.cmp_host(index.idx_SPr(ic2, jd2, site, 0));
430  real_t im = m_v2.cmp_host(index.idx_SPi(ic2, jd2, site, 0));
431  int iT = id2 + Nd2 * id + ND;
432  m_Tinv.set_host(index.idx_Gr(ic2, ic, site, iT), re);
433  m_Tinv.set_host(index.idx_Gi(ic2, ic, site, iT), -im);
434  }
435  }
436 
437  }
438 #pragma omp barrier
439 
440  }
441  }
442 
443  update_device(m_Tinv);
444 
445 #pragma omp barrier
446 
447 }
448 
449 //====================================================================
450 template<typename AFIELD>
451 void AFopr_CloverTerm<AFIELD>::get_csw(AFIELD& T)
452 {
453  if(T.check_size(m_T) == false){
454  vout.crucial(m_vl, "%s: in get_csw, incorrect AFIELD size.\n",
455  class_name.c_str());
456  exit(EXIT_FAILURE);
457  }
458 
459  copy(T, m_T);
460 
461 }
462 
463 //====================================================================
464 template<typename AFIELD>
465 void AFopr_CloverTerm<AFIELD>::get_csw_inv(AFIELD& T)
466 {
467  if(T.check_size(m_T) == false){
468  vout.crucial(m_vl, "%s: in get_csw_inv, incorrect AFIELD size.\n",
469  class_name.c_str());
470  exit(EXIT_FAILURE);
471  }
472 
473  // the following is temporary prescription.
474  // to be performed in set_config. [H.Matsufuru 16 May 2021]
475  Timer timer;
476 
477  timer.start();
478 
479  solve_csw_inv();
480 
481  timer.stop();
482  double elapsed_time = timer.elapsed_sec();
483  vout.detailed(m_vl, "%s: set_csw_inv: %11.6f [sec]\n",
484  class_name.c_str(), elapsed_time);
485 
486  copy(T, m_Tinv);
487 #pragma omp barrier
488 
489 }
490 
491 //====================================================================
492 template<typename AFIELD>
493 void AFopr_CloverTerm<AFIELD>::get_csw_inv(AFIELD& T, const int j)
494 {
495  int Ndf = m_T.nin();
496  int Nst = m_T.nvol();
497  int Nex = m_T.nex();
498  if(T.check_size(Ndf, Nst, 1) == false){
499  vout.crucial(m_vl, "%s: in get_csw_inv, incorrect AFIELD size.\n",
500  class_name.c_str());
501  exit(EXIT_FAILURE);
502  }
503 
504  if(j >= Nex){
505  vout.crucial(m_vl, "%s: in get_csw_inv, index is too large.\n",
506  class_name.c_str());
507  exit(EXIT_FAILURE);
508  }
509 
510  copy(T, 0, m_Tinv, j);
511 #pragma omp barrier
512 
513 }
514 
515 //====================================================================
516 template<typename AFIELD>
517 void AFopr_CloverTerm<AFIELD>::set_csw_dirac(Field& U)
518 { // this is for Dirac representation.
519 
520  // The clover term in the Dirac representation is as spin-space
521  // matrix
522  // [ P Q ]
523  // [ Q P ],
524  // where P and Q are 2x2 block matrices as
525  // P = [ iF(1,2) F(3,1) + iF(2,3) ]
526  // [-F(3,1) + iF(2,3) - iF(1,2) ]
527  // and
528  // Q = [ - iF(4,3) -F(4,2) - iF(4,1) ]
529  // [ F(4,2) - iF(4,1) iF(4,3) ]
530  // up to the coefficient.
531  // in the following what defined is
532  // [ P Q ] = [ T(0) T(1) T(2) T(3) ]
533  // [ T(4) T(5) T(6) T(7) ].
534 
536 
537 #pragma omp barrier
538 
539  AIndex_lex<real_t,AFIELD::IMPL> index_lex;
540 
541  convert_gauge(index_lex, m_U2, U);
542 
543  m_T.set(0.0);
544 #pragma omp barrier
545 
546  //- sigma23
547  set_fieldstrength(m_F2, m_U2, 1, 2);
548  xI(m_F2);
549 
550 #pragma omp barrier
551  axpy(m_T, 1, real_t(1.0), m_F2, 0);
552  axpy(m_T, 4, real_t(1.0), m_F2, 0);
553 #pragma omp barrier
554 
555  //- sigma31
556  set_fieldstrength(m_F2, m_U2, 2, 0);
557 #pragma omp barrier
558  axpy(m_T, 1, real_t( 1.0), m_F2, 0);
559  axpy(m_T, 4, real_t(-1.0), m_F2, 0);
560 #pragma omp barrier
561 
562  //- sigma12
563  set_fieldstrength(m_F2, m_U2, 0, 1);
564  xI(m_F2);
565 #pragma omp barrier
566  axpy(m_T, 0, real_t( 1.0), m_F2, 0);
567  axpy(m_T, 5, real_t(-1.0), m_F2, 0);
568 #pragma omp barrier
569 
570  //- sigma41
571  set_fieldstrength(m_F2, m_U2, 3, 0);
572  xI(m_F2);
573 #pragma omp barrier
574  axpy(m_T, 3, real_t(-1.0), m_F2, 0);
575  axpy(m_T, 6, real_t(-1.0), m_F2, 0);
576 #pragma omp barrier
577 
578  //- sigma42
579  set_fieldstrength(m_F2, m_U2, 3, 1);
580 #pragma omp barrier
581  axpy(m_T, 3, real_t(-1.0), m_F2, 0);
582  axpy(m_T, 6, real_t( 1.0), m_F2, 0);
583 #pragma omp barrier
584 
585  //- sigma43
586  set_fieldstrength(m_F2, m_U2, 3, 2);
587  xI(m_F2);
588 #pragma omp barrier
589  axpy(m_T, 2, real_t(-1.0), m_F2, 0);
590  axpy(m_T, 7, real_t( 1.0), m_F2, 0);
591 #pragma omp barrier
592 
593  scal(m_T, -m_CKs * m_cSW);
594 
595 #pragma omp barrier
596 
597  add_unit(m_T, 0, real_t(1.0));
598  add_unit(m_T, 5, real_t(1.0));
599 
600 #pragma omp barrier
601 
602  // check of norm
603  //real_t tt = m_T.norm2();
604  //vout.crucial("norm of T = %e\n", tt);
605 
606 #pragma omp barrier
607 
608 }
609 
610 //====================================================================
611 template<typename AFIELD>
612 void AFopr_CloverTerm<AFIELD>::set_csw_chiral(Field& U)
613 { // this is for Dirac representation.
614 
615  // The clover term in the chiral representation is as spin-space
616  // matrix
617  // [ P+Q 0 ]
618  // [ 0 P-Q ],
619  // where P and Q are 2x2 block matrices as
620  // P = [ iF(1,2) F(3,1) + iF(2,3) ]
621  // [-F(3,1) + iF(2,3) - iF(1,2) ]
622  // and
623  // Q = [ - iF(4,3) -F(4,2) - iF(4,1) ]
624  // [ F(4,2) - iF(4,1) iF(4,3) ]
625  // up to the coefficient.
626  // in the following what defined is
627  // [ T(0) T(1) ] = P + Q [ T(4) T(5) ] = P - Q
628  // [ T(2) T(3) ] [ T(6) T(7) ]
629  // T = 1 - kappa c_SW sigma F / 2
630 
631  // AFIELD U2(m_Ndf, m_Nst, 4), F2(m_Ndf, m_Nst, 1);
632 
633 #pragma omp barrier
634 
635  AIndex_lex<real_t,AFIELD::IMPL> index_lex;
636  // convert_gauge(index_lex, U2, U);
637 
638  //#pragma omp parallel
639  {
640  convert_gauge(index_lex, m_U2, U);
641 
642  m_T.set(0.0);
643  #pragma omp barrier
644 
645  //- sigma23
646  set_fieldstrength(m_F2, m_U2, 1, 2); // F_23
647  xI(m_F2);
648 #pragma omp barrier
649  axpy(m_T, 1, real_t(1.0), m_F2, 0);
650  axpy(m_T, 2, real_t(1.0), m_F2, 0);
651  axpy(m_T, 5, real_t(1.0), m_F2, 0);
652  axpy(m_T, 6, real_t(1.0), m_F2, 0);
653 #pragma omp barrier
654 
655  //- sigma31
656  set_fieldstrength(m_F2, m_U2, 2, 0); // F_31
657 #pragma omp barrier
658  axpy(m_T, 1, real_t( 1.0), m_F2, 0);
659  axpy(m_T, 2, real_t(-1.0), m_F2, 0);
660  axpy(m_T, 5, real_t( 1.0), m_F2, 0);
661  axpy(m_T, 6, real_t(-1.0), m_F2, 0);
662 #pragma omp barrier
663 
664  //- sigma12
665  set_fieldstrength(m_F2, m_U2, 0, 1); // F_12
666  xI(m_F2);
667 #pragma omp barrier
668  axpy(m_T, 0, real_t( 1.0), m_F2, 0);
669  axpy(m_T, 3, real_t(-1.0), m_F2, 0);
670  axpy(m_T, 4, real_t( 1.0), m_F2, 0);
671  axpy(m_T, 7, real_t(-1.0), m_F2, 0);
672 #pragma omp barrier
673 
674  //- sigma41
675  set_fieldstrength(m_F2, m_U2, 3, 0); // F_41
676  xI(m_F2);
677 #pragma omp barrier
678  axpy(m_T, 1, real_t(-1.0), m_F2, 0);
679  axpy(m_T, 2, real_t(-1.0), m_F2, 0);
680  axpy(m_T, 5, real_t( 1.0), m_F2, 0);
681  axpy(m_T, 6, real_t( 1.0), m_F2, 0);
682 #pragma omp barrier
683 
684  //- sigma42
685  set_fieldstrength(m_F2, m_U2, 3, 1); // F_42
686 #pragma omp barrier
687  axpy(m_T, 1, real_t(-1.0), m_F2, 0);
688  axpy(m_T, 2, real_t( 1.0), m_F2, 0);
689  axpy(m_T, 5, real_t( 1.0), m_F2, 0);
690  axpy(m_T, 6, real_t(-1.0), m_F2, 0);
691 #pragma omp barrier
692 
693  //- sigma43
694  set_fieldstrength(m_F2, m_U2, 3, 2); // F_43
695  xI(m_F2);
696 #pragma omp barrier
697  axpy(m_T, 0, real_t(-1.0), m_F2, 0);
698  axpy(m_T, 3, real_t( 1.0), m_F2, 0);
699  axpy(m_T, 4, real_t( 1.0), m_F2, 0);
700  axpy(m_T, 7, real_t(-1.0), m_F2, 0);
701 #pragma omp barrier
702 
703  scal(m_T, -m_CKs * m_cSW);
704 
705 #pragma omp barrier
706 
707  add_unit(m_T, 0, real_t(1.0));
708  add_unit(m_T, 3, real_t(1.0));
709  add_unit(m_T, 4, real_t(1.0));
710  add_unit(m_T, 7, real_t(1.0));
711 
712 #pragma omp barrier
713  }
714 
715 }
716 
717 //====================================================================
718 template<typename AFIELD>
719 void AFopr_CloverTerm<AFIELD>::set_fieldstrength(AFIELD& Fst, AFIELD& U,
720  int mu, int nu)
721 {
722  m_staple->upper(m_ut2, U, mu, nu);
723 
724  mult_Gnd(Fst, 0, U, mu, m_ut2, 0);
725  mult_Gdn(m_ut1, 0, m_ut2, 0, U, mu);
726 
727  m_staple->lower(m_ut2, U, mu, nu);
728 
729  multadd_Gnd(Fst, 0, U, mu, m_ut2, 0, real_t(-1.0));
730  multadd_Gdn(m_ut1, 0, m_ut2, 0, U, mu, real_t(-1.0));
731 
732  m_staple->shift_forward(m_ut2, 0, m_ut1, 0, mu);
733 
734  axpy(Fst, real_t(1.0), m_ut2);
735 
736 #pragma omp barrier
737  ah_G(Fst, 0);
738 
739 #pragma omp barrier
740  scal(Fst, real_t(0.25));
741 
742 #pragma omp barrier
743 
744 }
745 
746 //====================================================================
747 template<typename AFIELD>
748 void AFopr_CloverTerm<AFIELD>::set_mode(std::string mode)
749 {
750 #pragma omp barrier
751 
752  int ith = ThreadManager::get_thread_id();
753  if (ith == 0) m_mode = mode;
754 
755 #pragma omp barrier
756 
757 }
758 
759 //====================================================================
760 template<typename AFIELD>
761 std::string AFopr_CloverTerm<AFIELD>::get_mode() const
762 {
763  return m_mode;
764 }
765 
766 //====================================================================
767 template<typename AFIELD>
768 void AFopr_CloverTerm<AFIELD>::mult(AFIELD &v, const AFIELD &w)
769 {
770  if(m_mode == "D"){
771  D(v, w);
772  }else if(m_mode == "Dinv"){
773  mult_csw_inv(v, w);
774  }else if(m_mode == "H"){
775  return H(v, w);
776  }else{
777  vout.crucial(m_vl, "%s: mode undefined.\n", class_name.c_str());
778  exit(EXIT_FAILURE);
779  }
780 
781 }
782 
783 //====================================================================
784 template<typename AFIELD>
785 void AFopr_CloverTerm<AFIELD>::mult_dag(AFIELD &v, const AFIELD &w)
786 {
787  if(m_mode == "D"){
788  D(v, w);
789  }else if(m_mode == "Dinv"){
790  mult_csw_inv(v, w);
791  }else if(m_mode == "H"){
792  H(v, w);
793  }else{
794  vout.crucial(m_vl, "%s: mode undefined.\n", class_name.c_str());
795  exit(EXIT_FAILURE);
796  }
797 
798 }
799 
800 //====================================================================
801 template<typename AFIELD>
802 void AFopr_CloverTerm<AFIELD>::mult(AFIELD &v, const AFIELD &w,
803  const std::string mode)
804 {
805  if(mode == "D"){
806  D(v, w);
807  }else if(mode == "H"){
808  H(v, w);
809  }else{
810  vout.crucial(m_vl, "%s: illegal mode is given to mult with mode\n",
811  class_name.c_str());
812  exit(EXIT_FAILURE);
813  }
814 
815 }
816 
817 //====================================================================
818 template<typename AFIELD>
819 void AFopr_CloverTerm<AFIELD>::mult_gm5(AFIELD &v, const AFIELD &w)
820 {
821  real_t* vp = v.ptr(0);
822  real_t* wp = const_cast<AFIELD*>(&w)->ptr(0);
823 
824  mult_gm5(vp, wp);
825 
826 }
827 
828 //====================================================================
829 template<typename AFIELD>
830 void AFopr_CloverTerm<AFIELD>::D(AFIELD &v, const AFIELD &w)
831 {
832  mult_csw(v, w);
833 }
834 
835 //====================================================================
836 template<typename AFIELD>
837 void AFopr_CloverTerm<AFIELD>::H(AFIELD &v, const AFIELD &w)
838 {
839  mult_csw(m_v2, w);
840  mult_gm5(v, m_v2);
841 }
842 
843 //====================================================================
844 template<typename AFIELD>
845 void AFopr_CloverTerm<AFIELD>::mult_gm5(real_t *v, real_t *w)
846 {
847 #pragma omp barrier
848 
849  int ith, nth;
850  set_thread(ith, nth);
851 
852  if(ith == 0){
853  if(m_repr == DIRAC){
854  BridgeACC::mult_wilson_gm5_dirac(v, w, m_Nsize, NC);
855  }else{
856  BridgeACC::mult_wilson_gm5_chiral(v, w, m_Nsize, NC);
857  }
858  }
859 
860 #pragma omp barrier
861 
862 }
863 
864 //====================================================================
865 template<typename AFIELD>
866 void AFopr_CloverTerm<AFIELD>::mult_csw(AFIELD& v, const AFIELD& w)
867 {
868 #pragma omp barrier
869 
870  int ith, nth;
871  set_thread(ith, nth);
872 
873  if(ith == 0){
874  real_t *v2 = v.ptr(0);
875  real_t *v1 = const_cast<AFIELD*>(&w)->ptr(0);
876  real_t *u = m_T.ptr(0);
877 
878  if(m_repr == DIRAC){
879  BridgeACC::mult_csw_dirac(v2, u, v1, m_Nsize, 0);
880  }else{
881  BridgeACC::mult_csw_chiral(v2, u, v1, m_Nsize, 0);
882  }
883  }
884 
885 #pragma omp barrier
886 
887 }
888 
889 //====================================================================
890 template<typename AFIELD>
891 void AFopr_CloverTerm<AFIELD>::mult_csw_inv(AFIELD& v, const AFIELD& w)
892 {
893  real_t *v2 = v.ptr(0);
894  real_t *v1 = const_cast<AFIELD*>(&w)->ptr(0);
895  real_t *u = m_Tinv.ptr(0);
896 
897 #pragma omp barrier
898 
899  int ith, nth;
900  set_thread(ith, nth);
901 
902  if(ith == 0){
903  if(m_repr == DIRAC){
904  BridgeACC::mult_csw_dirac(v2, u, v1, m_Nsize, 0);
905  }else{
906  BridgeACC::mult_csw_chiral(v2, u, v1, m_Nsize, 0);
907  }
908  }
909 
910 #pragma omp barrier
911 
912 }
913 
914 //====================================================================
915 template<typename AFIELD>
916 void AFopr_CloverTerm<AFIELD>::multadd_csw(AFIELD& v, const AFIELD& w)
917 {
918  real_t *v2 = v.ptr(0);
919  real_t *v1 = const_cast<AFIELD*>(&w)->ptr(0);
920  real_t *u = m_T.ptr(0);
921 
922 #pragma omp barrier
923 
924  int ith, nth;
925  set_thread(ith, nth);
926 
927  if(ith == 0){
928  if(m_repr == DIRAC){
929  BridgeACC::mult_csw_dirac(v2, u, v1, m_Nsize, 1);
930  }else{
931  BridgeACC::mult_csw_chiral(v2, u, v1, m_Nsize, 1);
932  }
933  }
934 
935 #pragma omp barrier
936 
937 }
938 
939 //====================================================================
940 template<typename AFIELD>
941 double AFopr_CloverTerm<AFIELD>::flop_count()
942 {
943  // The following counting explicitly depends on the implementation.
944  // It will be recalculated when the code is modified.
945  // The present counting is based on rev.1107. [24 Aug 2014 H.Matsufuru]
946 
947  int Lvol = CommonParameters::Lvol();
948  double flop_site, flop;
949 
950  if (m_repr == DIRAC) {
951  flop_site = static_cast<double>(
952  m_Nc * m_Nd * (4 + 6 * (4 * m_Nc + 2) + 2 * (4 * m_Nc + 1))
953  + 8 * m_Nc * m_Nc * m_Nd * m_Nd); // <- clover term
954  } else if (m_repr == CHIRAL) {
955  flop_site = static_cast<double>(
956  m_Nc * m_Nd * (4 + 8 * (4 * m_Nc + 2))
957  + 8 * m_Nc * m_Nc * m_Nd * m_Nd); // <- clover term
958  } else {
959  // vout.crucial(m_vl, "%s: input repr is undefined.\n",
960  // class_name.c_str());
961  vout.crucial(m_vl, "%s: input repr is undefined.\n");
962  abort();
963  }
964 
965  flop = flop_site * static_cast<double>(Lvol);
966  if ((m_mode == "DdagD") || (m_mode == "DDdag")) flop *= 2.0;
967 
968  return flop;
969 }
970 
971 //============================================================END=====
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
BridgeACC::mult_Gnd
void mult_Gnd(double *restrict u, const int exu, double *restrict v, const int exv, double *restrict w, const int exw, const int nst)
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
AFopr_CloverTerm::set_csw_chiral
void set_csw_chiral(Field &u)
setting clover term (Chiral repr).
Definition: afopr_CloverTerm-tmpl.h:612
CommonParameters::Lvol
static long_t Lvol()
Definition: commonParameters.h:95
convert_gauge
void convert_gauge(INDEX &index, AFIELD &v, const Field &w)
Definition: afield-inc.h:223
AFopr_CloverTerm::solve_csw_inv_chiral
void solve_csw_inv_chiral()
Definition: afopr_CloverTerm-tmpl.h:379
BridgeACC::ah_G
void ah_G(double *u, const int ex, const int nst)
Definition: afield_Gauge_openacc-inc.h:643
BridgeACC::mult_wilson_gm5_chiral
void mult_wilson_gm5_chiral(double *RESTRICT v2, double *RESTRICT v1, int *Nsize, int Nc)
Definition: mult_Wilson_openacc-inc.h:49
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
AFopr_CloverTerm::multadd_csw
void multadd_csw(AFIELD &, const AFIELD &)
Definition: afopr_CloverTerm-tmpl.h:916
AFopr_CloverTerm::get_csw_inv
void get_csw_inv(AFIELD &T)
getting the inverse of clover term.
Definition: afopr_CloverTerm-tmpl.h:465
AFopr_CloverTerm::set_parameters
void set_parameters(const Parameters &params)
setting parameters by a Parameter object.
Definition: afopr_CloverTerm-tmpl.h:105
ThreadManager::get_num_threads
static int get_num_threads()
returns available number of threads.
Definition: threadManager.cpp:246
AFopr_CloverTerm::tidyup
void tidyup()
final tidy-up.
Definition: afopr_CloverTerm-tmpl.h:94
AFopr_CloverTerm::mult
void mult(AFIELD &, const AFIELD &)
multiplies fermion operator to a given field.
Definition: afopr_CloverTerm-tmpl.h:768
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
AFopr_CloverTerm::set_config_omp
void set_config_omp(Field *u)
Definition: afopr_CloverTerm-tmpl.h:199
Parameters
Class for parameters.
Definition: parameters.h:46
AIndex_lex
Definition: aindex_lex_base.h:17
AFopr_CloverTerm::set_mode
void set_mode(std::string mode)
setting mult mode.
Definition: afopr_CloverTerm-tmpl.h:748
update_device
void update_device(AField< REALTYPE, ACCEL > &v)
Definition: afield-inc.h:33
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
ASolver_CG
Definition: asolver_CG.h:16
AFopr_CloverTerm::solve_csw_inv_dirac
void solve_csw_inv_dirac()
Definition: afopr_CloverTerm-tmpl.h:321
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Bridge::BridgeIO::detailed
void detailed(const char *format,...)
Definition: bridgeIO.cpp:281
AFopr_CloverTerm
Definition: afopr_CloverTerm.h:31
BridgeACC::mult_wilson_gm5_dirac
void mult_wilson_gm5_dirac(double *RESTRICT v2, double *RESTRICT v1, int *Nsize, int Nc)
Definition: mult_Wilson_openacc-inc.h:22
Field::check_size
bool check_size(const int nin, const int nvol, const int nex) const
checking size parameters. [23 May 2016 H.Matsufuru]
Definition: field.h:135
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
AFopr_CloverTerm::H
void H(AFIELD &, const AFIELD &)
Definition: afopr_CloverTerm-tmpl.h:837
update_host
void update_host(AField< REALTYPE, ACCEL > &v)
Definition: afield-inc.h:26
AFopr_CloverTerm::set_fieldstrength
void set_fieldstrength(AFIELD &Fst, AFIELD &u, int, int)
setting field strength.
Definition: afopr_CloverTerm-tmpl.h:719
Timer
Definition: timer.h:31
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:300
AFopr_CloverTerm::set_csw
void set_csw(Field &u)
setting clover term.
Definition: afopr_CloverTerm-tmpl.h:255
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
AFopr_CloverTerm::flop_count
double flop_count()
returns floating operation counts.
Definition: afopr_CloverTerm-tmpl.h:941
Timer::start
void start()
Definition: timer.cpp:44
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
BridgeACC::multadd_Gdn
void multadd_Gdn(double *restrict u, const int exu, double *restrict v, const int exv, double *restrict w, const int exw, const double a, const int nst)
AFopr_CloverTerm::set_csw_dirac
void set_csw_dirac(Field &u)
setting clover term (Dirac repr).
Definition: afopr_CloverTerm-tmpl.h:517
AStaple_lex
Staple construction.
Definition: afopr_CloverTerm.h:27
AFopr_CloverTerm::set_config
void set_config(Field *u)
setting gauge configuration.
Definition: afopr_CloverTerm-tmpl.h:186
NC
#define NC
Definition: field_F_imp_SU2-inc.h:15
AFopr_CloverTerm::mult_gm5
void mult_gm5(AFIELD &, const AFIELD &)
multiplies gamma_5 matrix.
Definition: afopr_CloverTerm-tmpl.h:819
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
AFopr_CloverTerm::get_csw
void get_csw(AFIELD &T)
getting clover term.
Definition: afopr_CloverTerm-tmpl.h:451
ND
#define ND
Definition: field_F_imp_SU2-inc.h:18
SU_N::xI
Mat_SU_N xI(const Mat_SU_N &u)
Definition: mat_SU_N.h:586
real_t
double real_t
Definition: bridgeACC_AField_double.cpp:14
BridgeACC::mult_Gdn
void mult_Gdn(double *restrict u, const int exu, double *restrict v, const int exv, double *restrict w, const int exw, const int nst)
AFopr_CloverTerm::mult_csw
void mult_csw(AFIELD &, const AFIELD &)
Definition: afopr_CloverTerm-tmpl.h:866
AFopr_CloverTerm< Field >::real_t
Field ::real_t real_t
Definition: afopr_CloverTerm.h:34
Field::ptr
const double * ptr(const int jin, const int site, const int jex) const
Definition: field.h:153
CommonParameters::Nd
static int Nd()
Definition: commonParameters.h:116
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
BridgeACC::mult_csw_chiral
void mult_csw_chiral(double *RESTRICT v2, double *RESTRICT u, double *RESTRICT v1, int *Nsize, int flag)
Definition: mult_CloverTerm_openacc-inc.h:110
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
BridgeACC::mult_csw_dirac
void mult_csw_dirac(double *RESTRICT v2, double *RESTRICT u, double *RESTRICT v1, int *Nsize, int flag)
Definition: mult_CloverTerm_openacc-inc.h:15
AFopr_CloverTerm::mult_dag
void mult_dag(AFIELD &, const AFIELD &)
hermitian conjugate of mult.
Definition: afopr_CloverTerm-tmpl.h:785
AFopr_CloverTerm::solve_csw_inv
void solve_csw_inv()
solve the inverse of clover term.
Definition: afopr_CloverTerm-tmpl.h:271
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
BridgeACC::add_unit
void add_unit(double *u, const int ex, const double a, const int nst)
Definition: afield_Gauge_openacc-inc.h:767
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
Parameters::fetch_double
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
AFopr_CloverTerm::set_config_impl
void set_config_impl(Field *u)
Definition: afopr_CloverTerm-tmpl.h:209
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
AFopr_CloverTerm::get_mode
std::string get_mode() const
returns mult mode.
Definition: afopr_CloverTerm-tmpl.h:761
Field
Container of Field-type object.
Definition: field.h:46
Timer::elapsed_sec
double elapsed_sec() const
Definition: timer.cpp:107
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
AFopr_CloverTerm::init
void init(const Parameters &params)
initial setup.
Definition: afopr_CloverTerm-tmpl.h:16
AFopr_CloverTerm::D
void D(AFIELD &, const AFIELD &)
Definition: afopr_CloverTerm-tmpl.h:830
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
AFopr_CloverTerm::mult_csw_inv
void mult_csw_inv(AFIELD &, const AFIELD &)
Definition: afopr_CloverTerm-tmpl.h:891
Timer::stop
void stop()
Definition: timer.cpp:69
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
BridgeACC::multadd_Gnd
void multadd_Gnd(double *restrict u, const int exu, double *restrict v, const int exv, double *restrict w, const int exw, const double a, const int nst)
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:572
Timer::reset
void reset()
Definition: timer.cpp:97