Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_Overlap_5d.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Overlap_5d.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 using std::valarray;
21 
22 //- parameter entries
23 namespace {
24  void append_entry(Parameters& param)
25  {
26  param.Register_double("quark_mass", 0.0);
27  param.Register_double("domain_wall_height", 0.0);
28  param.Register_int("number_of_poles", 0);
29  param.Register_double("lower_bound", 0.0);
30  param.Register_double("upper_bound", 0.0);
31  param.Register_int("maximum_number_of_iteration", 0);
32  param.Register_double("convergence_criterion_squared", 0.0);
33  param.Register_int_vector("boundary_condition", std::valarray<int>());
34 
35  param.Register_string("verbose_level", "NULL");
36  }
37 
38 
39 #ifdef USE_PARAMETERS_FACTORY
40  bool init_param = ParametersFactory::Register("Fopr.Overlap_5d", append_entry);
41 #endif
42 }
43 //- end
44 
45 //- parameters class
46 Parameters_Fopr_Overlap_5d::Parameters_Fopr_Overlap_5d() { append_entry(*this); }
47 //- end
48 
49 //====================================================================
50 void Fopr_Overlap_5d::set_parameters(const Parameters& params)
51 {
52  const string str_vlevel = params.get_string("verbose_level");
53 
54  m_vl = vout.set_verbose_level(str_vlevel);
55 
56  //- fetch and check input parameters
57  double mq, M0;
58  int Np;
59  double x_min, x_max;
60  int Niter_ms;
61  double Stop_cond_ms;
62  valarray<int> bc;
63 
64  int err = 0;
65  err += params.fetch_double("quark_mass", mq);
66  err += params.fetch_double("domain_wall_height", M0);
67  err += params.fetch_int("number_of_poles", Np);
68  err += params.fetch_double("lower_bound", x_min);
69  err += params.fetch_double("upper_bound", x_max);
70  err += params.fetch_int("maximum_number_of_iteration", Niter_ms);
71  err += params.fetch_double("convergence_criterion_squared", Stop_cond_ms);
72  err += params.fetch_int_vector("boundary_condition", bc);
73 
74  if (err) {
75  vout.crucial(m_vl, "Fopr_Overlap_5d: fetch error, input parameter not found.\n");
76  abort();
77  }
78 
79  set_parameters(mq, M0, Np, x_min, x_max, Niter_ms, Stop_cond_ms, bc);
80 }
81 
82 
83 //====================================================================
84 void Fopr_Overlap_5d::set_parameters(const double mq, const double M0,
85  const int Np, const double x_min, const double x_max,
86  const int Niter_ms, const double Stop_cond_ms,
87  const std::valarray<int> bc)
88 {
89  int Ndim = CommonParameters::Ndim();
90 
91  //- print input parameters
92  vout.general(m_vl, "Parameters of 5D overlap fermion operator:\n");
93  vout.general(m_vl, " mq = %8.4f\n", mq);
94  vout.general(m_vl, " M0 = %8.4f\n", M0);
95  vout.general(m_vl, " Np = %4d\n", Np);
96  vout.general(m_vl, " x_min = %10.6f\n", x_min);
97  vout.general(m_vl, " x_max = %10.6f\n", x_max);
98  vout.general(m_vl, " Niter_ms = %8d\n", Niter_ms);
99  vout.general(m_vl, " Stop_cond_ms = %10.4e\n", Stop_cond_ms);
100  for (int mu = 0; mu < Ndim; ++mu) {
101  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
102  }
103 
104  //- range check
105  int err = 0;
106  err += ParameterCheck::non_zero(mq);
107  err += ParameterCheck::non_zero(M0);
108  err += ParameterCheck::non_zero(Np);
109  // NB. x_min,x_max == 0 is allowed.
110  err += ParameterCheck::non_zero(Niter_ms);
111  err += ParameterCheck::square_non_zero(Stop_cond_ms);
112 
113  if (err) {
114  vout.crucial(m_vl, "Fopr_Overlap: parameter range check failed.\n");
115  abort();
116  }
117 
118  assert(bc.size() == Ndim);
119 
120  //- store values
121  m_mq = mq;
122  m_M0 = M0;
123 
124  m_Np = Np;
125  m_x_min = x_min;
126  m_x_max = x_max;
127 
128  m_Niter_ms = Niter_ms;
129  m_Stop_cond_ms = Stop_cond_ms;
130 
131  m_boundary.resize(Ndim);
132  for (int mu = 0; mu < Ndim; ++mu) {
133  m_boundary[mu] = bc[mu];
134  }
135 
136  //- propagate parameters
137  //- Zolotarev coefficients and shift values
138  m_sigma.resize(m_Np);
139  m_cl.resize(2 * m_Np);
140  m_bl.resize(m_Np);
141 
142  //- Zolotarev coefficient defined
143  double bmax = m_x_max / m_x_min;
144  Math_Sign_Zolotarev sign_func(m_Np, bmax);
145  sign_func.get_sign_parameters(m_cl, m_bl);
146 
147  for (int i = 0; i < m_Np; i++) {
148  m_sigma[i] = m_cl[2 * i] * m_x_min * m_x_min;
149  }
150 
151  for (int i = 0; i < m_Np; i++) {
152  vout.general(m_vl, " %3d %12.4e %12.4e %12.4e\n",
153  i, m_cl[i], m_cl[i + m_Np], m_bl[i]);
154  }
155 
156  m_p_sqrt.resize(m_Np);
157  m_q_sqrt.resize(m_Np);
158 
159  double b_sum = 0.0;
160  for (int ip = 0; ip < m_Np; ++ip) {
161  int ik = m_Np - ip - 1;
162  double p_tmp = (m_M0 - 0.5 * m_mq) * m_bl[ip] *
163  (m_cl[2 * m_Np - 1] - m_cl[2 * ip]) * m_x_min;
164  m_p_sqrt[ik] = sqrt(p_tmp);
165  m_q_sqrt[ik] = sqrt(m_cl[2 * ip] * m_x_min * m_x_min);
166  b_sum = b_sum + m_bl[ip];
167  }
168 
169  m_p0_parameter = (m_M0 - 0.5 * m_mq) * b_sum / m_x_min;
170  m_R_parameter = m_M0 + 0.5 * m_mq;
171  m_h = 1.0;
172 
173  m_rl.resize(m_Np);
174  m_sl.resize(m_Np);
175  for (int ik = 0; ik < m_Np; ++ik) {
176  m_rl[ik] = -m_q_sqrt[ik];
177  m_sl[ik] = -m_p_sqrt[ik] / (1.0 + m_q_sqrt[ik] * m_q_sqrt[ik]);
178  }
179 
180  m_u0 = 0.0;
181  for (int ik = 0; ik < m_Np; ++ik) {
182  m_u0 = m_u0 - m_sl[ik] * m_p_sqrt[ik];
183  }
184 }
185 
186 
187 //====================================================================
188 void Fopr_Overlap_5d::set_lowmodes(int Nsbt, valarray<double> *ev,
189  valarray<Field> *vk)
190 {
191  m_Nsbt = Nsbt;
192 
193  vout.general(m_vl, " Nsbt = %d", Nsbt);
194 
195  if (m_Nsbt > 0) {
196  m_ev.resize(m_Nsbt);
197  m_vk.resize(2 * m_Nsbt);
198 
199  int Nin = (*vk)[0].nin();
200  int Nvol = (*vk)[0].nvol();
201  int Nvol2 = Nvol / 2;
202 
203  for (int k = 0; k < 2 * m_Nsbt; ++k) {
204  m_vk[k].reset(Nin, Nvol2, 1);
205  }
206 
207  Field vt(Nin, Nvol, 1);
208  Field vt_eo(Nin, Nvol2, 1);
209 
210  for (int k = 0; k < m_Nsbt; ++k) {
211  m_ev[k] = (*ev)[k];
212  vt = (*vk)[k];
213  m_index_eo.convertField(m_vk[k], vt, 0);
214  m_index_eo.convertField(m_vk[k + m_Nsbt], vt, 1);
215  }
216  }
217 
218  //- setting parameter for low-mode subtraction
219  m_prf.resize(m_Nsbt);
220  for (int k = 0; k < m_Nsbt; ++k) {
221  double sign_ev = m_ev[k] / fabs(m_ev[k]);
222  m_prf[k] = (m_M0 - 0.5 * m_mq) * sign_ev
223  - m_p0_parameter * m_ev[k];
224  }
225 
226  int Nsbt2 = 2 * m_Nsbt;
227 
228  m_u0c_e.resize(Nsbt2 * Nsbt2);
229  m_u0cinv_e.resize(Nsbt2 * Nsbt2);
230  m_u0c_o.resize(Nsbt2 * Nsbt2);
231  m_u0cinv_o.resize(Nsbt2 * Nsbt2);
232 
233  if (m_Nsbt > 0) {
234  Calc_Coeff_u0inv();
235  }
236 }
237 
238 
239 //====================================================================
240 const Field Fopr_Overlap_5d::DdagD_eo(const Field& f)
241 {
242  Field v(f);
243  Field w(f);
244 
245  v = DD_5d_eo(f, 1);
246  w = DD_5d_eo(v, -1);
247 
248  return w;
249 }
250 
251 
252 //====================================================================
253 const Field Fopr_Overlap_5d::DD_5d_eo(const Field& w, const int jd)
254 {
255  int Nin = w.nin();
256  int Nvol = w.nvol();
257  int Nex = w.nex();
258  Field t(Nin, Nvol, Nex), v(Nin, Nvol, Nex);
259 
260  if (jd == 1) {
261  Mopr_5d_eo(t, w, 1);
262  LUprecond(v, t, 1);
263  Mopr_5d_eo(t, v, 0);
264  LUprecond(v, t, 0);
265  v = -v + w;
266  } else if (jd == -1) {
267  LUprecond(t, w, 0);
268  Mopr_5d_eo(v, t, 1);
269  LUprecond(t, v, 1);
270  Mopr_5d_eo(v, t, 0);
271  v = -v + w;
272  } else {
273  vout.crucial(m_vl, "Fopr_overlap_5d: illegal jd.\n");
274  abort();
275  }
276 
277  return v;
278 }
279 
280 
281 //====================================================================
282 void Fopr_Overlap_5d::Mopr_5d_eo(Field& v, const Field& w,
283  const int ieo)
284 {
285  // ieo = 1: M_eo, 2: M_oe
286 
287  int ieo1 = ieo;
288  int ieo2 = 1 - ieo;
289 
290  int Nin = w.nin();
291  int Nvol2 = w.nvol();
292  int Nex = w.nex();
293 
294  Field z(Nin, Nvol2, 1);
295  Field z1(Nin, Nvol2, 1);
296  Field z2(Nin, Nvol2, 1);
297 
298  Field w1(Nin, Nvol2, 1);
299  Field v1(Nin, Nvol2, 1);
300 
301  w1.setpart_ex(0, w, 2 * m_Np);
302  Proj_H_eo(ieo1, ieo2, z1, w1);
303 
304  z2 = 0.0;
305 
306  for (int j = 0; j < m_Np; ++j) {
307  w1.setpart_ex(0, w, 2 * j);
308  v1 = (Field)m_fopr_w->Meo_gm5(w1, ieo);
309  v.setpart_ex(2 * j, v1, 0);
310 
311  w1.setpart_ex(0, w, 2 * j + 1);
312  z2 += m_p_sqrt[j] * w1;
313 
314  v1 = (Field)m_fopr_w->Meo_gm5(w1, ieo);
315  v1 *= -1.0;
316  v1 += m_p_sqrt[j] * z1;
317  v.setpart_ex(2 * j + 1, v1, 0);
318  }
319 
320  Proj_H_eo(ieo1, ieo2, v1, z2);
321  w1.setpart_ex(0, w, 2 * m_Np);
322  Proj_L_mult_eo(ieo1, ieo2, z1, w1);
323  v1 += z1;
324  z1 = (Field)m_fopr_w->Meo_gm5(w1, ieo);
325  v1 += m_p0_parameter * z1;
326  v.setpart_ex(2 * m_Np, v1, 0);
327 }
328 
329 
330 //====================================================================
331 void Fopr_Overlap_5d::LUprecond(Field& v, const Field& w,
332  const int ieo)
333 {
334  // ieo=0: 1/M_ee, 1: 1/M_oo
335 
336  Field t(w);
337 
338  int Nin = w.nin();
339  int Nvol2 = w.nvol();
340  int Nex = w.nex();
341 
342  Field vx(Nin, Nvol2, 1);
343  Field vy(Nin, Nvol2, 1);
344  Field vj(Nin, Nvol2, 1);
345  Field v_tmp(Nin, Nvol2, 1);
346 
347  // --- L^-1 ---
348  for (int ip = 0; ip < m_Np; ++ip) {
349  int jx = 2 * ip;
350  int jy = 2 * ip + 1;
351 
352  vx.setpart_ex(0, w, jx);
353  vy.setpart_ex(0, w, jy);
354 
355  v_tmp = m_fopr_w->mult_gm5(vx);
356  vy += -m_rl[ip] * v_tmp;
357 
358  t.setpart_ex(jx, vx, 0);
359  t.setpart_ex(jy, vy, 0);
360  }
361 
362  int j = 2 * m_Np;
363  vj = 0.0;
364  for (int ip = 0; ip < m_Np; ++ip) {
365  int jy = 2 * ip + 1;
366  vy.setpart_ex(0, t, jy);
367  vj += -m_sl[ip] * vy;
368  }
369  v_tmp = m_fopr_w->mult_gm5(vj);
370  Proj_H_eo(ieo, ieo, vj, v_tmp);
371  v_tmp.setpart_ex(0, w, j);
372  vj += v_tmp;
373  t.setpart_ex(j, vj, 0);
374 
375  // --- U^-1 ---
376  mult_u0inv(v_tmp, vj, ieo);
377  v.setpart_ex(j, v_tmp, 0);
378 
379  Proj_H_eo(ieo, ieo, vj, v_tmp);
380  t.setpart_ex(j, vj, 0);
381 
382  for (int ip = 0; ip < m_Np; ++ip) {
383  int jx = 2 * ip;
384  int jy = 2 * ip + 1;
385 
386  v_tmp.setpart_ex(0, t, jy);
387  vj.setpart_ex(0, t, j);
388  v_tmp += -m_p_sqrt[ip] * vj;
389  vy = m_fopr_w->mult_gm5(v_tmp);
390  vy *= -1.0 / (1.0 + m_q_sqrt[ip] * m_q_sqrt[ip]);
391  v.setpart_ex(jy, vy, 0);
392 
393  v_tmp.setpart_ex(0, t, jx);
394  v_tmp += m_q_sqrt[ip] * vy;
395  vx = m_fopr_w->mult_gm5(v_tmp);
396  v.setpart_ex(jx, vx, 0);
397  }
398 }
399 
400 
401 //====================================================================
402 void Fopr_Overlap_5d::Proj_H_eo(const int ieo1, const int ieo2,
403  Field& v, const Field& w)
404 {
405  if (ieo1 == ieo2) {
406  v = w;
407  } else {
408  v = 0.0;
409  }
410 
411  int j1, j2;
412  double prd_r, prd_i, vr, vi;
413 
414  for (int k = 0; k < m_Nsbt; ++k) {
415  j1 = k + ieo1 * m_Nsbt;
416  j2 = k + ieo2 * m_Nsbt;
417  innerprd_c(prd_r, prd_i, m_vk[j2], w);
418  add_c(v, m_vk[j1], -prd_r, -prd_i);
419  }
420 }
421 
422 
423 //====================================================================
424 void Fopr_Overlap_5d::Proj_L_mult_eo(const int ieo1, const int ieo2,
425  Field& v, const Field& w)
426 {
427  int j1, j2;
428  double prd_r, prd_i, v_r, v_i;
429 
430  v = 0.0;
431  for (int k = 0; k < m_Nsbt; ++k) {
432  j1 = k + ieo1 * m_Nsbt;
433  j2 = k + ieo2 * m_Nsbt;
434  innerprd_c(prd_r, prd_i, m_vk[j2], w);
435  prd_r *= m_prf[k];
436  prd_i *= m_prf[k];
437  add_c(v, m_vk[j1], prd_r, prd_i);
438  }
439 }
440 
441 
442 //====================================================================
443 void Fopr_Overlap_5d::Calc_Coeff_u0inv()
444 {
445  int Nsbt = m_Nsbt;
446  int Nsbt2 = 2 * m_Nsbt;
447 
448  valarray<dcomplex> c(Nsbt2 * Nsbt2);
449  valarray<dcomplex> cinv(Nsbt2 * Nsbt2);
450  valarray<dcomplex> vprd(Nsbt2 * Nsbt2);
451  valarray<dcomplex> W(Nsbt2 * Nsbt2);
452 
453  valarray<dcomplex> c_src(Nsbt2);
454  valarray<dcomplex> c_x(Nsbt2);
455 
456  int Nin = m_vk[0].nin();
457  int Nvol = m_vk[0].nvol();
458  Field w(Nin, Nvol, 1);
459  Field v1(Nin, Nvol, 1);
460  Field v2(Nin, Nvol, 1);
461 
462  double u0_a = m_R_parameter + m_p0_parameter + m_u0;
463  double u0inv_a = 1.0 / u0_a;
464 
465  double a_r, a_i;
466 
467  for (int ieo = 0; ieo < 2; ++ieo) {
468  //- Determine inner product of eigenvectors
469  for (int i = 0; i < Nsbt; ++i) {
470  for (int j = 0; j < i + 1; ++j) {
471  int i2 = i + ieo * Nsbt;
472  int j2 = j + ieo * Nsbt;
473  v1 = m_vk[i2];
474  v2 = m_vk[j2];
475  innerprd_c(a_r, a_i, v1, v2);
476  vprd[i + j * Nsbt2] = cmplx(a_r, a_i);
477  vprd[i + Nsbt + (j + Nsbt) * Nsbt2] = cmplx(a_r, a_i);
478  vprd[j + i * Nsbt2] = cmplx(a_r, -a_i);
479  vprd[j + Nsbt + (i + Nsbt) * Nsbt2] = cmplx(a_r, -a_i);
480  }
481  }
482 
483  for (int i = 0; i < Nsbt; ++i) {
484  for (int j = 0; j < i + 1; ++j) {
485  int i2 = i + ieo * Nsbt;
486  int j2 = j + ieo * Nsbt;
487  v1 = m_vk[i2];
488  w = m_vk[j2];
489  v2 = m_fopr_w->mult_gm5(w);
490  innerprd_c(a_r, a_i, v1, v2);
491  vprd[i + (j + Nsbt) * Nsbt2] = cmplx(a_r, a_i);
492  vprd[i + Nsbt + j * Nsbt2] = cmplx(a_r, a_i);
493  vprd[j + (i + Nsbt) * Nsbt2] = cmplx(a_r, -a_i);
494  vprd[j + Nsbt + i * Nsbt2] = cmplx(a_r, -a_i);
495  }
496  }
497 
498  //- definition of matrix c(i,j)
499  c = cmplx(0.0, 0.0);
500  for (int i = 0; i < Nsbt2; ++i) {
501  c[i + i * Nsbt2] = cmplx(-m_u0, 0.0);
502  }
503 
504  for (int i = 0; i < Nsbt; ++i) {
505  double prf = (m_M0 - 0.5 * m_mq) * (m_ev[i] / fabs(m_ev[i]))
506  - m_p0_parameter * m_ev[i];
507  c[i + (i + Nsbt) * Nsbt2] = cmplx(prf, 0.0);
508  }
509 
510  for (int i = 0; i < Nsbt; ++i) {
511  for (int j = Nsbt; j < Nsbt2; ++j) {
512  c[i + j * Nsbt2] += cmplx(m_u0, 0.0) * vprd[i + j * Nsbt2];
513  }
514  }
515 
516  //- Definition of matrix W(i,j)
517  W = cmplx(0.0, 0.0);
518  for (int i = 0; i < Nsbt2; ++i) {
519  W[i + i * Nsbt2] = cmplx(u0_a, 0.0);
520  }
521 
522  for (int i = 0; i < Nsbt2; ++i) {
523  for (int j = 0; j < Nsbt2; ++j) {
524  for (int k = 0; k < Nsbt2; ++k) {
525  W[i + j * Nsbt2] += c[i + k * Nsbt2] * vprd[k + j * Nsbt2];
526  }
527  }
528  }
529 
530  //- Solve cinv
531  for (int i = 0; i < Nsbt2; ++i) {
532  for (int j = 0; j < Nsbt2; ++j) {
533  c_src[j] = cmplx(-u0inv_a, 0.0) * c[j + i * Nsbt2];
534  }
535 
536  Solv_Coeff_u0inv(c_x, W, c_src);
537 
538  for (int j = 0; j < Nsbt2; ++j) {
539  cinv[j + i * Nsbt2] = c_x[j];
540  }
541  }
542 
543  if (ieo == 0) {
544  for (int i = 0; i < Nsbt2 * Nsbt2; ++i) {
545  m_u0c_e[i] = c[i];
546  m_u0cinv_e[i] = cinv[i];
547  }
548  } else {
549  for (int i = 0; i < Nsbt2 * Nsbt2; ++i) {
550  m_u0c_o[i] = c[i];
551  m_u0cinv_o[i] = cinv[i];
552  }
553  }
554  }
555 }
556 
557 
558 //====================================================================
559 void Fopr_Overlap_5d::innerprd_c(double& prd_r, double& prd_i,
560  const Field& v, const Field& w)
561 {
562  int size = w.size();
563 
564  assert(v.size() == size);
565 
566  prd_r = 0.0;
567  prd_i = 0.0;
568 
569  for (int i = 0; i < size; i += 2) {
570  prd_r += v.cmp(i) * w.cmp(i) + v.cmp(i + 1) * w.cmp(i + 1);
571  prd_i += v.cmp(i) * w.cmp(i + 1) - v.cmp(i + 1) * w.cmp(i);
572  }
573 
574  prd_r = Communicator::reduce_sum(prd_r);
575  prd_i = Communicator::reduce_sum(prd_i);
576 }
577 
578 
579 //====================================================================
580 void Fopr_Overlap_5d::add_c(Field& v, const Field& w,
581  const double a_r, const double a_i)
582 {
583  int size = w.size();
584 
585  assert(v.size() == size);
586 
587  double v_r, v_i;
588  for (int i = 0; i < size; i += 2) {
589  v_r = a_r * w.cmp(i) - a_i * w.cmp(i + 1);
590  v_i = a_r * w.cmp(i + 1) + a_i * w.cmp(i);
591  v.add(i, v_r);
592  v.add(i + 1, v_i);
593  }
594 }
595 
596 
597 //====================================================================
598 void Fopr_Overlap_5d::mult_u0inv(Field& v1, const Field& w1,
599  const int ieo)
600 {
601  if (m_Nsbt == 0) {
602  v1 = m_fopr_w->mult_gm5(w1);
603  v1 *= 1.0 / (m_R_parameter + m_p0_parameter + m_u0);
604  } else {
605  int Nin = w1.nin();
606  int Nvol2 = w1.nvol();
607  Field vt(Nin, Nvol2, 1);
608 
609  int Nsbt = m_Nsbt;
610  int Nsbt2 = 2 * Nsbt;
611  valarray<dcomplex> prd_vb(Nsbt2), coeff(Nsbt2);
612  valarray<dcomplex> u0c(Nsbt2 * Nsbt2), u0cinv(Nsbt2 * Nsbt2);
613 
614  if (ieo == 0) {
615  u0c = m_u0c_e;
616  u0cinv = m_u0cinv_e;
617  } else {
618  u0c = m_u0c_o;
619  u0cinv = m_u0cinv_o;
620  }
621 
622  double a_r, a_i;
623  for (int k = 0; k < Nsbt; ++k) {
624  innerprd_c(a_r, a_i, m_vk[k + ieo * Nsbt], w1);
625  prd_vb[k] = cmplx(a_r, a_i);
626  }
627  vt = m_fopr_w->mult_gm5(w1);
628  for (int k = 0; k < Nsbt; ++k) {
629  innerprd_c(a_r, a_i, m_vk[k + ieo * Nsbt], vt);
630  prd_vb[k + Nsbt] = cmplx(a_r, a_i);
631  }
632 
633  for (int i = 0; i < Nsbt2; ++i) {
634  coeff[i] = cmplx(0.0, 0.0);
635  for (int j = 0; j < Nsbt2; ++j) {
636  coeff[i] += u0cinv[i + j * Nsbt2] * prd_vb[j];
637  }
638  }
639 
640  double u0inv_a = 1.0 / (m_R_parameter + m_p0_parameter + m_u0);
641 
642  vt = u0inv_a * w1;
643 
644  for (int k = 0; k < Nsbt; ++k) {
645  add_c(vt, m_vk[k + ieo * Nsbt], real(coeff[k]), imag(coeff[k]));
646  }
647 
648  v1 = m_fopr_w->mult_gm5(vt);
649 
650  for (int k = 0; k < Nsbt; ++k) {
651  add_c(v1, m_vk[k + ieo * Nsbt],
652  real(coeff[k + Nsbt]), imag(coeff[k + Nsbt]));
653  }
654  }
655 }
656 
657 
658 //====================================================================
659 void Fopr_Overlap_5d::Solv_Coeff_u0inv(valarray<dcomplex>& c_x,
660  const valarray<dcomplex>& W, const valarray<dcomplex>& c_src)
661 {
662  //- This is an implementation of CG solver.
663  int Nsbt = m_Nsbt;
664  int Nsbt2 = 2 * m_Nsbt;
665 
666  assert(c_x.size() == Nsbt2);
667  assert(c_src.size() == Nsbt2);
668  assert(W.size() == (Nsbt2 * Nsbt2));
669 
670  valarray<dcomplex> x(Nsbt2);
671  valarray<dcomplex> r(Nsbt2);
672  valarray<dcomplex> p(Nsbt2);
673  valarray<dcomplex> s(Nsbt2);
674  valarray<dcomplex> vt(Nsbt2);
675 
676  int Niter = 100;
677  double Encg = 1.e-32;
678 
679  double ww = norm_c(c_src);
680  double snorm = 1.0 / ww;
681 
682  double rr, rr0;
683  int nconv = -1; // superficial initialization
684 
685  s = cmplx(0.0, 0.0);
686  for (int i = 0; i < Nsbt2; ++i) {
687  for (int j = 0; j < Nsbt2; ++j) {
688  s[i] += conj(W[j + i * Nsbt2]) * c_src[j];
689  }
690  }
691  x = s;
692  r = s;
693  mult_WdagW(s, W, x);
694 
695  r -= s;
696  p = r;
697  rr = norm_c(r);
698 
699  vout.detailed(m_vl, " init %16.8e\n", rr * snorm);
700 
701  if (rr * snorm < Encg) goto converged;
702 
703 
704  for (int iter = 0; iter < Niter; ++iter) {
705  mult_WdagW(s, W, p);
706 
707  double pap = innerprd_c(p, s);
708  double alpha = rr / pap;
709 
710  x += cmplx(alpha, 0.0) * p;
711  r -= cmplx(alpha, 0.0) * s;
712  rr0 = rr;
713  rr = norm_c(r);
714 
715  vout.detailed(m_vl, " %4d %16.8e\n", iter, rr * snorm);
716 
717  if (rr * snorm < Encg) {
718  nconv = iter;
719  goto converged;
720  }
721 
722  double beta = rr / rr0;
723 
724  p *= cmplx(beta, 0.0);
725  p += r;
726  }
727 
728  nconv = -1;
729 
730  vout.crucial(m_vl, "Fopr_Overlap_5d: NOT CONVERGED.\n");
731  abort();
732 
733 
734 converged:
735  vout.detailed(m_vl, " converged\n");
736 
737  s = cmplx(0.0, 0.0);
738  for (int i = 0; i < Nsbt2; ++i) {
739  for (int j = 0; j < Nsbt2; ++j) {
740  s[i] += W[i + j * Nsbt2] * x[j];
741  }
742  }
743  s -= c_src;
744 
745  double diff = norm_c(s);
746  diff *= snorm;
747 
748  c_x = x;
749 
750  vout.general(m_vl, " u0 solver: Nconv = %4d, diff = %12.4e\n", nconv, diff);
751 }
752 
753 
754 //====================================================================
755 void Fopr_Overlap_5d::mult_WdagW(valarray<dcomplex>& v2,
756  const valarray<dcomplex>& W,
757  const valarray<dcomplex>& v1)
758 {
759  int size = v1.size();
760 
761  assert(v2.size() == size);
762  assert(W.size() == (size * size));
763 
764  valarray<dcomplex> vt(size);
765 
766  vt = cmplx(0.0, 0.0);
767  for (int i = 0; i < size; ++i) {
768  for (int j = 0; j < size; ++j) {
769  vt[i] += W[i + j * size] * v1[j];
770  }
771  }
772 
773  v2 = cmplx(0.0, 0.0);
774  for (int i = 0; i < size; ++i) {
775  for (int j = 0; j < size; ++j) {
776  v2[i] += conj(W[j + i * size]) * vt[j];
777  }
778  }
779 }
780 
781 
782 //====================================================================
783 double Fopr_Overlap_5d::norm_c(const valarray<dcomplex>& v)
784 {
785  int size = v.size();
786 
787  double vv = 0.0;
788 
789  for (int i = 0; i < size; ++i) {
790  vv += real(v[i]) * real(v[i]) + imag(v[i]) * imag(v[i]);
791  }
792 
793  return vv;
794 }
795 
796 
797 //====================================================================
798 double Fopr_Overlap_5d::innerprd_c(const valarray<dcomplex>& v,
799  const valarray<dcomplex>& w)
800 {
801  int size = v.size();
802 
803  double vw = 0.0;
804 
805  for (int i = 0; i < size; ++i) {
806  vw += real(v[i]) * real(w[i]) + imag(v[i]) * imag(w[i]);
807  }
808 
809  return vw;
810 }
811 
812 
813 //====================================================================
814 //============================================================END=====