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

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