Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fopr_Rational_SF.cpp
Go to the documentation of this file.
1 
14 #include "fopr_Rational_SF.h"
15 
16 #ifdef USE_FACTORY_AUTOREGISTER
17 namespace {
18  bool init = Fopr_Rational_SF::register_factory();
19 }
20 #endif
21 
22 const std::string Fopr_Rational_SF::class_name = "Fopr_Rational_SF";
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  //- fetch and check input parameters
32  int Np, n_exp, d_exp;
33  double x_min, x_max;
34  int Niter;
35  double Stop_cond;
36 
37  int err = 0;
38  err += params.fetch_int("number_of_poles", Np);
39  err += params.fetch_int("exponent_numerator", n_exp);
40  err += params.fetch_int("exponent_denominator", d_exp);
41  err += params.fetch_double("lower_bound", x_min);
42  err += params.fetch_double("upper_bound", x_max);
43  err += params.fetch_int("maximum_number_of_iteration", Niter);
44  err += params.fetch_double("convergence_criterion_squared", Stop_cond);
45 
46  if (err) {
47  vout.crucial(m_vl, "Error at %s: input parameter not found.\n", class_name.c_str());
48  exit(EXIT_FAILURE);
49  }
50 
51 
52  set_parameters(Np, n_exp, d_exp, x_min, x_max, Niter, Stop_cond);
53 }
54 
55 
56 //====================================================================
57 void Fopr_Rational_SF::set_parameters(const int Np, const int n_exp, const int d_exp,
58  const double x_min, const double x_max,
59  const int Niter, const double Stop_cond)
60 {
61  //- print input parameters
62  vout.general(m_vl, "%s:\n", class_name.c_str());
63  vout.general(m_vl, " Np = %d\n", Np);
64  vout.general(m_vl, " n_exp = %d\n", n_exp);
65  vout.general(m_vl, " d_exp = %d\n", d_exp);
66  vout.general(m_vl, " x_min = %12.8f\n", x_min);
67  vout.general(m_vl, " x_max = %12.8f\n", x_max);
68  vout.general(m_vl, " Niter = %d\n", Niter);
69  vout.general(m_vl, " Stop_cond = %8.2e\n", Stop_cond);
70 
71  //- range check
72  int err = 0;
73  err += ParameterCheck::non_zero(Np);
74  err += ParameterCheck::non_zero(n_exp);
75  err += ParameterCheck::non_zero(d_exp);
76  // NB. x_min,x_max == 0 is allowed.
77  err += ParameterCheck::non_zero(Niter);
78  err += ParameterCheck::square_non_zero(Stop_cond);
79 
80  if (err) {
81  vout.crucial(m_vl, "Error at %s: parameter range check failed.\n", class_name.c_str());
82  exit(EXIT_FAILURE);
83  }
84 
85  //- store values
86  m_Np = Np;
87  m_n_exp = n_exp;
88  m_d_exp = d_exp;
89  m_x_min = x_min;
90  m_x_max = x_max;
91  m_Niter = Niter;
92  m_Stop_cond = Stop_cond;
93 
94  //- post-process
96 }
97 
98 
99 //====================================================================
101 {
102  const int Nin = m_fopr->field_nin();
103  const int Nvol = m_fopr->field_nvol();
104  const int Nex = m_fopr->field_nex();
105 
106  const int Nshift = m_Np;
107 
108  const double x_min2 = m_x_min * m_x_min;
109  const double x_max2 = m_x_max * m_x_max;
110 
111  m_cl.resize(m_Np);
112  m_bl.resize(m_Np);
113 
114  m_xq.resize(m_Np);
115  for (int i = 0; i < Nshift; ++i) {
116  m_xq[i].reset(Nin, Nvol, Nex);
117  }
118 
120 
121  Math_Rational rational;
122  rational.set_parameters(m_Np, m_n_exp, m_d_exp, x_min2, x_max2);
123  rational.get_parameters(m_a0, m_bl, m_cl);
124 
125  vout.general(m_vl, " a0 = %18.14f\n", m_a0);
126  for (int i = 0; i < m_Np; i++) {
127  vout.general(m_vl, " bl[%d] = %18.14f cl[%d] = %18.14f\n",
128  i, m_bl[i], i, m_cl[i]);
129  }
130 }
131 
132 
133 //====================================================================
134 
143 {
144  assert(v.nin() == b.nin());
145  assert(v.nvol() == b.nvol());
146  assert(v.nex() == b.nex());
147 
148  copy(v, b);
149 
150  vout.general(m_vl, " Shift solver in rational function\n");
151  vout.general(m_vl, " Number of shift values = %d\n", m_Np);
152 
153  Field_F_SF set_zero;
154  set_zero.set_boundary_zero(v);
155 
156  int Nconv;
157  double diff;
158 
159  m_fopr->set_mode("DdagD");
160  m_solver->solve(m_xq, m_cl, v, Nconv, diff);
161 
162  vout.general(m_vl, " diff(max) = %22.15e \n", diff);
163 
164  copy(v, b); // v = b;
165  scal(v, m_a0); // v *= m_a0;
166  for (int i = 0; i < m_Np; i++) {
167  axpy(v, m_bl[i], m_xq[i]); // v += m_bl[i] * m_xq[i];
168  }
169 
170  set_zero.set_boundary_zero(v);
171 }
172 
173 
174 //====================================================================
175 double Fopr_Rational_SF::func(const double x)
176 {
177  double y = m_a0;
178 
179  for (int k = 0; k < m_Np; ++k) {
180  y += m_bl[k] / (x + m_cl[k]);
181  }
182 
183  return y;
184 }
185 
186 
187 //====================================================================
188 //============================================================END=====
void scal(Field &x, const double a)
scal(x, a): x = a * x
Definition: field.cpp:433
void get_parameters(double &norm, std::vector< double > &res, std::vector< double > &pole)
BridgeIO vout
Definition: bridgeIO.cpp:503
double func(const double x)
void general(const char *format,...)
Definition: bridgeIO.cpp:197
void set_parameters(const Parameters &params)
virtual void set_mode(const std::string mode)
setting the mode of multiplication if necessary. Default implementation here is just to avoid irrelev...
Definition: fopr.h:94
Container of Field-type object.
Definition: field.h:45
int fetch_double(const string &key, double &value) const
Definition: parameters.cpp:327
int nvol() const
Definition: field.h:127
Multishift Conjugate Gradient solver.
static const std::string class_name
Class for parameters.
Definition: parameters.h:46
void copy(Field &y, const Field &x)
copy(y, x): y = x
Definition: field.cpp:532
std::vector< Field > m_xq
virtual int field_nin()=0
returns the on-site d.o.f. for which the fermion operator is defined.
Determionation of coefficients of rational approximation.
Definition: math_Rational.h:40
int square_non_zero(const double v)
int nin() const
Definition: field.h:126
void set_parameters(const Parameters &params)
void mult(Field &v, const Field &f)
Bridge::VerboseLevel m_vl
Definition: fopr.h:127
int fetch_int(const string &key, int &value) const
Definition: parameters.cpp:346
virtual int field_nex()=0
returns the external d.o.f. for which the fermion operator is defined.
int nex() const
Definition: field.h:128
void axpy(Field &y, const double a, const Field &x)
axpy(y, a, x): y := a * x + y
Definition: field.cpp:319
void crucial(const char *format,...)
Definition: bridgeIO.cpp:178
Shiftsolver_CG * m_solver
void set_boundary_zero(Field &f)
Definition: field_F_SF.h:53
int non_zero(const double v)
string get_string(const string &key) const
Definition: parameters.cpp:221
virtual int field_nvol()=0
returns the volume for which the fermion operator is defined.
std::vector< double > m_cl
A class generated to add a function for the SF.
Definition: field_F_SF.h:33
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131
std::vector< double > m_bl
void solve(std::vector< Field > &solution, const std::vector< double > &shift, const Field &source, int &Nconv, double &diff)