root_bkt_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_bkt_cern.h
24  \brief File defining \ref o2scl::root_bkt_cern
25 */
26 #ifndef O2SCL_ROOT_BKT_CERN_H
27 #define O2SCL_ROOT_BKT_CERN_H
28 
29 #include <o2scl/root.h>
30 
31 #ifndef DOXYGEN_NO_O2NS
32 namespace o2scl {
33 #endif
34 
35  /** \brief One-dimensional root-finding routine (CERNLIB)
36 
37  This class attempts to find \f$ x_1 \f$ and \f$ x_2 \f$ in \f$
38  [a,b] \f$ such that:
39  - \f$ f(x_1) f(x_2) \leq 0 \f$,
40  - \f$ |f(x_1)| \leq|f(x_2)| \f$, and
41  - \f$ | x_1-x_2| \leq 2~\mathrm{tol\_abs}~(1+|x_0|) \f$.
42 
43  The function solve_bkt() requires inputs \c x1 and \c x2 such
44  that the first condition, \f$ f(x_1) f(x_2) \leq 0 \f$, already
45  holds.
46 
47  The variable root::tol_abs defaults to \f$ 10^{-8} \f$ and
48  root::ntrial defaults to 200.
49  \comment
50  I'm not sure why I chose these values for tol_abs and ntrial above,
51  as it doesn't seem to me that these values are suggested in
52  CERNLIB. However, they're reasonable so I'll leave them in for
53  now.
54  \endcomment
55 
56  The function solve_bkt() will call the error handler if the root
57  is not initially bracketed. If \ref root::err_nonconv is true
58  (as it is by default), then the error handler will also be
59  called if the number function evaluations is greater than
60  root::ntrial.
61 
62  After a call to \ref solve_bkt(), \ref root::last_ntrial
63  contains the total number of iterations which were used
64 
65 
66  See the \ref onedsolve_subsect section of the User's guide for
67  general information about \o2 solvers.
68 
69  Based on the CERNLIB routines RZEROX and DZEROX, which was
70  based on \ref Bus75 and is documented at
71  http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/c200/top.html
72  */
73  template<class func_t=funct11> class root_bkt_cern :
74  public root_bkt<func_t> {
75 
76 #ifndef DOXYGEN_INTERNAL
77 
78  protected:
79 
80  /** \brief Internal storage for the mode
81 
82  This internal variable is actually defined to be smaller by
83  1 than the "mode" as it is defined in the CERNLIB
84  documentation in order to avoid needless subtraction in
85  solve_bkt().
86  */
87  int mode;
88 
89  /// FORTRAN-like function for sign
90  inline double sign(double a, double b) {
91  if (b>=0.0) return fabs(a);
92  return -fabs(a);
93  }
94 
95 #endif
96 
97  public:
98 
99  root_bkt_cern() {
100  mode=0;
101  this->tol_abs=1.0e-8;
102  this->ntrial=200;
103  }
104 
105  /** \brief Set mode of solution (1 or 2)
106 
107  - \c 1 should be used for simple functions where the cost is
108  inexpensive in comparison to one iteration of solve_bkt(),
109  or functions which have a pole near the root (this is the
110  default).
111  - \c 2 should be used for more time-consuming functions.
112 
113  If an integer other than \c 1 or \c 2 is specified,
114  the error handler is called.
115  */
116  int set_mode(int m) {
117  if (m!=1 && m!=2) {
118  O2SCL_ERR("Invalid mode in root_bkt_cern::set_mode().",
120  }
121  mode=m-1;
122  return 0;
123  }
124 
125  /// Return the type, \c "root_bkt_cern".
126  virtual const char *type() { return "root_bkt_cern"; }
127 
128  /** \brief Solve \c func in region \f$ x_1<x<x_2 \f$ returning
129  \f$ x_1 \f$.
130  */
131  virtual int solve_bkt(double &x1, double x2, func_t &func) {
132 
133  double im1[2]={2.0,3.0}, im2[2]={-1.0,3.0}, c=0.0, fa, fb;
134  double atl, a, b, mft;
135  double fc=0.0, d=0.0, fd=0.0, tol, h, hb, w, p, q, fdb, fda, f=0.0;
136  bool lmt[2];
137  int ie=0, loop, mf;
138  char ch;
139 
140  fb=func(x1);
141  fa=func(x2);
142 
143  if (fa*fb>0.0) {
144  O2SCL_ERR2("Endpoints don't bracket function in ",
145  "root_bkt_cern::solve_bkt().",exc_einval);
146  }
147 
148  atl=fabs(this->tol_abs);
149  b=x1;
150  a=x2;
151  lmt[1]=true;
152  mf=2;
153  loop=1;
154  do {
155  if (loop==1) {
156  c=a;
157  fc=fa;
158  ie=0;
159  } else if (loop==2) {
160  ie=0;
161  }
162  if (fabs(fc)<fabs(fb)) {
163  if (c!=a) {
164  d=a;
165  fd=fa;
166  }
167  a=b;
168  b=c;
169  c=a;
170  fa=fb;
171  fb=fc;
172  fc=fa;
173  }
174  tol=atl*(1.0+fabs(c));
175  h=0.5*(c+b);
176  hb=h-b;
177 
178  if (this->verbose>0) {
179  this->print_iter(c,fc,mf-2,fabs(hb),tol,"root_bkt_cern");
180  }
181 
182  if (fabs(hb)>tol) {
183  if (ie>im1[mode]) {
184  w=hb;
185  } else {
186  tol*=sign(1.0,hb);
187  p=(b-a)*fb;
188  lmt[0]=(ie<=1);
189  if (lmt[mode]) {
190  q=fa-fb;
191  lmt[1]=false;
192  } else {
193  fdb=(fd-fb)/(d-b);
194  fda=(fd-fa)/(d-a);
195  p*=fda;
196  q=fdb*fa-fda*fb;
197  }
198  if (p<0.0) {
199  p=-p;
200  q=-q;
201  }
202  if (ie==im2[mode]) p+=p;
203  if (p==0.0 || p<=q*tol) {
204  w=tol;
205  } else if (p<hb*q) {
206  w=p/q;
207  } else {
208  w=hb;
209  }
210  }
211  d=a;
212  a=b;
213  fd=fa;
214  fa=fb;
215  b+=w;
216  mf++;
217  if (mf>this->ntrial) {
218  this->last_ntrial=this->ntrial;
219  O2SCL_CONV2_RET("Too many function calls in ",
220  "root_bkt_cern::solve_bkt().",exc_emaxiter,
221  this->err_nonconv);
222  }
223 
224  fb=func(b);
225 
226  if (fb==0.0 || sign(1.0,fc)==sign(1.0,fb)) {
227  loop=1;
228  } else if (w==hb) {
229  loop=2;
230  } else {
231  ie++;
232  loop=3;
233  }
234  } else {
235  loop=0;
236  }
237  } while (loop>0);
238  x1=c;
239  this->last_ntrial=mf;
240  return o2scl::success;
241  }
242 
243  };
244 
245 #ifndef DOXYGEN_NO_O2NS
246 }
247 #endif
248 
249 #endif
One-dimensional root-finding routine (CERNLIB)
Definition: root_bkt_cern.h:73
int ntrial
Maximum number of iterations (default 100)
Definition: root.h:76
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
int mode
Internal storage for the mode.
Definition: root_bkt_cern.h:87
double tol_abs
The minimum allowable stepsize (default )
Definition: root.h:70
invalid argument supplied by user
Definition: err_hnd.h:59
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
One-dimensional bracketing solver [abstract base].
Definition: root.h:145
bool err_nonconv
If true, call the error handler if the solver does not converge (default true)
Definition: root.h:81
virtual int solve_bkt(double &x1, double x2, func_t &func)
Solve func in region returning .
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
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
int set_mode(int m)
Set mode of solution (1 or 2)
double sign(double a, double b)
FORTRAN-like function for sign.
Definition: root_bkt_cern.h:90
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
int last_ntrial
The number of iterations used in the most recent solve.
Definition: root.h:84
static const double x2[5]
Definition: inte_qng_gsl.h:66
virtual const char * type()
Return the type, "root_bkt_cern".
static const double x1[5]
Definition: inte_qng_gsl.h:48
int verbose
Output control (default 0)
Definition: root.h:73
Success.
Definition: err_hnd.h:47

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