root_cern.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file root_cern.h
24  \brief File defining \ref o2scl::root_cern
25 */
26 #ifndef O2SCL_ROOT_CERN_H
27 #define O2SCL_ROOT_CERN_H
28 
29 #include <string>
30 
31 #include <o2scl/string_conv.h>
32 #include <o2scl/root.h>
33 
34 namespace o2scl {
35 
36  /** \brief One-dimensional version of cern_mroot
37 
38  This one-dimensional root-finding routine, based on \ref
39  o2scl::mroot_cern, is probably slower than the more typical 1-d
40  routines, but also tends to converge for a larger class of
41  functions than \ref o2scl::root_bkt_cern, \ref
42  o2scl::root_brent_gsl, or \ref o2scl::root_stef. It has been
43  modified from \ref o2scl::mroot_cern and slightly optimized, but
44  has the same basic behavior.
45 
46  If \f$ x_i \f$ denotes the current iteration, and \f$
47  x^{\prime}_i \f$ denotes the previous iteration, then the
48  calculation is terminated if either (or both) of the following
49  tests is successful
50  \f[
51  1:\quad \mathrm{max} | f_i(x) | \leq \mathrm{tol\_rel}
52  \f]
53  \f[
54  2:\quad \mathrm{max} |x_i-x^{\prime}_i| \leq
55  \mathrm{tol\_abs} \times \mathrm{max} | x_i |
56  \f]
57 
58  \note This class has difficulty finding the root when
59  the desired root is near 0 (AWS 1/22/19)
60 
61  \note This code has not been checked to ensure that it cannot
62  fail to solve the equations without calling the error handler
63  and returning a non-zero value. Until then, the solution may
64  need to be checked explicitly by the caller.
65 
66  See the \ref onedsolve_subsect section of the User's guide for
67  general information about \o2 solvers.
68 
69  \future Double-check this class to make sure it cannot fail
70  while returning 0 for success.
71  */
72 #ifdef DOXYGEN_NO_O2NS
73  template<class func_t=funct> class root_cern : public root
74 #else
75  template<class func_t=funct> class root_cern :
76  public root<func_t>
77 #endif
78  {
79 
80  public:
81 
82  root_cern() {
83  info=0;
84  eps=0.1490116119384766e-07;
85  scale=10.0;
86  maxf=0;
87  }
88 
89  virtual ~root_cern() {}
90 
91  /** \brief Get the value of \c INFO from the last call to solve()
92  (default 0)
93 
94  The value of info is assigned according to the following list.
95  The values 1-8 are the standard behavior from CERNLIB.
96  0 - The function solve() has not been called.
97  1 - Test 1 was successful. \n
98  2 - Test 2 was successful. \n
99  3 - Both tests were successful. \n
100  4 - Number of iterations is greater than root_cern::maxf. \n
101  5 - Approximate (finite difference) Jacobian matrix is
102  singular. \n
103  6 - Iterations are not making good progress. \n
104  7 - Iterations are diverging. \n
105  8 - Iterations are converging, but either root_cern::tol_abs
106  is too small or the Jacobian is nearly singular
107  or the variables are badly scaled. \n
108  9 - Either root::tol_rel or root::tol_abs is not greater than zero.
109 
110  */
111  int get_info() { return info; }
112 
113  /** \brief Maximum number of function evaluations
114 
115  If \f$ \mathrm{maxf}\leq 0 \f$, then 200 (which is the CERNLIB
116  default) is used. The default value of \c maxf is zero which
117  then implies the default from CERNLIB.
118  */
119  int maxf;
120 
121  /// Return the type, \c "root_cern".
122  virtual const char *type() { return "root_cern"; }
123 
124  /** \brief The original scale parameter from CERNLIB (default 10.0)
125  */
126  double scale;
127 
128  /** \brief The smallest floating point number
129  (default \f$ \sim 1.49012 \times 10^{-8} \f$)
130 
131  The original prescription from CERNLIB for \c eps is
132  given below:
133  \verbatim
134  #if !defined(CERNLIB_DOUBLE)
135  PARAMETER (EPS = 0.84293 69702 17878 97282 52636 392E-07)
136  #endif
137  #if defined(CERNLIB_IBM)
138  PARAMETER (EPS = 0.14901 16119 38476 562D-07)
139  #endif
140  #if defined(CERNLIB_VAX)
141  PARAMETER (EPS = 0.37252 90298 46191 40625D-08)
142  #endif
143  #if (defined(CERNLIB_UNIX))&&(defined(CERNLIB_DOUBLE))
144  PARAMETER (EPS = 0.14901 16119 38476 600D-07)
145  #endif
146  \endverbatim
147 
148  \future This number should probably default to one of the
149  GSL tolerances.
150  */
151  double eps;
152 
153  /// Solve \c func using \c x as an initial guess, returning \c x.
154  virtual int solve(double &ux, func_t &func) {
155 
156  int it=0;
157  double fky;
158 
159  int lmaxf;
160  if (maxf<=0) lmaxf=200;
161  else lmaxf=maxf;
162 
163  info=0;
164 
165  if (this->tol_rel<=0.0 || this->tol_abs<=0.0) {
166  info=9;
167  std::string str="Invalid value of tol_rel ("+dtos(this->tol_rel)+
168  ") or tol_abs ("+dtos(this->tol_abs)+") in root_cern::solve().";
169  O2SCL_ERR(str.c_str(),exc_einval);
170  }
171 
172  int iflag=0, numf=0, nfcall=0, nier6=-1, nier7=-1, nier8=0;
173  double fnorm=0.0, difit=0.0, xnorm=0.0;
174  bool set=false;
175 
176  if (xnorm<fabs(ux)) {
177  xnorm=fabs(ux);
178  set=true;
179  }
180 
181  double delta=scale*xnorm;
182  if (set==false) delta=scale;
183 
184  double wmat, farr, w0arr, w1arr, w2arr;
185 
186  bool solve_done=false;
187  while (solve_done==false) {
188 
189  int nsing=1;
190  double fnorm1=fnorm;
191  double difit1=difit;
192  fnorm=0.0;
193 
194  // Compute step H for the divided difference which approximates
195  // the K-th row of the Jacobian matrix
196 
197  double h=eps*xnorm;
198  if (h==0.0) h=eps;
199 
200  wmat=h;
201  w1arr=ux;
202 
203  // Enter a subiteration
204 
205  iflag=0;
206 
207  farr=func(w1arr);
208 
209  fky=farr;
210  nfcall++;
211  numf=nfcall;
212  if (fnorm<fabs(fky)) fnorm=fabs(fky);
213 
214  // Compute the K-th row of the Jacobian matrix
215 
216  w2arr=w1arr+wmat;
217 
218  farr=func(w2arr);
219 
220  double fkz=farr;
221  nfcall++;
222  numf=nfcall;
223  w0arr=fkz-fky;
224 
225  farr=fky;
226 
227  // Compute the Householder transformation to reduce the K-th row
228  // of the Jacobian matrix to a multiple of the K-th unit vector
229 
230  double eta=0.0;
231  if (eta<fabs(w0arr)) eta=fabs(w0arr);
232 
233  if (eta!=0.0) {
234  nsing--;
235  double sknorm=0.0;
236 
237  w0arr/=eta;
238  sknorm+=w0arr*w0arr;
239 
240  sknorm=sqrt(sknorm);
241  if (w0arr<0.0) sknorm=-sknorm;
242  w0arr+=sknorm;
243 
244  // Apply the transformation
245 
246  w2arr=0.0;
247  w2arr+=w0arr*wmat;
248  double temp=w0arr/(sknorm*w0arr);
249  wmat-=temp*w2arr;
250 
251  // Compute the subiterate
252 
253  w0arr=sknorm*eta;
254  double temp2=fky/w0arr;
255  if (h*fabs(temp2)>delta)
256  temp2=(temp2>=0.0) ? fabs(delta/h) : -fabs(delta/h);
257  w1arr+=temp2*wmat;
258  }
259 
260  // Compute the norms of the iterate and correction vector
261 
262  xnorm=0.0;
263  difit=0.0;
264 
265  if (xnorm<fabs(w1arr)) xnorm=fabs(w1arr);
266  if (difit<fabs(ux-w1arr)) difit=fabs(ux-w1arr);
267  ux=w1arr;
268 
269  // Update the bound on the correction vector
270 
271  if(delta<scale*xnorm) delta=scale*xnorm;
272 
273  // Determine the progress of the iteration
274 
275  bool lcv=(fnorm<fnorm1 && difit<difit1 && nsing==0);
276  nier6++;
277  nier7++;
278  nier8++;
279  if (lcv) nier6=0;
280  if (fnorm<fnorm1 || difit<difit1) nier7=0;
281  if (difit>eps*xnorm) nier8=0;
282 
283  // Print iteration information
284 
285  if (this->verbose>0) {
286  this->print_iter(ux,farr,++it,fnorm,this->tol_rel,
287  "root_cern");
288  }
289 
290  // Tests for convergence
291 
292  if (fnorm<=this->tol_rel) info=1;
293  if (difit<=this->tol_abs*xnorm && lcv) info=2;
294  if (fnorm<=this->tol_rel && info==2) info=3;
295  if (info!=0) {
296  if (!std::isfinite(ux)) {
297  O2SCL_CONV2_RET("Solver converged to non-finite value ",
298  "in root_cern::solve().",exc_erange,
299  this->err_nonconv);
300  }
301  return 0;
302  }
303 
304  // Tests for termination
305 
306  if (numf>=lmaxf) {
307  info=4;
308  O2SCL_CONV_RET("Too many iterations in root_cern::solve().",
309  exc_emaxiter,this->err_nonconv);
310  }
311  if (nsing==1) {
312  info=5;
313  O2SCL_CONV2_RET("Jacobian matrix singular in ",
314  "root_cern::solve().",
315  exc_esing,this->err_nonconv);
316  }
317  if (nier6==5) {
318  info=6;
319  O2SCL_CONV_RET("No progress in root_cern::solve().",
320  exc_enoprog,this->err_nonconv);
321  }
322  if (nier7==3) {
323  info=7;
324  O2SCL_CONV_RET("Iterations diverging in root_cern::solve().",
325  exc_erunaway,this->err_nonconv);
326  }
327  if (nier8==4) {
328  info=8;
329  std::string s2="Variable tol_abs too small, J singular, ";
330  s2+="or bad scaling in root_cern::solve().";
331  O2SCL_CONV(s2.c_str(),exc_efailed,this->err_nonconv);
332  }
333 
334  // Exit if necessary
335 
336  if (info!=0) return exc_efailed;
337 
338  }
339 
340  if (!std::isfinite(ux)) {
341  O2SCL_CONV2_RET("Solver converged to non-finite value ",
342  "in root_cern::solve() (2).",exc_erange,
343  this->err_nonconv);
344  }
345  return 0;
346  }
347 
348 #ifndef DOXYGEN_INTERNAL
349 
350  protected:
351 
352  /// Internal storage for the value of \c info
353  int info;
354 
355 #endif
356 
357  };
358 
359 }
360 
361 #endif
o2scl::root< funct >::tol_rel
double tol_rel
The maximum value of the functions for success (default )
Definition: root.h:66
o2scl::exc_esing
@ exc_esing
apparent singularity detected
Definition: err_hnd.h:93
o2scl::exc_erange
@ exc_erange
output range error, e.g. exp(1e100)
Definition: err_hnd.h:55
O2SCL_CONV2_RET
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:303
o2scl::dtos
std::string dtos(const fp_t &x, int prec=6, bool auto_prec=false)
Convert a double to a string.
Definition: string_conv.h:77
o2scl::exc_efailed
@ exc_efailed
generic failure
Definition: err_hnd.h:61
o2scl::root< funct >::tol_abs
double tol_abs
The minimum allowable stepsize (default )
Definition: root.h:71
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::root_cern::eps
double eps
The smallest floating point number (default )
Definition: root_cern.h:151
o2scl::root_cern::get_info
int get_info()
Get the value of INFO from the last call to solve() (default 0)
Definition: root_cern.h:111
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::root< funct >::print_iter
virtual int print_iter(double x, double y, int iter, double value=0.0, double limit=0.0, std::string comment="")
Print out iteration information.
Definition: root.h:100
o2scl::root_cern::maxf
int maxf
Maximum number of function evaluations.
Definition: root_cern.h:119
o2scl::root< funct >::verbose
int verbose
Output control (default 0)
Definition: root.h:74
o2scl::root
One-dimensional solver [abstract base].
Definition: root.h:48
o2scl::root_cern::solve
virtual int solve(double &ux, func_t &func)
Solve func using x as an initial guess, returning x.
Definition: root_cern.h:154
o2scl::root_cern::scale
double scale
The original scale parameter from CERNLIB (default 10.0)
Definition: root_cern.h:126
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
O2SCL_CONV_RET
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:297
o2scl::exc_enoprog
@ exc_enoprog
iteration is not making progress toward solution
Definition: err_hnd.h:105
o2scl::root_cern::type
virtual const char * type()
Return the type, "root_cern".
Definition: root_cern.h:122
o2scl::root_cern
One-dimensional version of cern_mroot.
Definition: root_cern.h:75
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::exc_erunaway
@ exc_erunaway
iterative process is out of control
Definition: err_hnd.h:71
O2SCL_CONV
#define O2SCL_CONV(d, n, b)
Set a "convergence" error.
Definition: err_hnd.h:277
o2scl::root_cern::info
int info
Internal storage for the value of info.
Definition: root_cern.h:353
o2scl::root< funct >::err_nonconv
bool err_nonconv
If true, call the error handler if the solver does not converge (default true)
Definition: root.h:82

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).