Bridge++  Ver. 1.1.x
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fprop_Overlap_5d.cpp
Go to the documentation of this file.
1 
14 #include "fprop_Overlap_5d.h"
15 
16 #ifdef USE_PARAMETERS_FACTORY
17 #include "parameters_factory.h"
18 #endif
19 
20 //- parameter entry
21 namespace {
22  void append_entry(Parameters& param)
23  {
24  param.Register_double("quark_mass", 0.0);
25  param.Register_double("domain_wall_height", 0.0);
26  param.Register_int("number_of_poles", 0);
27  param.Register_double("lower_bound", 0.0);
28  param.Register_double("upper_bound", 0.0);
29  param.Register_int("maximum_number_of_iteration", 0);
30  param.Register_double("convergence_criterion_squared", 0.0);
31  param.Register_int_vector("boundary_condition", std::valarray<int>());
32 
33  param.Register_string("verbose_level", "NULL");
34  }
35 
36 
37 #ifdef USE_PARAMETERS_FACTORY
38  bool init_param = ParametersFactory::Register("Fprop_Overlap_5d", append_entry);
39 #endif
40 }
41 //- end
42 
43 //- parameters class
44 Parameters_Fprop_Overlap_5d::Parameters_Fprop_Overlap_5d() { append_entry(*this); }
45 //- end
46 
47 //====================================================================
48 void Fprop_Overlap_5d::set_parameters(const Parameters& params)
49 {
50  const string str_vlevel = params.get_string("verbose_level");
51 
52  m_vl = vout.set_verbose_level(str_vlevel);
53 
54  //- fetch and check input parameters
55  double mq, M0;
56  int Np;
57  double x_min, x_max;
58  int Niter_ms;
59  double Stop_cond_ms;
60  valarray<int> bc;
61 
62  int err = 0;
63  err += params.fetch_double("quark_mass", mq);
64  err += params.fetch_double("domain_wall_height", M0);
65  err += params.fetch_int("number_of_poles", Np);
66  err += params.fetch_double("lower_bound", x_min);
67  err += params.fetch_double("upper_bound", x_max);
68  err += params.fetch_int("maximum_number_of_iteration", Niter_ms);
69  err += params.fetch_double("convergence_criterion_squared", Stop_cond_ms);
70  err += params.fetch_int_vector("boundary_condition", bc);
71 
72  if (err) {
73  vout.crucial(m_vl, "Fprop_Overlap_5d: fetch error, input parameter not found.\n");
74  abort();
75  }
76 
77 
78  set_parameters(mq, M0, Np, x_min, x_max, Niter_ms, Stop_cond_ms, bc);
79 }
80 
81 
82 //====================================================================
83 void Fprop_Overlap_5d::set_parameters(const Parameters& params_overlap,
84  const Parameters& params_solver)
85 {
86  double mq = params_overlap.get_double("quark_mass");
87  double M0 = params_overlap.get_double("domain_wall_height");
88  int Np = params_overlap.get_int("number_of_poles");
89  double x_min = params_overlap.get_double("lower_bound");
90  double x_max = params_overlap.get_double("upper_bound");
91  int Niter_ms = params_overlap.get_int("maximum_number_of_iteration");
92  double Stop_cond_ms = params_overlap.get_double("convergence_criterion_squared");
93  std::valarray<int> bc = params_overlap.get_int_vector("boundary_condition");
94 
95  set_parameters(mq, M0, Np, x_min, x_max, Niter_ms, Stop_cond_ms, bc);
96  m_params_solver = params_solver;
97 }
98 
99 
100 //====================================================================
101 void Fprop_Overlap_5d::set_parameters(const double mq, const double M0,
102  const int Np, const double x_min, const double x_max,
103  const int Niter_ms, const double Stop_cond_ms,
104  const valarray<int> bc)
105 {
106  int Ndim = CommonParameters::Ndim();
107 
108  //- print input parameters
109  vout.general(m_vl, "Fprop_Overlap_5d paramters:\n");
110  vout.general(m_vl, " mq = %10.6f\n", mq);
111  vout.general(m_vl, " M0 = %10.6f\n", M0);
112  vout.general(m_vl, " Np = %4d\n", Np);
113  vout.general(m_vl, " x_min = %12.8f\n", x_min);
114  vout.general(m_vl, " x_max = %12.6f\n", x_max);
115  vout.general(m_vl, " Niter_ms = %6d\n", Niter_ms);
116  vout.general(m_vl, " Stop_cond_ms = %12.4e\n", Stop_cond_ms);
117  for (int mu = 0; mu < Ndim; ++mu) {
118  vout.general(m_vl, " boundary[%d] = %2d\n", mu, bc[mu]);
119  }
120 
121  //- range check
122  int err = 0;
123  err += ParameterCheck::non_zero(mq);
124  err += ParameterCheck::non_zero(M0);
125  err += ParameterCheck::non_zero(Np);
126  // NB. x_min,x_max == 0 is allowed.
127  err += ParameterCheck::non_zero(Niter_ms);
128  err += ParameterCheck::square_non_zero(Stop_cond_ms);
129 
130  if (err) {
131  vout.crucial(m_vl, "Fprop_Overlap_5d: parameter range check failed.\n");
132  abort();
133  }
134 
135  assert(bc.size() == Ndim);
136 
137  //- store values
138  m_mq = mq;
139  m_M0 = M0;
140  m_Np = Np;
141  m_x_min = x_min;
142  m_x_max = x_max;
143  m_Niter_ms = Niter_ms;
144  m_Stop_cond_ms = Stop_cond_ms;
145 
146  m_boundary.resize(Ndim);
147  assert(bc.size() == Ndim);
148  for (int mu = 0; mu < Ndim; ++mu) {
149  m_boundary[mu] = bc[mu];
150  }
151 
152  m_Nsbt = 0;
153 }
154 
155 
156 //====================================================================
157 void Fprop_Overlap_5d::set_env()
158 {
159  int Nvol = CommonParameters::Nvol();
160  int Ndim = CommonParameters::Ndim();
161 
162  //- setting even-odd envoronement
163  m_Ueo = new Field_G(Nvol, Ndim);
164  m_index_eo.convertField(*m_Ueo, *m_conf);
165 
166  //- setting Wilson fermion operator
167  double kappa = 0.5 / (4.0 - m_M0);
168  m_fopr_w_eo = new Fopr_Wilson_eo("Dirac");
169  m_fopr_w_eo->set_parameters(kappa, m_boundary);
170  // NB. U is converted to Ueo in fopr_eo->set_config
171  // m_fopr_w_eo->set_config(m_Ueo);
172  m_fopr_w_eo->set_config(m_conf);
173 
174  //- setting 5D overlap fermion operator
175  m_fopr_ov5d = new Fopr_Overlap_5d(m_fopr_w_eo);
176  m_fopr_ov5d->set_parameters(m_mq, m_M0, m_Np, m_x_min, m_x_max,
177  m_Niter_ms, m_Stop_cond_ms, m_boundary);
178  m_fopr_ov5d->set_lowmodes(m_Nsbt, m_ev, m_vk);
179 }
180 
181 
182 //====================================================================
183 void Fprop_Overlap_5d::calc_fprop(Field& xq, const Field& b)
184 {
185  set_env();
186 
187  double snorm = 1.0 / b.norm2();
188 
189  Field v1((Field)b);
190  int nconv, nconv1, nconv2;
191  double diff, diff1, diff2;
192 
193  v1 = m_fopr_w_eo->mult_gm5(b);
194  solve_overlap_5D_1(xq, v1, nconv, diff, snorm);
195  //v1 = m_fopr_ov->D(xq);
196  m_fopr_ov->set_mode("D");
197  m_fopr_ov->D(v1, xq);
198 
199  v1 -= b;
200  double vv = v1 * v1;
201  vv = vv * snorm;
202 
203  vout.general(m_vl, "5D overlap solver total:\n");
204  vout.general(m_vl, " Nconv_eo(ov) = %10d\n", nconv);
205  vout.general(m_vl, " Diff_eo(ov) = %16.8e\n", diff);
206  vout.general(m_vl, " Diff(ov) = %16.8e\n", vv);
207 
208  delete m_Ueo;
209  delete m_fopr_w_eo;
210  delete m_fopr_ov5d;
211 }
212 
213 
214 //====================================================================
215 void Fprop_Overlap_5d::calc_H2inv(Field& xq, const Field& b)
216 {
217  set_env();
218 
219  double snorm = 1.0 / b.norm2();
220 
221  Field v1((Field)b);
222  int nconv, nconv1, nconv2;
223  double diff, diff1, diff2;
224 
225  solve_overlap_5D_1(v1, (Field)b, nconv1, diff1, snorm);
226  solve_overlap_5D_1(xq, v1, nconv2, diff2, snorm);
227 
228  nconv = nconv1 + nconv2;
229  diff = diff1 + diff2;
230 
231  // v1 = m_fopr_ov->DdagD(xq);
232  m_fopr_ov->set_mode("DdagD");
233  m_fopr_ov->mult(v1, xq);
234 
235  v1 -= b;
236  double vv = v1 * v1;
237  vv = vv * snorm;
238 
239  vout.general(m_vl, " 5D overlap solver for Q_ov total:\n");
240  vout.general(m_vl, " 5D overlap solver total:\n");
241  vout.general(m_vl, " Nconv_eo(ov) = %10d\n", nconv);
242  vout.general(m_vl, " Diff_eo(ov) = %16.8e\n", diff);
243  vout.general(m_vl, " Diff(ov) = %16.8e\n", vv);
244 
245  delete m_Ueo;
246  delete m_fopr_w_eo;
247  delete m_fopr_ov5d;
248 }
249 
250 
251 //====================================================================
252 void Fprop_Overlap_5d::solve_overlap_5D_1(Field& xq, const Field& b,
253  int& nconv, double& diff, double& snorm)
254 {
255  // #### parameter setup ####
256  int Npl = 2 * m_Np + 1;
257 
258  int Nin = b.nin();
259  int Nvol = b.nvol();
260  int Nvol2 = Nvol / 2;
261 
262  const string str_solver_type = m_params_solver.get_string("solver_type");
263 
264  // #### object setup ####
265  Solver *solver = Solver::New(str_solver_type, m_fopr_ov5d);
266 
267  solver->set_parameters(m_params_solver);
268 
269  Field v(Nin, Nvol2, Npl);
270  Field s(Nin, Nvol2, Npl);
271  Field t(Nin, Nvol2, Npl);
272  Field p(Nin, Nvol2, Npl);
273  Field x(Nin, Nvol2, Npl);
274 
275  Field be(Nin, Nvol2, 1);
276  Field bo(Nin, Nvol2, 1);
277 
278  m_index_eo.convertField(be, b, 0);
279  m_index_eo.convertField(bo, b, 1);
280 
281  v = 0.0;
282  v.setpart_ex(2 * m_Np, be, 0);
283  m_fopr_ov5d->LUprecond(t, v, 0);
284 
285  v = 0.0;
286  v.setpart_ex(2 * m_Np, bo, 0);
287  m_fopr_ov5d->LUprecond(p, v, 1);
288  m_fopr_ov5d->Mopr_5d_eo(v, p, 0);
289  m_fopr_ov5d->LUprecond(p, v, 0);
290 
291  t -= p;
292 
293  s = m_fopr_ov5d->DD_5d_eo(t, -1);
294 
295 
296  solver->solve(x, s, nconv, diff);
297 
298 
299  //- Check
300  v = 0.0;
301  v.setpart_ex(2 * m_Np, be, 0);
302  m_fopr_ov5d->LUprecond(s, v, 0);
303 
304  v = 0.0;
305  v.setpart_ex(2 * m_Np, bo, 0);
306  m_fopr_ov5d->LUprecond(p, v, 1);
307  m_fopr_ov5d->Mopr_5d_eo(t, p, 0);
308  m_fopr_ov5d->LUprecond(p, t, 0);
309 
310  s -= p;
311 
312  p = m_fopr_ov5d->DD_5d_eo(x, 1);
313  p -= s;
314  diff = p * p;
315  diff *= snorm;
316 
317  //- Solution
318  Field xqe(Nin, Nvol2, 1);
319  Field xqo(Nin, Nvol2, 1);
320  Field xt(Nin, Nvol2, 1);
321 
322  xqe.setpart_ex(0, x, 2 * m_Np);
323 
324  v = 0.0;
325  v.setpart_ex(2 * m_Np, bo, 0);
326  m_fopr_ov5d->LUprecond(p, v, 1);
327  xqo.setpart_ex(0, p, 2 * m_Np);
328 
329  m_fopr_ov5d->Mopr_5d_eo(p, x, 1);
330  m_fopr_ov5d->LUprecond(t, p, 1);
331  xt.setpart_ex(0, t, 2 * m_Np);
332  xqo -= xt;
333 
334  m_index_eo.reverseField(xq, xqe, 0);
335  m_index_eo.reverseField(xq, xqo, 1);
336 
337  vout.general(m_vl, " solver_eo1(ov):\n");
338  vout.general(m_vl, " Nconv_eo1(ov) = %10d\n", nconv);
339  vout.general(m_vl, " Diff_eo1(ov) = %16.8e\n", diff);
340 
341 
342  // #### tidy up ####
343  delete solver;
344 }
345 
346 
347 //====================================================================
348 //============================================================END=====