Bridge++  Ver. 2.0.2
afopr_Overlap-tmpl.h
Go to the documentation of this file.
1 
14 #include "Fopr/afopr_Overlap.h"
15 
16 #include "ResourceManager/threadManager.h"
17 
18 #ifdef USE_FACTORY_AUTOREGISTER
19 namespace {
20  bool init = AFopr_Overlap<AFIELD>::register_factory();
21 }
22 #endif
23 
24 template<typename AFIELD>
25 const std::string AFopr_Overlap<AFIELD>::class_name = "Fopr_Overlap";
26 
27 //====================================================================
28 template<typename AFIELD>
29 void AFopr_Overlap<AFIELD>::init(const Parameters& params)
30 {
31  ThreadManager::assert_single_thread(class_name);
32 
33  m_vl = CommonParameters::Vlevel();
34 
35  vout.general(m_vl, "%s: construction\n", class_name.c_str());
36  vout.increase_indent();
37 
38  m_is_initial_step = true;
39 
40  std::string kernel_type;
41  int err = params.fetch_string("kernel_type", kernel_type);
42  if (err > 0) {
43  vout.crucial(m_vl, "Error at %s: kernel_type is not specified.\n",
44  class_name.c_str());
45  exit(EXIT_FAILURE);
46  }
47  m_kernel_type = kernel_type;
48 
49  double M0;
50  params.fetch_double("domain_wall_height", M0);
51 
52  Parameters params_kernel = params;
53  double kappa = 1.0 / (8.0 - 2.0 * M0);
54  params_kernel.set_double("hopping_parameter", kappa);
55 
56  m_fopr_w = AFOPR::New(kernel_type, params_kernel);
57  m_kernel_created = true;
58 
59  m_Nin = m_fopr_w->field_nin();
60  m_Nvol = m_fopr_w->field_nvol();
61  m_Nex = m_fopr_w->field_nex();
62 
63  m_sign = new AFopr_Sign<AFIELD>(m_fopr_w, params);
64 
65  const int Ndim = CommonParameters::Ndim();
66  m_boundary.resize(Ndim);
67 
68  m_w0.reset(m_Nin, m_Nvol, m_Nex);
69  m_w1.reset(m_Nin, m_Nvol, m_Nex);
70  m_w2.reset(m_Nin, m_Nvol, m_Nex);
71 
72  set_parameters(params);
73 
74  m_is_initial_step = false;
75 
76  vout.decrease_indent();
77  vout.general(m_vl, "%s: construction finished.\n",
78  class_name.c_str());
79 }
80 
81 
82 //====================================================================
83 template<typename AFIELD>
84 void AFopr_Overlap<AFIELD>::init(AFOPR *fopr, const Parameters& params)
85 {
86  ThreadManager::assert_single_thread(class_name);
87 
88  m_vl = CommonParameters::Vlevel();
89 
90  vout.general(m_vl, "%s: construction\n", class_name.c_str());
91  vout.increase_indent();
92 
93  m_is_initial_step = true;
94 
95  m_fopr_w = fopr;
96  m_kernel_created = false;
97 
98  m_kernel_type = "unknown";
99  m_repr = "unknown";
100 
101  m_Nin = m_fopr_w->field_nin();
102  m_Nvol = m_fopr_w->field_nvol();
103  m_Nex = m_fopr_w->field_nex();
104 
105  m_sign = new AFopr_Sign<AFIELD>(m_fopr_w, params);
106 
107  const int Ndim = CommonParameters::Ndim();
108  m_boundary.resize(Ndim);
109 
110  m_w0.reset(m_Nin, m_Nvol, m_Nex);
111  m_w1.reset(m_Nin, m_Nvol, m_Nex);
112  m_w2.reset(m_Nin, m_Nvol, m_Nex);
113 
114  set_parameters(params);
115 
116  m_is_initial_step = false;
117 
118  vout.decrease_indent();
119  vout.general(m_vl, "%s: construction finished.\n",
120  class_name.c_str());
121 }
122 
123 
124 //====================================================================
125 template<typename AFIELD>
126 void AFopr_Overlap<AFIELD>::init()
127 {
128  ThreadManager::assert_single_thread(class_name);
129 
130  m_vl = CommonParameters::Vlevel();
131 
132  vout.general(m_vl, "%s: construction (obsolete)\n",
133  class_name.c_str());
134  vout.increase_indent();
135 
136  m_is_initial_step = true;
137  m_kernel_created = false;
138 
139  m_kernel_type = "unknown";
140  m_repr = "unknown";
141 
142  m_Nin = m_fopr_w->field_nin();
143  m_Nvol = m_fopr_w->field_nvol();
144  m_Nex = m_fopr_w->field_nex();
145 
146  m_sign = new AFopr_Sign<AFIELD>(m_fopr_w);
147 
148  const int Ndim = CommonParameters::Ndim();
149  m_boundary.resize(Ndim);
150 
151  m_w0.reset(m_Nin, m_Nvol, m_Nex);
152  m_w1.reset(m_Nin, m_Nvol, m_Nex);
153  m_w2.reset(m_Nin, m_Nvol, m_Nex);
154 
155  m_is_initial_step = false;
156 
157  vout.decrease_indent();
158  vout.general(m_vl, "%s: construction finished.\n",
159  class_name.c_str());
160 }
161 
162 
163 //====================================================================
164 template<typename AFIELD>
165 void AFopr_Overlap<AFIELD>::tidyup()
166 {
167  delete m_sign;
168  if (m_kernel_created == true) delete m_fopr_w;
169 }
170 
171 
172 //====================================================================
173 template<typename AFIELD>
174 void AFopr_Overlap<AFIELD>::set_parameters(const Parameters& params)
175 {
176 #pragma omp barrier
177  int ith = ThreadManager::get_thread_id();
178 
179  std::string vlevel;
180  if (!params.fetch_string("verbose_level", vlevel)) {
181  if (ith == 0) m_vl = vout.set_verbose_level(vlevel);
182  }
183 #pragma omp barrier
184 
185  //- fetch and check input parameters
186  double mq, M0, x_min, x_max;
187  int Np;
188  std::vector<int> bc;
189 
190  int Niter;
191  double Stop_cond;
192 
193  int err = 0;
194  err += params.fetch_double("quark_mass", mq);
195  err += params.fetch_double("domain_wall_height", M0);
196  err += params.fetch_int("number_of_poles", Np);
197  err += params.fetch_double("lower_bound", x_min);
198  err += params.fetch_double("upper_bound", x_max);
199  err += params.fetch_int("maximum_number_of_iteration", Niter);
200  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
201  err += params.fetch_int_vector("boundary_condition", bc);
202 
203  if (err) {
204  vout.crucial(m_vl, "Error at %s: input parameter not found.\n",
205  class_name.c_str());
206  exit(EXIT_FAILURE);
207  }
208 
209  std::string repr, kernel;
210  if (!params.fetch_string("gamma_matrix_type", repr)) m_repr = repr;
211  if (!params.fetch_string("kernel_type", kernel)) m_kernel_type = kernel;
212 
213  set_parameters(real_t(mq), real_t(M0),
214  Np, real_t(x_min), real_t(x_max),
215  Niter, real_t(Stop_cond), bc);
216 
217  if (!m_is_initial_step) {
218  Parameters params_kernel = params;
219  double kappa = 1.0 / (8.0 - 2.0 * M0);
220  params_kernel.set_double("hopping_parameter", kappa);
221  m_fopr_w->set_parameters(params_kernel);
222 
223  m_sign->set_parameters(params);
224  }
225 }
226 
227 
228 //====================================================================
229 template<typename AFIELD>
230 void AFopr_Overlap<AFIELD>::set_parameters(const real_t mq,
231  const real_t M0,
232  const int Np,
233  const real_t x_min,
234  const real_t x_max,
235  const int Niter,
236  const real_t Stop_cond,
237  const std::vector<int> bc)
238 {
239 #pragma omp barrier
240 
241  int ith = ThreadManager::get_thread_id();
242 
243  //- range check
244  int err = 0;
245  err += ParameterCheck::non_zero(mq);
246  err += ParameterCheck::non_zero(M0);
247  err += ParameterCheck::non_zero(Np);
248  // NB. x_min,x_max == 0 is allowed.
249  err += ParameterCheck::non_zero(Niter);
250  err += ParameterCheck::square_non_zero(Stop_cond);
251 
252  if (err) {
253  vout.crucial(m_vl, "Error at %s: parameter range check failed.\n",
254  class_name.c_str());
255  exit(EXIT_FAILURE);
256  }
257 
258  const int Ndim = CommonParameters::Ndim();
259  assert(bc.size() == Ndim);
260 
261  if (ith == 0) {
262  m_mq = mq;
263  m_M0 = M0;
264  m_Np = Np;
265  m_x_min = x_min;
266  m_x_max = x_max;
267  m_Niter = Niter;
268  m_Stop_cond = Stop_cond;
269  m_boundary = bc;
270  }
271 #pragma omp barrier
272 
273  //- propagate parameters
274  // m_sign->set_parameters(m_Np, m_x_min, m_x_max, m_Niter, m_Stop_cond);
275 
276  //- print input parameters
277  vout.general(m_vl, "%s parameters:\n", class_name.c_str());
278  vout.general(m_vl, " mq = %12.8f\n", m_mq);
279  vout.general(m_vl, " M0 = %12.8f\n", m_M0);
280  vout.general(m_vl, " Np = %4d\n", m_Np);
281  vout.general(m_vl, " x_min = %12.8f\n", m_x_min);
282  vout.general(m_vl, " x_max = %12.8f\n", m_x_max);
283  vout.general(m_vl, " Niter = %6d\n", m_Niter);
284  vout.general(m_vl, " Stop_cond = %8.2e\n", m_Stop_cond);
285  for (int mu = 0; mu < Ndim; ++mu) {
286  vout.general(m_vl, " boundary[%d] = %2d\n", mu, m_boundary[mu]);
287  }
288 }
289 
290 
291 //====================================================================
292 template<typename AFIELD>
293 void AFopr_Overlap<AFIELD>::get_parameters(Parameters& params) const
294 {
295  params.set_double("quark_mass", double(m_mq));
296  params.set_double("domain_wall_height", double(m_M0));
297  params.set_int("number_of_poles", m_Np);
298  params.set_double("lower_bound", double(m_x_min));
299  params.set_double("upper_bound", double(m_x_max));
300  params.set_int("maximum_number_of_iteration", m_Niter);
301  params.set_double("convergence_criterion_squared", double(m_Stop_cond));
302  params.set_int_vector("boundary_condition", m_boundary);
303  params.set_string("kernel_type", m_kernel_type);
304  params.set_string("gamma_matrix_type", m_repr);
305 
306  params.set_string("verbose_level", vout.get_verbose_level(m_vl));
307 }
308 
309 
310 //====================================================================
311 template<typename AFIELD>
312 void AFopr_Overlap<AFIELD>::set_config(Field *U)
313 {
314  m_sign->set_config(U);
315 }
316 
317 
318 //====================================================================
319 template<typename AFIELD>
320 void AFopr_Overlap<AFIELD>::set_lowmodes(const int Nsbt,
321  std::vector<real_t> *ev,
322  std::vector<AFIELD> *vk)
323 {
324 #pragma omp barrier
325 
326  int ith = ThreadManager::get_thread_id();
327 
328  if (ith == 0) {
329  m_Nsbt = Nsbt;
330  m_ev = ev;
331  m_vk = vk;
332  }
333 
334 #pragma omp barrier
335 
336  m_sign->set_lowmodes(m_Nsbt, m_ev, m_vk);
337 }
338 
339 
340 //====================================================================
341 template<typename AFIELD>
342 void AFopr_Overlap<AFIELD>::set_mode(const std::string mode)
343 {
344 #pragma omp barrier
345 
346  int ith = ThreadManager::get_thread_id();
347  if (ith == 0) m_mode = mode;
348 
349  m_sign->set_mode(mode);
350 
351 #pragma omp barrier
352 }
353 
354 
355 //====================================================================
356 template<typename AFIELD>
357 void AFopr_Overlap<AFIELD>::mult(AFIELD& v, const AFIELD& w)
358 {
359  if (m_mode == "D") {
360  D(v, w);
361  } else if (m_mode == "Ddag") {
362  Ddag(v, w);
363  } else if (m_mode == "DdagD") {
364  DdagD(v, w);
365  } else if (m_mode == "H") {
366  H(v, w);
367  } else {
368  vout.crucial(m_vl, "Error at %s: mode undefined.\n",
369  class_name.c_str());
370  exit(EXIT_FAILURE);
371  }
372 }
373 
374 
375 //====================================================================
376 template<typename AFIELD>
377 void AFopr_Overlap<AFIELD>::mult_dag(AFIELD& v, const AFIELD& w)
378 {
379  if (m_mode == "D") {
380  Ddag(v, w);
381  } else if (m_mode == "Ddag") {
382  D(v, w);
383  } else if (m_mode == "DdagD") {
384  DdagD(v, w);
385  } else if (m_mode == "H") {
386  H(v, w);
387  } else {
388  vout.crucial(m_vl, "Error at %s: mode undefined.\n",
389  class_name.c_str());
390  exit(EXIT_FAILURE);
391  }
392 }
393 
394 
395 //====================================================================
396 template<typename AFIELD>
397 void AFopr_Overlap<AFIELD>::mult_gm5(AFIELD& v, const AFIELD& w)
398 {
399  m_fopr_w->mult_gm5(v, w);
400 }
401 
402 
403 //====================================================================
404 template<typename AFIELD>
405 void AFopr_Overlap<AFIELD>::DdagD(AFIELD& v, const AFIELD& b)
406 {
407  D(m_w1, b);
408  m_fopr_w->mult_gm5(m_w2, m_w1);
409 
410  D(m_w1, m_w2);
411  m_fopr_w->mult_gm5(v, m_w1);
412 }
413 
414 
415 //====================================================================
416 template<typename AFIELD>
417 void AFopr_Overlap<AFIELD>::H(AFIELD& v, const AFIELD& b)
418 {
419  D(m_w1, b);
420  m_fopr_w->mult_gm5(v, m_w1);
421 }
422 
423 
424 //====================================================================
425 template<typename AFIELD>
426 void AFopr_Overlap<AFIELD>::Ddag(AFIELD& v, const AFIELD& b)
427 {
428  m_fopr_w->mult_gm5(m_w1, b);
429  D(m_w2, m_w1);
430  m_fopr_w->mult_gm5(v, m_w2);
431 }
432 
433 
434 //====================================================================
435 template<typename AFIELD>
436 void AFopr_Overlap<AFIELD>::D(AFIELD& v, const AFIELD& b)
437 {
438  assert(b.check_size(m_Nin, m_Nvol, m_Nex));
439  assert(v.check_size(m_Nin, m_Nvol, m_Nex));
440 
441 #pragma omp barrier
442 
443  real_t f1 = m_M0 + 0.5 * m_mq;
444  real_t f2 = m_M0 - 0.5 * m_mq;
445 
446  m_sign->mult(m_w0, b);
447 
448  m_fopr_w->mult_gm5(v, m_w0);
449 
450  scal(v, f2);
451  axpy(v, f1, b);
452 
453 #pragma omp barrier
454 }
455 
456 
457 //====================================================================
458 template<typename AFIELD>
459 double AFopr_Overlap<AFIELD>::flop_count()
460 {
461  flop_count(m_mode);
462 }
463 
464 
465 //====================================================================
466 template<typename AFIELD>
467 double AFopr_Overlap<AFIELD>::flop_count(const std::string mode)
468 {
469  int NPE = CommonParameters::NPE();
470 
471  double gflop_sign = m_sign->flop_count();
472 
473  double flop_blas = double(m_Nin * m_Nex * 3);
474  // scal(1 flop per component) + axpy(2 flops per component)
475 
476  double gflop = flop_blas * double(m_Nvol) * double(NPE) * 1.0e-9;
477 
478  gflop += gflop_sign;
479 
480  if (mode == "DdagD") gflop *= 2.0;
481 
482  return gflop;
483 }
484 
485 
486 //============================================================END=====
AFopr_Overlap::real_t
AFIELD::real_t real_t
Definition: afopr_Overlap.h:49
AFopr_Overlap::mult_dag
void mult_dag(AFIELD &, const AFIELD &)
hermitian conjugate of mult.
Definition: afopr_Overlap-tmpl.h:377
AFopr_Overlap::get_parameters
void get_parameters(Parameters &params) const
gets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Overlap-tmpl.h:293
Parameters::set_string
void set_string(const string &key, const string &value)
Definition: parameters.cpp:39
AFopr
Definition: afopr.h:48
AFopr_Sign
Sign function of a given fermion operator.
Definition: afopr_Sign.h:56
AFopr_Overlap::mult
void mult(AFIELD &, const AFIELD &)
multiplies fermion operator to a given field.
Definition: afopr_Overlap-tmpl.h:357
CommonParameters::Ndim
static int Ndim()
Definition: commonParameters.h:117
Parameters
Class for parameters.
Definition: parameters.h:46
Parameters::set_double
void set_double(const string &key, const double value)
Definition: parameters.cpp:33
Bridge::BridgeIO::decrease_indent
void decrease_indent()
Definition: bridgeIO.h:86
Bridge::BridgeIO::increase_indent
void increase_indent()
Definition: bridgeIO.h:85
Field::check_size
bool check_size(const int nin, const int nvol, const int nex) const
checking size parameters. [23 May 2016 H.Matsufuru]
Definition: field.h:135
AFopr_Overlap::init
void init()
Definition: afopr_Overlap-tmpl.h:126
axpy
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:380
real_t
double real_t
Definition: bridgeQXS_Clover_coarse_double.cpp:16
AFopr_Overlap::set_lowmodes
void set_lowmodes(const int Nsbt, std::vector< real_t > *, std::vector< AFIELD > *)
Definition: afopr_Overlap-tmpl.h:320
AFopr_Overlap::tidyup
void tidyup()
Definition: afopr_Overlap-tmpl.h:165
AFopr_Overlap::D
void D(AFIELD &v, const AFIELD &w)
Definition: afopr_Overlap-tmpl.h:436
AFopr_Overlap::set_parameters
void set_parameters(const Parameters &params)
sets parameters by a Parameter object: to be implemented in a subclass.
Definition: afopr_Overlap-tmpl.h:174
AFopr_Overlap::mult_gm5
void mult_gm5(AFIELD &, const AFIELD &)
multiplies gamma_5 matrix.
Definition: afopr_Overlap-tmpl.h:397
Parameters::fetch_int_vector
int fetch_int_vector(const string &key, vector< int > &value) const
Definition: parameters.cpp:429
ParameterCheck::square_non_zero
int square_non_zero(const double v)
Definition: parameterCheck.cpp:43
threadManager.h
AFopr_Overlap::H
void H(AFIELD &v, const AFIELD &w)
Definition: afopr_Overlap-tmpl.h:417
AFopr_Overlap::set_mode
void set_mode(const std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: afopr_Overlap-tmpl.h:342
AFopr_Overlap::flop_count
double flop_count()
returns the number of floating point operations.
Definition: afopr_Overlap-tmpl.h:459
AFopr_Overlap::DdagD
void DdagD(AFIELD &v, const AFIELD &w)
Definition: afopr_Overlap-tmpl.h:405
CommonParameters::NPE
static int NPE()
Definition: commonParameters.h:101
Parameters::set_int_vector
void set_int_vector(const string &key, const vector< int > &value)
Definition: parameters.cpp:45
CommonParameters::Vlevel
static Bridge::VerboseLevel Vlevel()
Definition: commonParameters.h:122
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:133
ParameterCheck::non_zero
int non_zero(const double v)
Definition: parameterCheck.cpp:32
AFopr_Overlap::Ddag
void Ddag(AFIELD &v, const AFIELD &w)
Definition: afopr_Overlap-tmpl.h:426
Parameters::set_int
void set_int(const string &key, const int value)
Definition: parameters.cpp:36
scal
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:261
AFopr_Overlap
Overlap fermion operator.
Definition: afopr_Overlap.h:45
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:180
afopr_Overlap.h
Field
Container of Field-type object.
Definition: field.h:46
AFopr::field_nin
virtual int field_nin()=0
returns the on-site degree of freedom of the fermion field.
ThreadManager::get_thread_id
static int get_thread_id()
returns thread id.
Definition: threadManager.cpp:253
AFopr_Overlap::set_config
void set_config(Field *U)
sets the gauge configuration.
Definition: afopr_Overlap-tmpl.h:312
Parameters::fetch_int
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:200
ThreadManager::assert_single_thread
static void assert_single_thread(const std::string &class_name)
assert currently running on single thread.
Definition: threadManager.cpp:372
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:512
Bridge::BridgeIO::get_verbose_level
static std::string get_verbose_level(const VerboseLevel vl)
Definition: bridgeIO.cpp:154