root_bkt_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_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=funct> 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
o2scl::root_bkt_cern
One-dimensional root-finding routine (CERNLIB)
Definition: root_bkt_cern.h:73
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::root_bkt_cern::set_mode
int set_mode(int m)
Set mode of solution (1 or 2)
Definition: root_bkt_cern.h:116
o2scl::root< funct, funct, double >::last_ntrial
int last_ntrial
The number of iterations used in the most recent solve.
Definition: root.h:85
o2scl::root< funct, funct, double >::tol_abs
double tol_abs
The minimum allowable stepsize (default )
Definition: root.h:71
o2scl::root_bkt
One-dimensional bracketing solver [abstract base].
Definition: root.h:151
O2SCL_ERR2
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
o2scl::root_bkt_cern::solve_bkt
virtual int solve_bkt(double &x1, double x2, func_t &func)
Solve func in region returning .
Definition: root_bkt_cern.h:131
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::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
o2scl::root_bkt_cern::type
virtual const char * type()
Return the type, "root_bkt_cern".
Definition: root_bkt_cern.h:126
o2scl::root< funct, funct, double >::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< funct, funct, double >::verbose
int verbose
Output control (default 0)
Definition: root.h:74
o2scl_inte_qng_coeffs::x2
static const double x2[5]
Definition: inte_qng_gsl.h:66
o2scl::success
@ success
Success.
Definition: err_hnd.h:47
o2scl_inte_qng_coeffs::x1
static const double x1[5]
Definition: inte_qng_gsl.h:48
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::root_bkt_cern::sign
double sign(double a, double b)
FORTRAN-like function for sign.
Definition: root_bkt_cern.h:90
o2scl::root_bkt_cern::mode
int mode
Internal storage for the mode.
Definition: root_bkt_cern.h:87
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::root< funct, funct, double >::err_nonconv
bool err_nonconv
If true, call the error handler if the solver does not converge (default true)
Definition: root.h:82
o2scl::root< funct, funct, double >::ntrial
int ntrial
Maximum number of iterations (default 100)
Definition: root.h:77

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