Bridge++  Ver.2.1.3
forceSmear_APE_SF.cpp
Go to the documentation of this file.
1 
12 
13 #ifdef USE_FACTORY_AUTOREGISTER
14 namespace {
15  bool init = ForceSmear_APE_SF::register_factory();
16 }
17 #endif
18 
19 const std::string ForceSmear_APE_SF::class_name = "ForceSmear_APE_SF";
20 
21 //====================================================================
23 {
25 
26  std::string vlevel;
27  if (!params.fetch_string("verbose_level", vlevel)) {
28  m_vl = vout.set_verbose_level(vlevel);
29  } else {
31  }
32 
33  vout.general(m_vl, "%s: construction\n", class_name.c_str());
35 
38 
39  m_rho.resize(m_Ndim * m_Ndim);
40  for (int mu = 0; mu < m_Ndim; ++mu) {
41  for (int nu = 0; nu < m_Ndim; ++nu) {
42  m_rho[index(mu, nu)] = 0.0;
43  }
44  }
45 
46  set_parameters(params);
47 
48  m_shift = new ShiftField_lex();
49 
50  m_U.resize(m_Ndim);
51  m_iTheta.resize(m_Ndim);
52 
53  m_vt1 = new Field_G(m_Nvol, 1);
54  m_vt2 = new Field_G(m_Nvol, 1);
55  m_vt3 = new Field_G(m_Nvol, 1);
56  m_v1 = new Field_G(m_Nvol, 1);
57  m_v2 = new Field_G(m_Nvol, 1);
58 
59  m_ct1 = new Field_G(m_Nvol, 1);
60  m_ct2 = new Field_G(m_Nvol, 1);
61  m_xi = new Field_G(m_Nvol, 1);
62 
64  vout.general(m_vl, "%s: construction finished.\n",
65  class_name.c_str());
66 
67 }
68 
69 //====================================================================
71 {
72  delete m_shift;
73 
74  delete m_vt1;
75  delete m_vt2;
76  delete m_vt3;
77  delete m_v1;
78  delete m_v2;
79 
80  delete m_ct1;
81  delete m_ct2;
82  delete m_xi;
83 }
84 
85 //====================================================================
87 {
88  std::string vlevel;
89  if (!params.fetch_string("verbose_level", vlevel)) {
90  m_vl = vout.set_verbose_level(vlevel);
91  }
92 
93  //- fetch and check input parameters
94  double rho1;
95  std::vector<double> phi, phipr;
96 
97  int err = 0;
98  err += params.fetch_double("rho_uniform", rho1);
99  err += params.fetch_double_vector("phi", phi);
100  err += params.fetch_double_vector("phipr", phipr);
101 
102  if (err) {
103  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
104  class_name.c_str());
105  exit(EXIT_FAILURE);
106  }
107 
108  set_parameters(rho1, phi, phipr);
109 }
110 
111 
112 //====================================================================
113 void ForceSmear_APE_SF::set_parameters(const double rho1,
114  const std::vector<double>& phi,
115  const std::vector<double>& phipr)
116 {
117  int ith = ThreadManager::get_thread_id();
118  if (ith == 0) {
119 
120  for (int mu = 0; mu < m_Ndim; ++mu) {
121  for (int nu = 0; nu < m_Ndim; ++nu) {
122  m_rho[index(mu, nu)] = rho1;
123  }
124  }
125 
126  m_phi.resize(3);
127  m_phipr.resize(3);
128  for (int i = 0; i < 3; ++i) {
129  m_phi[i] = phi[i];
130  m_phipr[i] = phipr[i];
131  }
132 
133  const int Lx = CommonParameters::Lx();
134 
135  double Lx_inv = 1.0 / double(Lx);
136  double c0r = cos(phi[0] * Lx_inv);
137  double c0i = sin(phi[0] * Lx_inv);
138  double c1r = cos(phi[1] * Lx_inv);
139  double c1i = sin(phi[1] * Lx_inv);
140  double c2r = cos(phi[2] * Lx_inv);
141  double c2i = sin(phi[2] * Lx_inv);
142 
143  m_wk.zero();
144  m_wk.set(0, 0, c0r, c0i);
145  m_wk.set(1, 1, c1r, c1i);
146  m_wk.set(2, 2, c2r, c2i);
147 
148  c0r = cos(phipr[0] * Lx_inv);
149  c0i = sin(phipr[0] * Lx_inv);
150  c1r = cos(phipr[1] * Lx_inv);
151  c1i = sin(phipr[1] * Lx_inv);
152  c2r = cos(phipr[2] * Lx_inv);
153  c2i = sin(phipr[2] * Lx_inv);
154 
155  m_wkpr.zero();
156  m_wkpr.set(0, 0, c0r, c0i);
157  m_wkpr.set(1, 1, c1r, c1i);
158  m_wkpr.set(2, 2, c2r, c2i);
159  }
160 #pragma omp barrier
161 
162  vout.general(m_vl, "%s: parameters\n", class_name.c_str());
163  vout.general(m_vl, " rho = %8.4f\n", rho1);
164 
165  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
166  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
167  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
168  vout.general(m_vl, " phipr1= %12.6f\n", phipr[0]);
169  vout.general(m_vl, " phipr2= %12.6f\n", phipr[1]);
170  vout.general(m_vl, " phipr3= %12.6f\n", phipr[2]);
171 
172 }
173 
174 
175 //====================================================================
176 void ForceSmear_APE_SF::set_parameters(const std::vector<double>& rho,
177  const std::vector<double>& phi,
178  const std::vector<double>& phipr)
179 {
180  int ith = ThreadManager::get_thread_id();
181  if (ith == 0) {
182 
183  assert(rho.size() == m_Ndim * m_Ndim);
184 
185  for (int mu = 0; mu < m_Ndim; ++mu) {
186  for (int nu = 0; nu < m_Ndim; ++nu) {
187  m_rho[index(mu, nu)] = rho[index(mu, nu)];
188  }
189  }
190 
191  m_phi.resize(3);
192  m_phipr.resize(3);
193  for (int i = 0; i < 3; ++i) {
194  m_phi[i] = phi[i];
195  m_phipr[i] = phipr[i];
196  }
197 
198  const int Lx = CommonParameters::Lx();
199 
200  double Lx_inv = 1.0 / double(Lx);
201  double c0r = cos(phi[0] * Lx_inv);
202  double c0i = sin(phi[0] * Lx_inv);
203  double c1r = cos(phi[1] * Lx_inv);
204  double c1i = sin(phi[1] * Lx_inv);
205  double c2r = cos(phi[2] * Lx_inv);
206  double c2i = sin(phi[2] * Lx_inv);
207  m_wk.zero();
208  m_wk.set(0, 0, c0r, c0i);
209  m_wk.set(1, 1, c1r, c1i);
210  m_wk.set(2, 2, c2r, c2i);
211 
212  c0r = cos(phipr[0] * Lx_inv);
213  c0i = sin(phipr[0] * Lx_inv);
214  c1r = cos(phipr[1] * Lx_inv);
215  c1i = sin(phipr[1] * Lx_inv);
216  c2r = cos(phipr[2] * Lx_inv);
217  c2i = sin(phipr[2] * Lx_inv);
218 
219  m_wkpr.zero();
220  m_wkpr.set(0, 0, c0r, c0i);
221  m_wkpr.set(1, 1, c1r, c1i);
222  m_wkpr.set(2, 2, c2r, c2i);
223  }
224 #pragma omp barrier
225 
226  vout.general(m_vl, "%s:\n", class_name.c_str());
227  for (int mu = 0; mu < m_Ndim; ++mu) {
228  vout.general(m_vl, " rho[%d] = %8.4f\n", mu, rho[mu]);
229  }
230  vout.general(m_vl, " phi1 = %12.6f\n", phi[0]);
231  vout.general(m_vl, " phi2 = %12.6f\n", phi[1]);
232  vout.general(m_vl, " phi3 = %12.6f\n", phi[2]);
233  vout.general(m_vl, " phipr1= %12.6f\n", phipr[0]);
234  vout.general(m_vl, " phipr2= %12.6f\n", phipr[1]);
235  vout.general(m_vl, " phipr3= %12.6f\n", phipr[2]);
236 
237 }
238 
239 
240 //====================================================================
242 {
243  params.set_double_vector("rho", m_rho);
244  params.set_double_vector("phi", m_phi);
245  params.set_double_vector("phipr", m_phipr);
246 
247  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
248 }
249 
250 
251 //====================================================================
252 
265 //====================================================================
267  const Field_G& Sigmap,
268  const Field_G& U)
269 {
270 #pragma omp barrier
271 
272  const int Nc = CommonParameters::Nc();
273 
274  assert(Sigmap.nin() == (2 * Nc * Nc));
275  assert(Sigmap.nvol() == m_Nvol);
276  assert(Sigmap.nex() == m_Ndim);
277 
278  for (int mu = 0; mu < m_Ndim; ++mu) {
279  copy(m_U[mu], 0, U, mu);
280 #pragma omp barrier
281 
282  if (mu != 3) Field_SF::set_boundary_wk(m_U[mu], m_wk);
283 #pragma omp barrier
284  }
285 
286  for (int mu = 0; mu < m_Ndim; ++mu) {
287  m_ct1->set(0.0);
288 #pragma omp barrier
289 
290  for (int nu = 0; nu < m_Ndim; ++nu) {
291  if (nu == mu) continue;
292  double rho = m_rho[index(mu, nu)];
293 
294  staple(*m_ct2, m_U[mu], m_U[nu], mu, nu);
295  axpy(*m_ct1, 0, rho, *m_ct2, 0);
296 #pragma omp barrier
297  }
298 
299  copy(*m_ct2, 0, Sigmap, mu);
300 #pragma omp barrier
301 
302  double alpha = m_rho[mu + m_Ndim * mu];
303 
305  alpha, *m_ct2, *m_ct1, m_U[mu]);
306  copy(Sigma, mu, *m_xi, 0);
307 #pragma omp barrier
308  }
309 
310  for (int mu = 0; mu < m_Ndim; ++mu) {
311  for (int nu = 0; nu < m_Ndim; ++nu) {
312  if (nu == mu) continue;
313 
314  double rho = m_rho[mu + m_Ndim * nu];
315 
316  force_each(*m_ct2, m_U[mu], m_U[nu],
317  m_iTheta[mu], m_iTheta[nu], mu, nu);
318  axpy(Sigma, mu, rho, *m_ct2, 0);
319 #pragma omp barrier
320  }
321 
322  }
324 #pragma omp barrier
325 
326 }
327 
328 
329 //====================================================================
330 
345  const Field_G& V_mu, const Field_G& V_nu,
346  const Field_G& iTheta_mu,
347  const Field_G& iTheta_nu,
348  const int mu, const int nu)
349 {
350 #pragma omp barrier
351 
352  Sigma_mu.set(0.0);
353 #pragma omp barrier
354 
355  //- The 1st block
356  m_shift->backward(*m_vt1, V_nu, mu);
357  if (mu == 3) Field_SF::set_boundary_wkpr(*m_vt1, m_wkpr);
358 #pragma omp barrier
359 
360  m_shift->backward(*m_vt2, V_mu, nu);
361  if (nu == 3) Field_SF::set_boundary_wkpr(*m_vt2, m_wkpr);
362 #pragma omp barrier
363 
364  mult_Field_Gnd(*m_vt3, 0, *m_vt1, 0, *m_vt2, 0);
365  multadd_Field_Gnd(Sigma_mu, 0, *m_vt3, 0, iTheta_nu, 0, 1.0);
366 #pragma omp barrier
367 
368  //- The 2nd block
369  mult_Field_Gdn(*m_vt3, 0, iTheta_mu, 0, V_nu, 0);
370  mult_Field_Gdn(*m_vt2, 0, *m_vt1, 0, *m_vt3, 0);
371  m_shift->forward(*m_vt3, *m_vt2, nu);
372  axpy(Sigma_mu, 1.0, *m_vt3);
373 #pragma omp barrier
374 
375  //- The 3rd block
376  mult_Field_Gdn(*m_vt3, 0, V_mu, 0, iTheta_nu, 0);
377  mult_Field_Gdn(*m_vt2, 0, *m_vt1, 0, *m_vt3, 0);
378  m_shift->forward(*m_vt3, *m_vt2, nu);
379  axpy(Sigma_mu, 1.0, *m_vt3);
380 #pragma omp barrier
381 
382  //- The 4th block
383  m_shift->backward(*m_vt1, iTheta_nu, mu);
384  m_shift->backward(*m_vt2, V_mu, nu);
385  if (nu == 3) Field_SF::set_boundary_wkpr(*m_vt2, m_wkpr);
386 #pragma omp barrier
387 
388  mult_Field_Gnd(*m_vt3, 0, *m_vt1, 0, *m_vt2, 0);
389  multadd_Field_Gnd(Sigma_mu, 0, *m_vt3, 0, V_nu, 0, 1.0);
390 #pragma omp barrier
391 
392  //- The 5th block
393  mult_Field_Gdd(*m_vt2, 0, *m_vt1, 0, V_mu, 0);
394  mult_Field_Gnn(*m_vt3, 0, *m_vt2, 0, V_nu, 0);
395  m_shift->forward(*m_vt2, *m_vt3, nu);
396  axpy(Sigma_mu, 1.0, *m_vt2);
397 #pragma omp barrier
398 
399  //- The 6th block
400  m_shift->backward(*m_vt1, V_nu, mu);
401  if (mu == 3) Field_SF::set_boundary_wkpr(*m_vt1, m_wkpr);
402 #pragma omp barrier
403 
404  m_shift->backward(*m_vt2, iTheta_mu, nu);
405  mult_Field_Gnd(*m_vt3, 0, *m_vt1, 0, *m_vt2, 0);
406  multadd_Field_Gnd(Sigma_mu, 0, *m_vt3, 0, V_nu, 0, 1.0);
407 #pragma omp barrier
408 
409 }
410 
411 
412 //====================================================================
414  const Field_G& u_mu, const Field_G& u_nu,
415  const int mu, const int nu)
416 {
417 #pragma omp barrier
418 
419  //- upper direction
420  m_shift->backward(*m_v1, u_mu, nu);
421  if (nu == 3) Field_SF::set_boundary_wkpr(*m_v1, m_wkpr);
422 #pragma omp barrier
423 
424  mult_Field_Gnn(*m_v2, 0, u_nu, 0, *m_v1, 0);
425 #pragma omp barrier
426 
427  m_shift->backward(*m_v1, u_nu, mu);
428  if (mu == 3) Field_SF::set_boundary_wkpr(*m_v1, m_wkpr);
429 #pragma omp barrier
430 
431  mult_Field_Gnd(c, 0, *m_v2, 0, *m_v1, 0);
432 #pragma omp barrier
433 
434  //- lower direction
435  m_shift->backward(*m_v2, u_nu, mu);
436  if (mu == 3) Field_SF::set_boundary_wkpr(*m_v2, m_wkpr);
437 #pragma omp barrier
438 
439  mult_Field_Gnn(*m_v1, 0, u_mu, 0, *m_v2, 0);
440  mult_Field_Gdn(*m_v2, 0, u_nu, 0, *m_v1, 0);
441  m_shift->forward(*m_v1, *m_v2, nu);
442  axpy(c, 0, 1.0, *m_v1, 0);
443 #pragma omp barrier
444 
445 }
446 
447 //============================================================END=====
Field_SF::set_boundary_wkpr
void set_boundary_wkpr(Field_G &u, const Mat_SU_N &wkpr)
Definition: field_SF.cpp:63
ForceSmear_APE_SF::m_U
std::vector< Field_G > m_U
Definition: forceSmear_APE_SF.h:48
ForceSmear_APE_SF::m_v1
Field_G * m_v1
Definition: forceSmear_APE_SF.h:54
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
ForceSmear_APE_SF::m_phi
std::vector< double > m_phi
SF boundary condition at t=0.
Definition: forceSmear_APE_SF.h:40
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
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.cpp:518
ForceSmear_APE_SF::m_xi
Field_G * m_xi
Definition: forceSmear_APE_SF.h:62
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.cpp:508
Field::nex
int nex() const
Definition: field.h:128
ForceSmear_APE_SF::m_wk
Mat_SU_N m_wk
SF boundary condition at t=0.
Definition: forceSmear_APE_SF.h:51
ForceSmear_APE_SF::force_udiv
void force_udiv(Field_G &Sigma, const Field_G &Sigma_p, const Field_G &U)
Definition: forceSmear_APE_SF.cpp:266
CommonParameters::Nvol
static int Nvol()
Definition: commonParameters.h:109
Field_SF::set_boundary_wk
void set_boundary_wk(Field_G &u, const Mat_SU_N &wk)
Definition: field_SF.cpp:32
ForceSmear_APE_SF::get_parameters
void get_parameters(Parameters &params) const
Definition: forceSmear_APE_SF.cpp:241
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:381
ForceSmear_APE_SF::m_v2
Field_G * m_v2
Definition: forceSmear_APE_SF.h:55
ForceSmear_APE_SF::init
void init(const Parameters &params)
Definition: forceSmear_APE_SF.cpp:22
ForceSmear_APE_SF::set_parameters
void set_parameters(const Parameters &params)
Definition: forceSmear_APE_SF.cpp:86
Parameters::set_double_vector
void set_double_vector(const string &key, const vector< double > &value)
Definition: parameters.cpp:42
Field::nin
int nin() const
Definition: field.h:126
AProjection::force_recursive
virtual void force_recursive(Field_G &Xi, Field_G &iTheta, const double alpha, const Field_G &Sigmap, const Field_G &C, const Field_G &U)=0
determination of fields for force calculation
ForceSmear_APE_SF::m_rho
std::vector< double > m_rho
Definition: forceSmear_APE_SF.h:39
ForceSmear_APE_SF::m_ct2
Field_G * m_ct2
Definition: forceSmear_APE_SF.h:61
copy
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:213
SU_N::Mat_SU_N::set
void set(int c, const double &re, const double &im)
Definition: mat_SU_N.h:137
ForceSmear_APE_SF::m_ct1
Field_G * m_ct1
Definition: forceSmear_APE_SF.h:60
ForceSmear_APE_SF::m_Nvol
int m_Nvol
Definition: forceSmear_APE_SF.h:38
CommonParameters::Lx
static int Lx()
Definition: commonParameters.h:91
CommonParameters::Nc
static int Nc()
Definition: commonParameters.h:115
ForceSmear_APE_SF::m_vt3
Field_G * m_vt3
Definition: forceSmear_APE_SF.h:58
ForceSmear_APE_SF::class_name
static const std::string class_name
Definition: forceSmear_APE_SF.h:35
SU_N::Mat_SU_N::zero
Mat_SU_N & zero()
Definition: mat_SU_N.h:429
ForceSmear_APE_SF::m_vt1
Field_G * m_vt1
Definition: forceSmear_APE_SF.h:56
ForceSmear_APE_SF::m_iTheta
std::vector< Field_G > m_iTheta
Definition: forceSmear_APE_SF.h:49
ForceSmear_APE_SF::m_vt2
Field_G * m_vt2
Definition: forceSmear_APE_SF.h:57
ForceSmear_APE_SF::m_wkpr
Mat_SU_N m_wkpr
SF boundary condition at t=Nt.
Definition: forceSmear_APE_SF.h:52
threadManager.h
ForceSmear_APE_SF::m_proj
Projection * m_proj
Definition: forceSmear_APE_SF.h:44
Field::nvol
int nvol() const
Definition: field.h:127
ForceSmear_APE_SF::m_shift
ShiftField_lex * m_shift
Definition: forceSmear_APE_SF.h:46
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
Field_SF::set_boundary_spatial_link_zero
void set_boundary_spatial_link_zero(Field_G &u)
Definition: field_SF.cpp:151
ForceSmear_APE_SF::m_Ndim
int m_Ndim
Definition: forceSmear_APE_SF.h:38
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
forceSmear_APE_SF.h
ForceSmear_APE_SF::m_phipr
std::vector< double > m_phipr
SF boundary condition at t=Nt.
Definition: forceSmear_APE_SF.h:41
ForceSmear_APE_SF::index
int index(int mu, int nu)
Definition: forceSmear_APE_SF.h:93
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:195
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
ForceSmear_APE_SF::m_vl
Bridge::VerboseLevel m_vl
Definition: forceSmear_APE_SF.h:42
mult_Field_Gdd
void mult_Field_Gdd(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:212
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
Bridge::BridgeIO::crucial
void crucial(const char *format,...)
Definition: bridgeIO.cpp:242
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
ForceSmear_APE_SF::tidyup
void tidyup()
Definition: forceSmear_APE_SF.cpp:70
Field_G
SU(N) gauge field.
Definition: field_G.h:38
Parameters::fetch_double_vector
int fetch_double_vector(const string &key, vector< double > &value) const
Definition: parameters.cpp:410
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:262
mult_Field_Gnd
void mult_Field_Gnd(Field_G &W, const int ex, const Field_G &U1, const int ex1, const Field_G &U2, const int ex2)
Definition: field_G_imp.cpp:173
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
ForceSmear_APE_SF::force_each
void force_each(Field_G &, const Field_G &, const Field_G &, const Field_G &, const Field_G &, const int mu, const int nu)
Definition: forceSmear_APE_SF.cpp:344
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:572
ForceSmear_APE_SF::staple
void staple(Field_G &, const Field_G &, const Field_G &, const int mu, const int nu)
Definition: forceSmear_APE_SF.cpp:413
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:216