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