Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
projection_Stout_SU3.cpp
Go to the documentation of this file.
1 
14 #include "projection_Stout_SU3.h"
15 
16 #ifdef USE_FACTORY_AUTOREGISTER
17 namespace {
18  bool init = Projection_Stout_SU3::register_factory();
19 }
20 #endif
21 
22 const std::string Projection_Stout_SU3::class_name = "Projection_Stout_SU3";
23 
24 //====================================================================
26 {
27  const string str_vlevel = params.get_string("verbose_level");
28 
29  m_vl = vout.set_verbose_level(str_vlevel);
30 }
31 
32 
33 //====================================================================
35 {
36  m_flop = 0;
37  m_time = 0.0;
38 
39  assert(CommonParameters::Nc() == NC);
40 }
41 
42 
43 //====================================================================
45 {
46  const double gflops = 1.0e-9 * double(m_flop) / m_time;
47 
48  vout.general(m_vl, "%s:\n", class_name.c_str());
49  vout.general(m_vl, " total time: %f\n", m_time);
50  vout.general(m_vl, " total flop: %d\n", m_flop);
51  vout.general(m_vl, " GFlops : %f\n", gflops);
52 }
53 
54 
55 //====================================================================
57  const double alpha,
58  const Field_G& Cst, const Field_G& Uorg)
59 {
60  // int id = 31;
61  // KEK_FopCountStart(id);
62  const double time0 = Communicator::get_time();
63 
64  // in stout projection, parameter alpha is dummy.
65 
66  const int Nex = Uorg.nex();
67  const int Nvol = Uorg.nvol();
68  const int NinG = Uorg.nin();
69 
70  assert(Cst.nex() == Nex);
71  assert(Cst.nvol() == Nvol);
72  assert(U.nex() == Nex);
73  assert(U.nvol() == Nvol);
74 
75  Mat_SU_N iQ0(NC);
76  iQ0.unit();
77 
78  for (int mu = 0; mu < Nex; ++mu) {
79  for (int site = 0; site < Nvol; ++site) {
80  Mat_SU_N ut(NC);
81  Uorg.mat(ut, site, mu);
82 
83  Mat_SU_N ct(NC);
84  Cst.mat(ct, site, mu);
85 
86  Mat_SU_N iQ1(NC);
87  iQ1.mult_nd(ct, ut);
88  iQ1.at();
89 
90  Mat_SU_N iQ2(NC);
91  iQ2.mult_nn(iQ1, iQ1);
92 
93  Mat_SU_N iQ3(NC);
94  iQ3.mult_nn(iQ1, iQ2);
95 
96  Mat_SU_N e_iQ(NC);
97 
98  double norm = iQ1.norm2();
99  if (norm > 1.0e-10) {
100  double u, w;
101  set_uw(u, w, iQ2, iQ3);
102 
103  dcomplex f0, f1, f2;
104  set_fj(f0, f1, f2, u, w);
105 
106  for (int cc = 0; cc < NC * NC; ++cc) {
107  dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
108  + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
109  - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
110  e_iQ.set_r(cc, real(qt));
111  e_iQ.set_i(cc, imag(qt));
112  }
113  } else {
114  // vout.general(m_vl,"project: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
115  e_iQ.unit();
116  }
117 
118  Mat_SU_N ut2(NC);
119  ut2.mult_nn(e_iQ, ut);
120  U.set_mat(site, mu, ut2);
121  }
122  }
123 
124  /*
125  unsigned long count;
126  double time;
127  KEK_FopCountFinish(id,&count,&time);
128  m_time += time;
129  m_flop += count;
130  */
131  const double time1 = Communicator::get_time();
132  m_time += time1 - time0;
133 }
134 
135 
136 //====================================================================
138 {
139  const int Nvol = iQ.nvol();
140  const int Nex = iQ.nex();
141 
142  Mat_SU_N iQ0(NC);
143 
144  iQ0.unit();
145 
146  for (int mu = 0; mu < Nex; ++mu) {
147  for (int site = 0; site < Nvol; ++site) {
148  Mat_SU_N iQ1 = iQ.mat(site, mu);
149  Mat_SU_N iQ2 = iQ1 * iQ1;
150  Mat_SU_N iQ3 = iQ1 * iQ2;
151 
152  double norm = iQ1.norm2();
153  if (norm > 1.0e-10) {
154  double u, w;
155  set_uw(u, w, iQ2, iQ3);
156 
157  dcomplex f0, f1, f2;
158  set_fj(f0, f1, f2, u, w);
159 
160  for (int cc = 0; cc < NC * NC; ++cc) {
161  dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
162  + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
163  - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
164  e_iQ.set_ri(cc, site, mu, real(qt), imag(qt));
165  }
166  } else {
167  // vout.general(m_vl,"exp_iQ: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
168  e_iQ.set_mat(site, mu, iQ0);
169  }
170  }
171  }
172 }
173 
174 
175 //====================================================================
176 
189 //====================================================================
191  const double alpha, const Field_G& Sigmap,
192  const Field_G& Cst, const Field_G& Uorg)
193 {
194  // in stout projection, parameter alpha is dummy.
195 
196  // int id = 31;
197  // KEK_FopCountStart(id);
198  const double time0 = Communicator::get_time();
199 
200  const int Nvol = CommonParameters::Nvol();
201  // int NinG = 2 * NC * NC;
202  const int Nex = Xi.nex();
203 
204  assert(Xi.nvol() == Nvol);
205  assert(iTheta.nvol() == Nvol);
206  assert(Sigmap.nvol() == Nvol);
207  assert(Cst.nvol() == Nvol);
208  assert(Uorg.nvol() == Nvol);
209  assert(iTheta.nex() == Nex);
210  assert(Sigmap.nex() == Nex);
211  assert(Cst.nex() == Nex);
212  assert(Uorg.nex() == Nex);
213 
214  Mat_SU_N iQ0(NC);
215  iQ0.unit();
216 
217  for (int mu = 0; mu < Nex; ++mu) {
218  for (int site = 0; site < Nvol; ++site) {
220  Mat_SU_N C_tmp(NC);
221  Cst.mat(C_tmp, site, mu);
222 
224  Mat_SU_N U_tmp(NC);
225  Uorg.mat(U_tmp, site, mu);
226 
227  // Sigmap_tmp \f$=\Sigma_\mu'(x)\f$
228  Mat_SU_N Sigmap_tmp(NC);
229  Sigmap.mat(Sigmap_tmp, site, mu);
230 
232  Mat_SU_N iQ1(NC);
233  iQ1.mult_nd(C_tmp, U_tmp);
234  iQ1.at();
235 
236  Mat_SU_N iQ2(NC);
237  iQ2.mult_nn(iQ1, iQ1);
238 
239  Mat_SU_N iQ3(NC);
240  iQ3.mult_nn(iQ1, iQ2);
241 
242  // In order to aviod 1Q1=0
243  Mat_SU_N e_iQ(NC), iGamma(NC);
244 
245  double norm = iQ1.norm2();
246  if (norm > 1.0e-10) {
247  double u, w;
248  set_uw(u, w, iQ2, iQ3);
249 
250  dcomplex f0, f1, f2;
251  set_fj(f0, f1, f2, u, w);
252 
253  for (int cc = 0; cc < NC * NC; ++cc) {
254  dcomplex qt = f0 * cmplx(iQ0.r(cc), iQ0.i(cc))
255  + f1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
256  - f2 * cmplx(iQ2.r(cc), iQ2.i(cc));
257  e_iQ.set(cc, real(qt), imag(qt));
258  }
259 
260  double xi0 = func_xi0(w);
261  double xi1 = func_xi1(w);
262  double u2 = u * u;
263  double w2 = w * w;
264  double cos_w = cos(w);
265 
266  dcomplex emiu = cmplx(cos(u), -sin(u));
267  dcomplex e2iu = cmplx(cos(2.0 * u), sin(2.0 * u));
268 
269  dcomplex r01 = cmplx(2.0 * u, 2.0 * (u2 - w2)) * e2iu
270  + emiu * cmplx(16.0 * u * cos_w + 2.0 * u * (3.0 * u2 + w2) * xi0,
271  -8.0 * u2 * cos_w + 2.0 * (9.0 * u2 + w2) * xi0);
272 
273  dcomplex r11 = cmplx(2.0, 4.0 * u) * e2iu
274  + emiu * cmplx(-2.0 * cos_w + (3.0 * u2 - w2) * xi0,
275  2.0 * u * cos_w + 6.0 * u * xi0);
276 
277  dcomplex r21 = cmplx(0.0, 2.0) * e2iu
278  + emiu * cmplx(-3.0 * u * xi0, cos_w - 3.0 * xi0);
279 
280  dcomplex r02 = cmplx(-2.0, 0.0) * e2iu
281  + emiu * cmplx(-8.0 * u2 * xi0,
282  2.0 * u * (cos_w + xi0 + 3.0 * u2 * xi1));
283 
284  dcomplex r12 = emiu * cmplx(2.0 * u * xi0,
285  -cos_w - xi0 + 3.0 * u2 * xi1);
286 
287  dcomplex r22 = emiu * cmplx(xi0, -3.0 * u * xi1);
288 
289  double fden = 1.0 / (2 * (9.0 * u2 - w2) * (9.0 * u2 - w2));
290 
291  dcomplex b10 = cmplx(2.0 * u, 0.0) * r01 + cmplx(3.0 * u2 - w2, 0.0) * r02
292  - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f0;
293  dcomplex b11 = cmplx(2.0 * u, 0.0) * r11 + cmplx(3.0 * u2 - w2, 0.0) * r12
294  - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f1;
295  dcomplex b12 = cmplx(2.0 * u, 0.0) * r21 + cmplx(3.0 * u2 - w2, 0.0) * r22
296  - cmplx(30.0 * u2 + 2.0 * w2, 0.0) * f2;
297 
298  dcomplex b20 = r01 - cmplx(3.0 * u, 0.0) * r02 - cmplx(24.0 * u, 0.0) * f0;
299  dcomplex b21 = r11 - cmplx(3.0 * u, 0.0) * r12 - cmplx(24.0 * u, 0.0) * f1;
300  dcomplex b22 = r21 - cmplx(3.0 * u, 0.0) * r22 - cmplx(24.0 * u, 0.0) * f2;
301 
302  b10 *= cmplx(fden, 0.0);
303  b11 *= cmplx(fden, 0.0);
304  b12 *= cmplx(fden, 0.0);
305  b20 *= cmplx(fden, 0.0);
306  b21 *= cmplx(fden, 0.0);
307  b22 *= cmplx(fden, 0.0);
308 
309  Mat_SU_N B1(NC), B2(NC);
310  for (int cc = 0; cc < NC * NC; ++cc) {
311  dcomplex qt1 = b10 * cmplx(iQ0.r(cc), iQ0.i(cc))
312  + b11 * cmplx(iQ1.i(cc), -iQ1.r(cc))
313  - b12 * cmplx(iQ2.r(cc), iQ2.i(cc));
314  B1.set(cc, real(qt1), imag(qt1));
315 
316  dcomplex qt2 = b20 * cmplx(iQ0.r(cc), iQ0.i(cc))
317  + b21 * cmplx(iQ1.i(cc), -iQ1.r(cc))
318  - b22 * cmplx(iQ2.r(cc), iQ2.i(cc));
319  B2.set(cc, real(qt2), imag(qt2));
320  }
321 
322  Mat_SU_N USigmap(NC);
323  USigmap.mult_nn(U_tmp, Sigmap_tmp);
324 
325  Mat_SU_N tmp1(NC);
326  tmp1.mult_nn(USigmap, B1);
327 
328  Mat_SU_N tmp2(NC);
329  tmp2.mult_nn(USigmap, B2);
330 
331  dcomplex tr1 = cmplx(tmp1.r(0) + tmp1.r(4) + tmp1.r(8),
332  tmp1.i(0) + tmp1.i(4) + tmp1.i(8));
333  dcomplex tr2 = cmplx(tmp2.r(0) + tmp2.r(4) + tmp2.r(8),
334  tmp2.i(0) + tmp2.i(4) + tmp2.i(8));
335 
336  Mat_SU_N iQUS(NC);
337  iQUS.mult_nn(iQ1, USigmap);
338 
339  Mat_SU_N iUSQ(NC);
340  iUSQ.mult_nn(USigmap, iQ1);
341 
342  for (int cc = 0; cc < NC * NC; ++cc) {
343  dcomplex qt = tr1 * cmplx(iQ1.i(cc), -iQ1.r(cc))
344  - tr2 * cmplx(iQ2.r(cc), iQ2.i(cc))
345  + f1 * cmplx(USigmap.r(cc), USigmap.i(cc))
346  + f2 * cmplx(iQUS.i(cc), -iQUS.r(cc))
347  + f2 * cmplx(iUSQ.i(cc), -iUSQ.r(cc));
348  iGamma.set(cc, -imag(qt), real(qt));
349  }
350  } else {
351  // vout.general(m_vl,"force_recursive: |iQ1|^2 too small: %lf. Set e_iQ=1.\n",norm);
352  iGamma.zero();
353  e_iQ.unit();
354  }
355 
357  iGamma.at();
358 
359  Mat_SU_N iTheta_tmp(NC);
360  iTheta_tmp.mult_nn(iGamma, U_tmp);
362  iTheta.set_mat(site, mu, iTheta_tmp);
363 
364  Mat_SU_N Xi_tmp(NC);
365  Xi_tmp.mult_nn(Sigmap_tmp, e_iQ);
366  Xi_tmp.multadd_dn(C_tmp, iGamma);
368  Xi.set_mat(site, mu, Xi_tmp);
369  }
370  }
371 
372  /*
373  unsigned long count;
374  double time;
375  KEK_FopCountFinish(id,&count,&time);
376  m_time += time;
377  m_flop += count;
378  */
379  const double time1 = Communicator::get_time();
380  m_time += time1 - time0;
381 }
382 
383 
384 //====================================================================
385 void Projection_Stout_SU3::set_fj(dcomplex& f0, dcomplex& f1, dcomplex& f2,
386  const double& u, const double& w)
387 {
388  const double xi0 = func_xi0(w);
389  const double u2 = u * u;
390  const double w2 = w * w;
391  const double cos_w = cos(w);
392 
393  const double cos_u = cos(u);
394  const double sin_u = sin(u);
395 
396  const dcomplex emiu = cmplx(cos_u, -sin_u);
397  const dcomplex e2iu = cmplx(cos_u * cos_u - sin_u * sin_u, 2.0 * sin_u * cos_u);
398 
399  const dcomplex h0 = e2iu * cmplx(u2 - w2, 0.0)
400  + emiu * cmplx(8.0 * u2 * cos_w, 2.0 * u * (3.0 * u2 + w2) * xi0);
401  const dcomplex h1 = cmplx(2 * u, 0.0) * e2iu
402  - emiu * cmplx(2.0 * u * cos_w, -(3.0 * u2 - w2) * xi0);
403  const dcomplex h2 = e2iu - emiu * cmplx(cos_w, 3.0 * u * xi0);
404 
405  const double fden = 1.0 / (9.0 * u2 - w2);
406 
407  //- output
408  f0 = h0 * fden;
409  f1 = h1 * fden;
410  f2 = h2 * fden;
411 }
412 
413 
414 //====================================================================
415 void Projection_Stout_SU3::set_uw(double& u, double& w,
416  const Mat_SU_N& iQ2, const Mat_SU_N& iQ3)
417 {
418  const double c0 = -(iQ3.i(0, 0) + iQ3.i(1, 1) + iQ3.i(2, 2)) / 3.0;
419  const double c1 = -0.5 * (iQ2.r(0, 0) + iQ2.r(1, 1) + iQ2.r(2, 2));
420  const double c13r = sqrt(c1 / 3.0);
421  const double c0max = 2.0 * c13r * c13r * c13r;
422 
423  const double theta = acos(c0 / c0max);
424 
425  //- output
426  u = c13r * cos(theta / 3.0);
427  w = sqrt(c1) * sin(theta / 3.0);
428 }
429 
430 
431 //====================================================================
432 double Projection_Stout_SU3::func_xi0(const double w)
433 {
434  if (w == 0.0) {
435  return 1.0;
436  } else {
437  return sin(w) / w;
438  }
439 }
440 
441 
442 //====================================================================
443 double Projection_Stout_SU3::func_xi1(const double w)
444 {
445  if (w < 0.25) {
446  const double w2 = w * w;
447  const static double c0 = -1.0 / 3.0;
448  const static double c1 = 1.0 / 30.0;
449  const static double c2 = -1.0 / 840.0;
450  const static double c3 = 1.0 / 45360.0;
451  const static double c4 = -1.0 / 3991680.0;
452 
453  return c0 + w2 * (c1 + w2 * (c2 + w2 * (c3 + w2 * c4)));
454  } else {
455  return (w * cos(w) - sin(w)) / (w * w * w);
456  }
457 }
458 
459 
460 //====================================================================
462 {
463  // brute force version of exponentiation: for check
464 
465  const static int Nprec = 32;
466 
467  const int Nvol = iQ.nvol();
468  const int Nex = iQ.nex();
469 
470  for (int ex = 0; ex < Nex; ++ex) {
471  for (int site = 0; site < Nvol; ++site) {
472  Mat_SU_N u0(NC);
473  u0.unit();
474 
475  Mat_SU_N u1(NC);
476  u1.unit();
477 
478  Mat_SU_N h1 = iQ.mat(site, ex);
479 
480  for (int iprec = 0; iprec < Nprec; ++iprec) {
481  double exf = 1.0 / (Nprec - iprec);
482  Mat_SU_N u2 = h1 * u1;
483  u2 *= exf;
484  u1 = u2;
485  u1 += u0;
486  }
487 
488  u1.reunit();
489  e_iQ.set_mat(site, ex, u1);
490  }
491  }
492 }
493 
494 
495 //====================================================================
496 //============================================================END=====
BridgeIO vout
Definition: bridgeIO.cpp:503
void force_recursive(Field_G &Xi, Field_G &iTheta, const double alpha, const Field_G &Sigmap, const Field_G &C, const Field_G &U)
determination of fields for force calculation
double i(int c) const
Definition: mat_SU_N.h:115
unsigned long int m_flop
void general(const char *format,...)
Definition: bridgeIO.cpp:197
Mat_SU_N & zero()
Definition: mat_SU_N.h:383
double func_xi1(const double w)
int nvol() const
Definition: field.h:127
void set_uw(double &u, double &w, const Mat_SU_N &iQ2, const Mat_SU_N &iQ3)
Class for parameters.
Definition: parameters.h:46
Mat_SU_N & at()
antihermitian traceless
Definition: mat_SU_N.h:329
double func_xi0(const double w)
void multadd_dn(const Mat_SU_N &u1, const Mat_SU_N &u2)
Definition: mat_SU_N.h:249
Mat_SU_N & reunit()
Definition: mat_SU_N.h:71
void set_i(int c, const double &im)
Definition: mat_SU_N.h:121
int nin() const
Definition: field.h:126
void project(Field_G &U, const double alpha, const Field_G &C, const Field_G &Uorg)
projection U = P[alpha, C, Uorg]
void set_r(int c, const double &re)
Definition: mat_SU_N.h:120
#define NC
SU(N) gauge field.
Definition: field_G.h:38
int nex() const
Definition: field.h:128
void mult_nd(const Mat_SU_N &u1, const Mat_SU_N &u2)
Definition: mat_SU_N.h:197
void set_fj(dcomplex &f0, dcomplex &f1, dcomplex &f2, const double &u, const double &w)
double norm2()
Definition: mat_SU_N.h:155
void set_parameters(const Parameters &param)
void exp_iQ(Field_G &e_iQ, const Field_G &iQ)
void exp_iQ_bf(Field_G &e_iQ, const Field_G &iQ)
static const std::string class_name
static double get_time()
obtain a wall-clock time.
Mat_SU_N & unit()
Definition: mat_SU_N.h:373
double r(int c) const
Definition: mat_SU_N.h:114
void set(int c, double re, const double &im)
Definition: mat_SU_N.h:133
string get_string(const string &key) const
Definition: parameters.cpp:221
void set_mat(const int site, const int mn, const Mat_SU_N &U)
Definition: field_G.h:160
Mat_SU_N mat(const int site, const int mn=0) const
Definition: field_G.h:114
void mult_nn(const Mat_SU_N &u1, const Mat_SU_N &u2)
Definition: mat_SU_N.h:163
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131
Bridge::VerboseLevel m_vl
Definition: projection.h:33
void set_ri(const int cc, const int site, const int mn, const double re, const double im)
Definition: field_G.h:107