Bridge++  Ver.2.1.3
staple_SF.cpp
Go to the documentation of this file.
1 
10 #include "lib/Measurements/Gauge/staple_SF.h"
11 #include "lib/ResourceManager/threadManager.h"
12 #include "lib/Field/field_thread-inc.h"
13 
30 const std::string Staple_SF::class_name = "Staple_SF";
31 
32 //====================================================================
33 void Staple_SF::init(const Parameters& params)
34 {
35  ThreadManager::assert_single_thread(class_name);
36 
37  std::string vlevel;
38  if (!params.fetch_string("verbose_level", vlevel)) {
39  m_vl = vout.set_verbose_level(vlevel);
40  } else {
41  m_vl = CommonParameters::Vlevel();
42  }
43 
44  vout.general(m_vl, "%s: construction\n", class_name.c_str());
45  vout.increase_indent();
46 
47  m_initialized = 0;
48 
49  set_parameters(params);
50 
51  const int Nvol = CommonParameters::Nvol();
52 
53  m_Umu = new Field_G(Nvol, 1);
54  m_Unu = new Field_G(Nvol, 1);
55  m_vt = new Field_G(Nvol, 1);
56  m_wt = new Field_G(Nvol, 1);
57 
58  m_stpl_upper = new Field_G(Nvol, 1);
59  m_stpl_lower = new Field_G(Nvol, 1);
60 
61 
62  vout.decrease_indent();
63  vout.general(m_vl, "%s: construction finished.\n",
64  class_name.c_str());
65 
66 }
67 
68 //====================================================================
69 void Staple_SF::init()
70 {
71  ThreadManager::assert_single_thread(class_name);
72 
73  m_vl = CommonParameters::Vlevel();
74 
75  vout.general(m_vl, "%s (obsolete): construction\n",
76  class_name.c_str());
77  vout.increase_indent();
78 
79 
80  m_initialized = 0;
81 
82  const int Nvol = CommonParameters::Nvol();
83 
84  m_Umu = new Field_G(Nvol, 1);
85  m_Unu = new Field_G(Nvol, 1);
86  m_vt = new Field_G(Nvol, 1);
87  m_wt = new Field_G(Nvol, 1);
88 
89  m_stpl_upper = new Field_G(Nvol, 1);
90  m_stpl_lower = new Field_G(Nvol, 1);
91 
92  vout.decrease_indent();
93  vout.general(m_vl, "%s: construction finished.\n",
94  class_name.c_str());
95 
96 }
97 
98 
99 //====================================================================
100 void Staple_SF::tidyup()
101 {
102  delete m_Umu;
103  delete m_Unu;
104  delete m_vt;
105  delete m_wt;
106 
107  delete m_stpl_upper;
108  delete m_stpl_lower;
109 }
110 
111 //====================================================================
112 void Staple_SF::set_parameters(const Parameters& params)
113 {
114  std::string vlevel;
115  if (!params.fetch_string("verbose_level", vlevel)) {
116  m_vl = vout.set_verbose_level(vlevel);
117  }
118 
119  //- fetch and check input parameters
120  std::vector<double> phi, phipr, p_omega;
121 
122  int err = 0;
123  err += params.fetch_double_vector("phi", phi);
124  err += params.fetch_double_vector("phipr", phipr);
125  //err += params.fetch_double_vector("p_omega", p_omega);
126 
127  if (err) {
128  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
129  class_name.c_str());
130  exit(EXIT_FAILURE);
131  }
132 
133  if (params.fetch_double_vector("p_omega", p_omega)) {
134  vout.general("%s: p_omega is not given: set to default values.\n",
135  class_name.c_str());
136  p_omega.resize(3);
137  p_omega[0] = 1.0;
138  p_omega[1] = -0.5;
139  p_omega[2] = -0.5;
140  }
141 
142  set_parameters(phi, phipr, p_omega); // call std::vector version
143 }
144 
145 
146 //====================================================================
147 void Staple_SF::get_parameters(Parameters& params) const
148 {
149  params.set_double_vector("phi", m_phi);
150  params.set_double_vector("phipr", m_phipr);
151  params.set_double_vector("p_omega", m_p_omega);
152 
153  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
154 }
155 
156 
157 //====================================================================
158 
174 void Staple_SF::set_parameters(const std::vector<double>& phi,
175  const std::vector<double>& phipr,
176  const std::vector<double>& p_omega)
177 {
178 #pragma omp barrier
179 
180  assert(phi.size() == 3);
181  assert(phipr.size() == 3);
182  assert(p_omega.size() == 3);
183 
184  int ith = ThreadManager::get_thread_id();
185  if (ith == 0) {
186 
187  m_i_omega0.zero();
188  m_i_omega0.set(0, 0, 0.0, p_omega[0]);
189  m_i_omega0.set(1, 1, 0.0, p_omega[1]);
190  m_i_omega0.set(2, 2, 0.0, p_omega[2]);
191 
192  m_initialized = 1;
193 
194  m_phi.resize(3);
195  m_phi[0] = phi[0];
196  m_phi[1] = phi[1];
197  m_phi[2] = phi[2];
198 
199  m_phipr.resize(3);
200  m_phipr[0] = phipr[0];
201  m_phipr[1] = phipr[1];
202  m_phipr[2] = phipr[2];
203 
204  m_p_omega.resize(3);
205  m_p_omega[0] = p_omega[0];
206  m_p_omega[1] = p_omega[1];
207  m_p_omega[2] = p_omega[2];
208 
209  Field_SF::set_boundary_matrix(m_wk, m_phi);
210  Field_SF::set_boundary_matrix(m_wkpr, m_phipr);
211  }
212 #pragma omp barrier
213 
214  //- print input parameters
215  vout.general(m_vl, "%s: parameters\n", class_name.c_str());
216  vout.general(m_vl, " phi1 = %12.6f\n", m_phi[0]);
217  vout.general(m_vl, " phi2 = %12.6f\n", m_phi[1]);
218  vout.general(m_vl, " phi3 = %12.6f\n", m_phi[2]);
219  vout.general(m_vl, " phipr1 = %12.6f\n", m_phipr[0]);
220  vout.general(m_vl, " phipr2 = %12.6f\n", m_phipr[1]);
221  vout.general(m_vl, " phipr3 = %12.6f\n", m_phipr[2]);
222  vout.general(m_vl, " p_omega1 = %12.6f\n", m_p_omega[0]);
223  vout.general(m_vl, " p_omega2 = %12.6f\n", m_p_omega[1]);
224  vout.general(m_vl, " p_omega3 = %12.6f\n", m_p_omega[2]);
225 
226 }
227 
228 //====================================================================
229 void Staple_SF::set_parameters(const std::vector<double>& phi,
230  const std::vector<double>& phipr)
231 {
232 #pragma omp barrier
233 
234  int ith = ThreadManager::get_thread_id();
235  if (ith == 0) {
236 
237  m_i_omega0.zero();
238  m_i_omega0.set(0, 0, 0.0, 1.0);
239  m_i_omega0.set(1, 1, 0.0, -0.5);
240  m_i_omega0.set(2, 2, 0.0, -0.5);
241 
242  m_initialized = 1;
243 
244  m_phi.resize(3);
245  m_phi[0] = phi[0];
246  m_phi[1] = phi[1];
247  m_phi[2] = phi[2];
248 
249  m_phipr.resize(3);
250  m_phipr[0] = phipr[0];
251  m_phipr[1] = phipr[1];
252  m_phipr[2] = phipr[2];
253 
254  Field_SF::set_boundary_matrix(m_wk, m_phi);
255  Field_SF::set_boundary_matrix(m_wkpr, m_phipr);
256  }
257 #pragma omp barrier
258 
259 }
260 
261 //====================================================================
262 
306 double Staple_SF::sf_coupling_plaq(const Field_G& U, const double ct)
307 {
308 #pragma omp barrier
309 
310  if (!m_initialized) {
311  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
312  class_name.c_str());
313  exit(EXIT_FAILURE);
314  }
315 
316  const int Nc = CommonParameters::Nc();
317  const int Ndim = CommonParameters::Ndim();
318 
319  const int Nx = CommonParameters::Nx();
320  const int Ny = CommonParameters::Ny();
321  const int Nz = CommonParameters::Nz();
322  const int Nt = CommonParameters::Nt();
323  const int NPEt = CommonParameters::NPEt();
324 
325  const int Lx = CommonParameters::Lx();
326 
327  double plaq = 0.0;
328  double plaqt0 = 0.0;
329  double plaqtT = 0.0;
330 
331  const int Nxyz = Nx * Ny * Nz;
332 
333  int ith, nth, is, ns;
334  set_threadtask(ith, nth, is, ns, Nxyz);
335 
336  for (int nu = 0; nu < Ndim - 1; nu++) {
337  upper(*m_stpl_upper, U, 3, nu);
338 
339  for (int ixyz = is; ixyz < ns; ixyz++) {
340  int x = ixyz % Nx;
341  int y = (ixyz/Nx) % Ny;
342  int z = ixyz/(Nx * Ny);
343 
344  // boundary
345  if (Communicator::ipe(3) == 0) {
346  int t = 0;
347  int site = m_index.site(x, y, z, t);
348 
349  Mat_SU_N up(Nc);
350  up = m_stpl_upper->mat(site) * U.mat_dag(site, 3);
351  double scr = ReTr(m_i_omega0 * up);
352 
353  /*
354  up.unit();
355  up *= staple.mat(site);
356  up *= U->mat_dag(site,3);
357  up *= m_i_omega0;
358  scr = ReTr( up );
359  */
360  plaq -= scr;
361  plaqt0 += scr;
362  }
363 #pragma omp barrier
364 
365  // boundary
366  if (Communicator::ipe(3) == NPEt - 1) {
367  int t = Nt - 1;
368  int site = m_index.site(x, y, z, t);
369 
370  Mat_SU_N up(Nc);
371  up = m_stpl_upper->mat_dag(site) * U.mat(site, 3);
372  double scr = ReTr(m_i_omega0 * up);
373 
374  /*
375  up.unit();
376  up *= staple.mat_dag(site);
377  up *= U->mat(site,3);
378  up *= m_i_omega0;
379  scr = ReTr( up );
380  */
381  plaq += scr;
382  plaqtT += scr;
383  }
384 #pragma omp barrier
385  }
386  }
387 #pragma omp barrier
388 
389  // plaq = Communicator::reduce_sum(plaq);
390  // plaqt0 = Communicator::reduce_sum(plaqt0);
391  // plaqtT = Communicator::reduce_sum(plaqtT);
392  ThreadManager::reduce_sum_global(plaq, ith, nth);
393  ThreadManager::reduce_sum_global(plaqt0, ith, nth);
394  ThreadManager::reduce_sum_global(plaqtT, ith, nth);
395 
396  plaq *= ct / Lx;
397  plaqt0 *= ct / Lx;
398  plaqtT *= ct / Lx;
399 
400  vout.general(m_vl, "SF_delSg_plaq, from 0, from T = %.8f %.8f %.8f\n",
401  plaq, plaqt0, plaqtT);
402 
403  return plaq;
404 }
405 
406 
407 //====================================================================
408 
507 double Staple_SF::sf_coupling_rect(const Field_G& U, const double ctr)
508 {
509 #pragma omp barrier
510 
511  if (!m_initialized) {
512  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
513  class_name.c_str());
514  exit(EXIT_FAILURE);
515  }
516 
517  const int Nc = CommonParameters::Nc();
518  const int Ndim = CommonParameters::Ndim();
519  const int Nvol = CommonParameters::Nvol();
520 
521  const int Nx = CommonParameters::Nx();
522  const int Ny = CommonParameters::Ny();
523  const int Nz = CommonParameters::Nz();
524  const int Nt = CommonParameters::Nt();
525  const int NPEt = CommonParameters::NPEt();
526 
527  const int Lx = CommonParameters::Lx();
528 
529  const int Nxyz = Nx * Ny * Nz;
530 
531  const int nu = 3;
532 
533  double rect01 = 0.0;
534  double rect02 = 0.0;
535  double rect03 = 0.0;
536  double rectt1 = 0.0;
537  double rectt2 = 0.0;
538  double rectt3 = 0.0;
539 
540  int ith, nth, is, ns;
541  set_threadtask(ith, nth, is, ns, Nxyz);
542 
543 
544  for (int mu = 0; mu < Ndim - 1; mu++) {
545  // rect01
546  // <---<---+
547  // | |
548  // t=0 x--->---+
549  // omega0
550 
551  // rect02
552  // <---<---+
553  // | |
554  // t=0 +---x---+
555  // omega0
556 
557  // rectt1
558  // omega0
559  // t=Nt x--->---+
560  // | |
561  // +---<---+
562 
563  // rectt2
564  // omega0
565  // t=Nt +---x---+
566  // | |
567  // <---<---+
568 
569  upper(*m_stpl_upper, U, nu, mu);
570 
571  copy(*m_Umu, 0, U, mu);
572 
573  copy(*m_Unu, 0, U, nu);
574 
575  m_shift.backward(*m_vt, *m_stpl_upper, mu);
576 
577  m_shift.backward(*m_wt, *m_Umu, nu);
578 
579  for (int ixyz = is; ixyz < ns; ixyz++) {
580  int x = ixyz % Nx;
581  int y = (ixyz/Nx) % Ny;
582  int z = ixyz/(Nx * Ny);
583  int t = 0;
584 
585  if (Communicator::ipe(3) == 0) {
586  int site = m_index.site(x, y, z, t);
587 
588  Mat_SU_N wmat(Nc);
589  wmat = m_wk * m_vt->mat(site);
590 
591  Mat_SU_N cmat(Nc);
592  cmat = wmat * m_wt->mat_dag(site);
593 
594  wmat = cmat * m_Unu->mat_dag(site);
595  rect01 += ReTr(m_i_omega0 * wmat);
596 
597  wmat = m_i_omega0 * m_vt->mat(site);
598  cmat = wmat * m_wt->mat_dag(site);
599  wmat = cmat * m_Unu->mat_dag(site);
600  rect02 += ReTr(m_wk * wmat);
601  }
602 
603  t = Nt - 1;
604  if (Communicator::ipe(3) == NPEt - 1) {
605  int site = m_index.site(x, y, z, t);
606 
607  Mat_SU_N cmat(Nc);
608  cmat = m_i_omega0 * m_wkpr;
609 
610  Mat_SU_N wmat(Nc);
611  wmat = cmat * m_vt->mat_dag(site);
612 
613  cmat = m_Unu->mat(site) * wmat;
614  rectt1 += ReTr(cmat * U.mat_dag(site, mu));
615 
616  cmat = m_i_omega0 * m_vt->mat_dag(site);
617  wmat = m_wkpr * cmat;
618  cmat = m_Unu->mat(site) * wmat;
619  rectt2 += ReTr(cmat * U.mat_dag(site, mu));
620  }
621  }
622 #pragma omp barrier
623 
624  // rect03
625  // +---+
626  // | |
627  // v ^
628  // | |
629  // t=0 x--->
630  // omega0
631 
632  upper(*m_stpl_upper, U, mu, nu);
633 
634  m_shift.backward(*m_vt, *m_Unu, mu);
635  m_shift.backward(*m_wt, *m_stpl_upper, nu);
636 
637  for (int ixyz = is; ixyz < ns; ixyz++) {
638  int x = ixyz % Nx;
639  int y = (ixyz/Nx) % Ny;
640  int z = ixyz/(Nx * Ny);
641  {
642  {
643  int t = 0;
644  if (Communicator::ipe(3) == 0) {
645  int site = m_index.site(x, y, z, t);
646 
647  Mat_SU_N wmat(Nc);
648  wmat = m_wt->mat(site) * m_vt->mat_dag(site);
649 
650  Mat_SU_N cmat(Nc);
651  cmat = m_Unu->mat(site) * wmat;
652 
653  wmat = m_wk * cmat.dag();
654  rect03 += ReTr(m_i_omega0 * wmat);
655  }
656  }
657  }
658  }
659 #pragma omp barrier
660 
661  // rectt3
662  // omega0
663  // t=Nt x--->
664  // | |
665  // ^ v
666  // | |
667  // +---+
668 
669  lower(*m_stpl_lower, U, mu, nu);
670 
671  m_shift.backward(*m_vt, *m_Unu, mu);
672 
673  for (int ixyz = is; ixyz < ns; ixyz++) {
674  int x = ixyz % Nx;
675  int y = (ixyz/Nx) % Ny;
676  int z = ixyz/(Nx * Ny);
677  int t = Nt - 1;
678 
679  if (Communicator::ipe(3) == NPEt - 1) {
680  int site = m_index.site(x, y, z, t);
681 
682  Mat_SU_N wmat(Nc);
683  wmat = m_i_omega0 * m_wkpr;
684 
685  Mat_SU_N cmat(Nc);
686  cmat = wmat * m_vt->mat_dag(site);
687 
688  wmat = cmat * m_stpl_lower->mat_dag(site);
689  rectt3 += ReTr(m_Unu->mat(site) * wmat);
690  }
691  }
692 #pragma omp barrier
693  }
694  // rect01 = Communicator::reduce_sum(rect01);
695  // rect02 = Communicator::reduce_sum(rect02);
696  // rect03 = Communicator::reduce_sum(rect03);
697  // rectt1 = Communicator::reduce_sum(rectt1);
698  // rectt2 = Communicator::reduce_sum(rectt2);
699  // rectt3 = Communicator::reduce_sum(rectt3);
700  ThreadManager::reduce_sum_global(rect01, ith, nth);
701  ThreadManager::reduce_sum_global(rect02, ith, nth);
702  ThreadManager::reduce_sum_global(rect03, ith, nth);
703  ThreadManager::reduce_sum_global(rectt1, ith, nth);
704  ThreadManager::reduce_sum_global(rectt2, ith, nth);
705  ThreadManager::reduce_sum_global(rectt3, ith, nth);
706 
707  rect01 *= ctr / Lx;
708  rect02 *= ctr / Lx;
709  rect03 /= Lx;
710  rectt1 *= ctr / Lx;
711  rectt2 *= ctr / Lx;
712  rectt3 /= Lx;
713 
714  double rect = -rect01 - rect02 - rect03 + rectt1 + rectt2 + rectt3;
715 
716  vout.general(m_vl, "SF_delSg_rect, at 01, 02, 03, at T1, T2, T3 = %.8f %.8f %.8f %.8f %.8f %.8f %.8f\n",
717  rect, rect01, rect02, rect03, rectt1, rectt2, rectt3);
718 
719  return rect;
720 }
721 
722 
723 //====================================================================
724 
733 double Staple_SF::plaquette(const Field_G& U)
734 {
735  if (!m_initialized) {
736  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
737  class_name.c_str());
738  exit(EXIT_FAILURE);
739  }
740 
741  return(plaq_s(U) + plaq_t(U));
742 }
743 
744 
745 //====================================================================
746 
761 double Staple_SF::plaquette_ct(const Field_G& U, const double ct)
762 {
763  if (!m_initialized) {
764  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
765  class_name.c_str());
766  exit(EXIT_FAILURE);
767  }
768 
769  return(plaq_s(U) + plaq_t_ct(U, ct));
770 }
771 
772 
773 //====================================================================
774 
783 double Staple_SF::plaq_s(const Field_G& U)
784 {
785 #pragma omp barrier
786 
787  const int Nvol = CommonParameters::Nvol();
788 
789  int ith, nth, is, ns;
790  set_threadtask(ith, nth, is, ns, Nvol);
791 
792  double plaq = 0.0;
793 
794  upper(*m_stpl_upper, U, 0, 1);
795  for (int site = is; site < ns; site++) {
796  plaq += ReTr(U.mat(site, 0) * m_stpl_upper->mat_dag(site)); // P_xy
797  }
798 
799  upper(*m_stpl_upper, U, 1, 2);
800  for (int site = is; site < ns; site++) {
801  plaq += ReTr(U.mat(site, 1) * m_stpl_upper->mat_dag(site)); // P_yz
802  }
803 
804  upper(*m_stpl_upper, U, 2, 0);
805  for (int site = is; site < ns; site++) {
806  plaq += ReTr(U.mat(site, 2) * m_stpl_upper->mat_dag(site)); // P_zx
807  }
808 
809  ThreadManager::reduce_sum_global(plaq, ith, nth);
810 
811  return plaq;
812 }
813 
814 //qqqqq
815 //====================================================================
816 
825 double Staple_SF::plaq_t(const Field_G& U)
826 {
827 #pragma omp barrier
828 
829  if (!m_initialized) {
830  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
831  class_name.c_str());
832  exit(EXIT_FAILURE);
833  }
834 
835  const int Ndim = CommonParameters::Ndim();
836  const int Nvol = CommonParameters::Nvol();
837 
838  int ith, nth, is, ns;
839  set_threadtask(ith, nth, is, ns, Nvol);
840 
841  double plaq = 0.0;
842 
843  for (int nu = 0; nu < Ndim - 1; nu++) {
844  lower(*m_stpl_lower, U, 3, nu);
845  // staple = upper(U,3,nu);
846  for (int site = is; site < ns; site++) {
847  plaq += ReTr(U.mat(site, 3) * m_stpl_lower->mat_dag(site)); // P_tk
848  }
849 #pragma omp barrier
850  }
851 
852  //plaq = Communicator::reduce_sum(plaq);
853  ThreadManager::reduce_sum_global(plaq, ith, nth);
854 
855  return plaq;
856 }
857 
858 
859 //====================================================================
860 
875 double Staple_SF::plaq_t_ct(const Field_G& U, const double ct)
876 {
877 #pragma omp barrier
878 
879  if (!m_initialized) {
880  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
881  class_name.c_str());
882  exit(EXIT_FAILURE);
883  }
884 
885  const int Ndim = CommonParameters::Ndim();
886  const int Nt = CommonParameters::Nt();
887  const int Nvol = CommonParameters::Nvol();
888  const int NPEt = CommonParameters::NPEt();
889 
890  //Field_G staple;
891  double plaq = 0.0;
892 
893  int ith, nth, is, ns;
894  set_threadtask(ith, nth, is, ns, Nvol);
895 
896  for (int nu = 0; nu < Ndim - 1; nu++) {
897  lower(*m_stpl_lower, U, 3, nu);
898  // If the node is at the boundary the temporal plaquette is multiplied with ct.
899  if (Communicator::ipe(3) == 0) {
900  Field_SF::mult_ct_boundary(*m_stpl_lower, 0, ct);
901  }
902  if (Communicator::ipe(3) == NPEt - 1) {
903  Field_SF::mult_ct_boundary(*m_stpl_lower, Nt - 1, ct);
904  }
905 #pragma omp barrier
906 
907  for (int site = is; site < ns; site++) {
908  plaq += ReTr(U.mat(site, 3) * m_stpl_lower->mat_dag(site)); // P_tk
909  }
910 #pragma omp barrier
911  }
912 
913  // plaq = Communicator::reduce_sum(plaq);
914  ThreadManager::reduce_sum_global(plaq, ith, nth);
915 
916  return plaq;
917 }
918 
919 
920 //====================================================================
921 
940 void Staple_SF::staple(Field_G& W, const Field_G& U, const int mu)
941 {
942 #pragma omp barrier
943 
944  if (!m_initialized) {
945  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
946  class_name.c_str());
947  exit(EXIT_FAILURE);
948  }
949 
950  const int Ndim = CommonParameters::Ndim();
951 
952  W.set(0.0);
953 #pragma omp barrier
954 
955  for (int nu = 0; nu < Ndim; nu++) {
956  if (nu != mu) {
957  // Field_G c_tmp;
958 
959  upper(*m_stpl_upper, U, mu, nu);
960  axpy(W, 1.0, *m_stpl_upper);
961 
962  lower(*m_stpl_lower, U, mu, nu);
963  axpy(W, 1.0, *m_stpl_lower);
964 #pragma omp barrier
965  }
966  }
967 }
968 
969 
970 //====================================================================
971 
990 void Staple_SF::staple_ct(Field_G& W, const Field_G& U, const int mu, const double ct)
991 {
992 #pragma omp barrier
993 
994  const int Ndim = CommonParameters::Ndim();
995  const int Nt = CommonParameters::Nt();
996  const int NPEt = CommonParameters::NPEt();
997 
998  W.set(0.0);
999 #pragma omp barrier
1000 
1001  for (int nu = 0; nu < Ndim; nu++) {
1002  if (nu != mu) {
1003  //Field_G staple_upper;
1004  //Field_G staple_lower;
1005 
1006  upper(*m_stpl_upper, U, mu, nu);
1007  lower(*m_stpl_lower, U, mu, nu);
1008 
1009  if (Communicator::ipe(3) == 0) {
1010  if (mu == 3) {
1011  Field_SF::mult_ct_boundary(*m_stpl_upper, 0, ct);
1012  Field_SF::mult_ct_boundary(*m_stpl_lower, 0, ct);
1013  }
1014  if (nu == 3) {
1015  Field_SF::mult_ct_boundary(*m_stpl_lower, 1, ct);
1016  }
1017  }
1018 
1019  if (Communicator::ipe(3) == NPEt - 1) {
1020  if (mu == 3) {
1021  Field_SF::mult_ct_boundary(*m_stpl_upper, Nt - 1, ct);
1022  Field_SF::mult_ct_boundary(*m_stpl_lower, Nt - 1, ct);
1023  }
1024  if (nu == 3) {
1025  Field_SF::mult_ct_boundary(*m_stpl_upper, Nt - 1, ct);
1026  }
1027  }
1028 #pragma omp barrier
1029 
1030  axpy(W, 1.0, *m_stpl_upper);
1031  axpy(W, 1.0, *m_stpl_lower);
1032 #pragma omp barrier
1033  }
1034  }
1035 }
1036 
1037 
1038 //====================================================================
1039 void Staple_SF::upper(Field_G& c, const Field_G& U, const int mu, const int nu)
1040 {
1041 #pragma omp barrier
1042 
1043  if (!m_initialized) {
1044  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
1045  class_name.c_str());
1046  exit(EXIT_FAILURE);
1047  }
1048 
1049  const int Nvol = CommonParameters::Nvol();
1050 
1051  // (1) mu (2)
1052  // +-->--+
1053  // nu | |
1054  // i+ +
1055 
1056  copy(*m_Umu, 0, U, mu);
1057  copy(*m_Unu, 0, U, nu);
1058 #pragma omp barrier
1059 
1060  if (mu != 3) Field_SF::set_boundary_wk(*m_Umu, m_wk);
1061  if (nu != 3) Field_SF::set_boundary_wk(*m_Unu, m_wk);
1062 #pragma omp barrier
1063 
1064  m_shift.backward(*m_vt, *m_Unu, mu);
1065  m_shift.backward(c, *m_Umu, nu);
1066 
1067  if (mu == 3) Field_SF::set_boundary_wkpr(*m_vt, m_wkpr);
1068  if (nu == 3) Field_SF::set_boundary_wkpr(c, m_wkpr);
1069 #pragma omp barrier
1070 
1071  mult_Field_Gnd(*m_wt, 0, c, 0, *m_vt, 0);
1072  mult_Field_Gnn(c, 0, *m_Unu, 0, *m_wt, 0);
1073 #pragma omp barrier
1074 
1075  if (mu != 3) Field_SF::set_boundary_zero(c);
1076 #pragma omp barrier
1077 }
1078 
1079 
1080 //====================================================================
1081 void Staple_SF::lower(Field_G& c, const Field_G& U, const int mu, const int nu)
1082 {
1083 #pragma omp barrier
1084 
1085  if (!m_initialized) {
1086  vout.crucial(m_vl, "Error at %s: Parameter is not initialized.\n",
1087  class_name.c_str());
1088  exit(EXIT_FAILURE);
1089  }
1090 
1091  const int Nvol = CommonParameters::Nvol();
1092 
1093  // + +
1094  // nu | |
1095  // i+-->--+
1096  // (1) mu (2)
1097 
1098  copy(*m_Umu, 0, U, mu);
1099  copy(*m_Unu, 0, U, nu);
1100 #pragma omp barrier
1101 
1102  if (mu != 3) Field_SF::set_boundary_wk(*m_Umu, m_wk);
1103  if (nu != 3) Field_SF::set_boundary_wk(*m_Unu, m_wk);
1104 #pragma omp barrier
1105 
1106  m_shift.backward(*m_wt, *m_Unu, mu);
1107  if (mu == 3) Field_SF::set_boundary_wkpr(*m_wt, m_wkpr);
1108 #pragma omp barrier
1109 
1110  mult_Field_Gnn(*m_vt, 0, *m_Umu, 0, *m_wt, 0);
1111  mult_Field_Gdn(*m_wt, 0, *m_Unu, 0, *m_vt, 0);
1112 #pragma omp barrier
1113 
1114  m_shift.forward(c, *m_wt, nu);
1115 
1116  if (mu != 3) Field_SF::set_boundary_zero(c);
1117 #pragma omp barrier
1118 }
1119 
1120 
1121 //====================================================================
1122 
1140 void Staple_SF::print_plaquette(const Field_G& U)
1141 {
1142 #pragma omp barrier
1143 
1144  const int Lx = CommonParameters::Lx();
1145  const int Ly = CommonParameters::Ly();
1146  const int Lz = CommonParameters::Lz();
1147  const int Lt = CommonParameters::Lt();
1148 
1149  const double plaq = plaquette(U);
1150  const double plaq2 = plaq + 3 * 3 * Lx * Ly * Lz;
1151 
1152  vout.general(m_vl, "plaq_SF without boundary spatial plaq = %.8f\n",
1153  plaq / (3 * Lx * Ly * Lz * (6 * Lt - 3)));
1154  vout.general(m_vl, "plaq_SF with boundary spatial plaq = %.8f\n",
1155  plaq2 / (3 * 6 * Lx * Ly * Lz * Lt));
1156 
1157 #pragma omp barrier
1158 }
1159 
1160 
1161 //============================================================END=====
Staple_SF::init
void init()
Definition: staple_SF.cpp:69
Staple_SF::m_index
Index_lex m_index
Definition: staple_SF.h:41
CommonParameters::Ny
static int Ny()
Definition: commonParameters.h:106
Field_SF::set_boundary_wkpr
void set_boundary_wkpr(Field_G &u, const Mat_SU_N &wkpr)
Definition: field_SF.cpp:63
CommonParameters::Nz
static int Nz()
Definition: commonParameters.h:107
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
Staple_SF::set_parameters
void set_parameters(const Parameters &params)
Definition: staple_SF.cpp:112
Staple_SF::staple_ct
void staple_ct(Field_G &, const Field_G &, const int, const double ct)
Definition: staple_SF.cpp:990
ShiftField_lex::forward
void forward(Field &, const Field &, const int mu)
Definition: shiftField_lex.cpp:79
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Field::set
void set(const int jin, const int site, const int jex, double v)
Definition: field.h:175
Parameters
Class for parameters.
Definition: parameters.h:46
Field_G::mat_dag
Mat_SU_N mat_dag(const int site, const int mn=0) const
Definition: field_G.h:127
Staple_SF::tidyup
void tidyup()
Definition: staple_SF.cpp:100
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:518
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Staple_SF::m_stpl_lower
Field_G * m_stpl_lower
Definition: staple_SF.h:56
CommonParameters::Ly
static int Ly()
Definition: commonParameters.h:92
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Staple_SF::m_shift
ShiftField_lex m_shift
Definition: staple_SF.h:43
Staple_SF::plaquette
double plaquette(const Field_G &)
Definition: staple_SF.cpp:733
Staple_SF::sf_coupling_rect
double sf_coupling_rect(const Field_G &, const double ctr)
Definition: staple_SF.cpp:507
Field_SF::set_boundary_wk
void set_boundary_wk(Field_G &u, const Mat_SU_N &wk)
Definition: field_SF.cpp:32
Staple_SF::m_wkpr
Mat_SU_N m_wkpr
Definition: staple_SF.h:45
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
Staple_SF::m_Umu
Field_G * m_Umu
Definition: staple_SF.h:50
Parameters::set_double_vector
void set_double_vector(const string &key, const vector< double > &value)
Definition: parameters.cpp:42
Staple_SF::staple
void staple(Field_G &, const Field_G &, const int)
Definition: staple_SF.cpp:940
Staple_SF::upper
void upper(Field_G &, const Field_G &, const int, const int)
Definition: staple_SF.cpp:1039
Staple_SF::m_phi
std::vector< double > m_phi
Definition: staple_SF.h:48
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
Staple_SF::plaquette_ct
double plaquette_ct(const Field_G &, const double ct)
Definition: staple_SF.cpp:761
Staple_SF::m_i_omega0
Mat_SU_N m_i_omega0
Definition: staple_SF.h:45
SU_N::Mat_SU_N::set
void set(int c, const double &re, const double &im)
Definition: mat_SU_N.h:137
Field_SF::mult_ct_boundary
void mult_ct_boundary(Field_G &u, const int t, const double ct)
Definition: field_SF.cpp:185
CommonParameters::Nx
static int Nx()
Definition: commonParameters.h:105
Field_SF::set_boundary_zero
void set_boundary_zero(Field_G &u)
Definition: field_SF.cpp:96
Staple_SF::m_initialized
int m_initialized
Definition: staple_SF.h:46
CommonParameters::Lt
static int Lt()
Definition: commonParameters.h:94
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
Staple_SF::sf_coupling_plaq
double sf_coupling_plaq(const Field_G &, const double ct)
Definition: staple_SF.cpp:306
SU_N::Mat_SU_N::zero
Mat_SU_N & zero()
Definition: mat_SU_N.h:429
CommonParameters::Lz
static int Lz()
Definition: commonParameters.h:93
CommonParameters::Nt
static int Nt()
Definition: commonParameters.h:108
SU_N::Mat_SU_N
Definition: mat_SU_N.h:36
ThreadManager::reduce_sum_global
static void reduce_sum_global(dcomplex &value, const int i_thread, const int Nthread)
global reduction with summation: dcomplex values are assumed thread local.
Definition: threadManager.cpp:288
threadManager.h
Staple_SF::m_phipr
std::vector< double > m_phipr
Definition: staple_SF.h:48
Index_lex::site
int site(const int &x, const int &y, const int &z, const int &t) const
Definition: index_lex.h:55
Staple_SF::m_p_omega
std::vector< double > m_p_omega
Definition: staple_SF.h:48
Field_SF::set_boundary_matrix
void set_boundary_matrix(Mat_SU_N &wk, const std::vector< double > &phi)
Definition: field_SF.cpp:216
Staple_SF::plaq_t_ct
double plaq_t_ct(const Field_G &, const double ct)
Definition: staple_SF.cpp:875
Staple_SF::m_wk
Mat_SU_N m_wk
Definition: staple_SF.h:45
staple_SF.h
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Staple_SF::lower
void lower(Field_G &, const Field_G &, const int, const int)
Definition: staple_SF.cpp:1081
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
Staple_SF::class_name
static const std::string class_name
Definition: staple_SF.h:35
mult_Field_Gnn
void mult_Field_Gnn(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:95
ShiftField_lex::backward
void backward(Field &, const Field &, const int mu)
Definition: shiftField_lex.cpp:59
Staple_SF::m_vt
Field_G * m_vt
Definition: staple_SF.h:52
Staple_SF::m_Unu
Field_G * m_Unu
Definition: staple_SF.h:51
Staple_SF::m_vl
Bridge::VerboseLevel m_vl
Definition: staple_SF.h:38
Communicator::ipe
static int ipe(const int dir)
logical coordinate of current proc.
Definition: communicator.cpp:105
Staple_SF::print_plaquette
void print_plaquette(const Field_G &)
Definition: staple_SF.cpp:1140
Parameters::fetch_string
int fetch_string(const string &key, string &value) const
Definition: parameters.cpp:378
field_thread-inc.h
CommonParameters::NPEt
static int NPEt()
Definition: commonParameters.h:100
ixyz
int ixyz
Definition: mult_Domainwall_eo_t_dirac_openacc-inc.h:13
SU_N::Mat_SU_N::dag
Mat_SU_N & dag()
Definition: mat_SU_N.h:329
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
SU_N::ReTr
double ReTr(const Mat_SU_N &m)
Definition: mat_SU_N.h:534
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
Field_G::mat
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:114
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Staple_SF::plaq_t
double plaq_t(const Field_G &)
Definition: staple_SF.cpp:825
Parameters::fetch_double_vector
int fetch_double_vector(const string &key, vector< double > &value) const
Definition: parameters.cpp:410
Staple_SF::m_stpl_upper
Field_G * m_stpl_upper
Definition: staple_SF.h:55
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
Staple_SF::get_parameters
void get_parameters(Parameters &params) const
Definition: staple_SF.cpp:147
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
Staple_SF::m_wt
Field_G * m_wt
Definition: staple_SF.h:53
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Staple_SF::plaq_s
double plaq_s(const Field_G &)
Definition: staple_SF.cpp:783
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