root_brent_gsl.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 /* roots/brent.c
24  *
25  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky,
26  * Brian Gough
27  *
28  * This program is free software; you can redistribute it and/or modify
29  * it under the terms of the GNU General Public License as published by
30  * the Free Software Foundation; either version 3 of the License, or (at
31  * your option) any later version.
32  *
33  * This program is distributed in the hope that it will be useful, but
34  * WITHOUT ANY WARRANTY; without even the implied warranty of
35  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
36  * General Public License for more details.
37  *
38  * You should have received a copy of the GNU General Public License
39  * along with this program; if not, write to the Free Software
40  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41  * 02110-1301, USA.
42  */
43 
44 #ifndef O2SCL_ROOT_BRENT_GSL_H
45 #define O2SCL_ROOT_BRENT_GSL_H
46 
47 /** \file root_brent_gsl.h
48  \brief File defining \ref o2scl::root_brent_gsl
49 */
50 
51 #include <limits>
52 #include <gsl/gsl_math.h>
53 #include <gsl/gsl_roots.h>
54 #include <o2scl/funct.h>
55 #include <o2scl/root.h>
56 
57 #ifndef DOXYGEN_NO_O2NS
58 namespace o2scl {
59 #endif
60 
61  /** \brief One-dimensional root-finding (GSL)
62 
63  This class finds the root of a user-specified function. If \ref
64  test_form is 0, then solve_bkt() stops when the size of the
65  bracket is smaller than \ref root::tol_abs. If \ref test_form is
66  1, then the function stops when the residual is less than \ref
67  root::tol_rel. If test_form is 2, then both tests are applied.
68 
69  See the \ref onedsolve_subsect section of the User's guide for
70  general information about \o2 solvers. An example demonstrating
71  the usage of this class is given in
72  <tt>examples/ex_fptr.cpp</tt> and the \ref ex_fptr_sect .
73 
74  \future There is some duplication in the variables \c x_lower,
75  \c x_upper, \c a, and \c b, which could be removed. Some
76  better variable names would also be helpful.
77 
78  \future Create a meaningful enum list for \ref
79  o2scl::root_brent_gsl::test_form.
80 
81  \comment
82  Note that I used \c instead of \ref to refer to variables above
83  since the variables a and b are protected, and not available if
84  compiling the documentation without the internal portion.
85 
86  Also, at first I got confused that GSL didn't require
87  lower<upper, but it turns out that this is indeed a requirement
88  in GSL, but I didn't see it because it was in roots/fsolver.c
89  rather than in roots/brent.c . Thus, everything looks fine now.
90  \endcomment
91  */
92  template<class func_t=funct, class fp_t=double> class root_brent_gsl :
93  public root_bkt<func_t,func_t,fp_t> {
94 
95  protected:
96 
97  /** \brief Floating point-type agnostic version of
98  \c gsl_root_test_interval() .
99  */
100  int test_interval(fp_t xx_lower, fp_t xx_upper, fp_t epsabs,
101  fp_t epsrel) {
102  fp_t abs_lower, abs_upper;
103 
104  if (xx_lower<0.0) abs_lower=-xx_lower;
105  else abs_lower=xx_lower;
106  if (xx_upper<0.0) abs_upper=-xx_upper;
107  else abs_upper=xx_upper;
108 
109  fp_t min_abs, tolerance;
110  if (epsrel<0.0) {
111  O2SCL_ERR2("Relative tolerance is negative in ",
112  "root_brent_gsl::test_interval().",o2scl::exc_ebadtol);
113  }
114  if (epsabs<0.0) {
115  O2SCL_ERR2("Absolute tolerance is negative in ",
116  "root_brent_gsl::test_interval().",o2scl::exc_ebadtol);
117  }
118  if (xx_lower>xx_upper) {
119  O2SCL_ERR2("Lower bound larger than upper bound in ",
120  "root_brent_gsl::test_interval().",o2scl::exc_einval);
121  }
122 
123  if ((xx_lower>0.0 && xx_upper>0.0) ||
124  (xx_lower<0.0 && xx_upper<0.0)) {
125  if (abs_lower<abs_upper) min_abs=abs_lower;
126  else min_abs=abs_upper;
127  } else {
128  min_abs=0;
129  }
130 
131  tolerance=epsabs+epsrel*min_abs;
132 
133  // AWS: I could combine this if statement but this form ensures
134  // success in case floating point problems imply that
135  // xx_lower<xx_upper and xx_lower>=xx_upper are both false.
136 
137  if (xx_lower<xx_upper) {
138  if (xx_upper-xx_lower<tolerance) {
139  return o2scl::success;
140  }
141  } else {
142  if (xx_lower-xx_upper<tolerance) {
143  return o2scl::success;
144  }
145  }
146 
147  return o2scl::gsl_continue;
148  }
149 
150  public:
151 
152  root_brent_gsl() {
153  test_form=0;
154  }
155 
156  /// Return the type, \c "root_brent_gsl".
157  virtual const char *type() { return "root_brent_gsl"; }
158 
159  /** \brief Perform an iteration
160 
161  This function currently always returns \ref success.
162  */
163  int iterate(func_t &f) {
164 
165  fp_t tol, m, two=2;
166 
167  int ac_equal=0;
168 
169  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
170  ac_equal=1;
171  c=a;
172  fc=fa;
173  d=b-a;
174  e=b-a;
175  }
176 
177  if (o2scl::o2abs(fc) < o2scl::o2abs(fb)) {
178  ac_equal=1;
179  a=b;
180  b=c;
181  c=a;
182  fa=fb;
183  fb=fc;
184  fc=fa;
185  }
186 
187  tol=o2scl::o2abs(b)*std::numeric_limits<fp_t>::epsilon()/two;
188  m=(c-b)/two;
189 
190  if (fb == 0) {
191  root=b;
192  x_lower=b;
193  x_upper=b;
194 
195  return o2scl::success;
196  }
197  if (o2scl::o2abs(m) <= tol) {
198  root=b;
199 
200  if (b < c) {
201  x_lower=b;
202  x_upper=c;
203  } else {
204  x_lower=c;
205  x_upper=b;
206  }
207 
208  return o2scl::success;
209  }
210 
211  if (o2scl::o2abs(e) < tol || o2scl::o2abs(fa) <= o2scl::o2abs(fb)) {
212  // [GSL] Use bisection
213  d=m;
214  e=m;
215  } else {
216 
217  // [GSL] Use inverse cubic interpolation
218  fp_t p, q, r;
219  fp_t s=fb/fa;
220 
221  if (ac_equal) {
222  p=2*m*s;
223  q=1-s;
224  } else {
225  q=fa/fc;
226  r=fb/fc;
227  p=s*(2*m*q*(q-r)-(b-a)*(r-1));
228  q=(q-1)*(r-1)*(s-1);
229  }
230 
231  if (p > 0) {
232  q=-q;
233  } else {
234  p=-p;
235  }
236  fp_t dtmp;
237  if (3*m*q-o2scl::o2abs(tol*q)<o2scl::o2abs(e*q)) dtmp=3*m*q-o2scl::o2abs(tol*q);
238  else dtmp=o2scl::o2abs(e*q);
239  if (2*p<dtmp) {
240  e=d;
241  d=p/q;
242  } else {
243  // [GSL] Interpolation failed, fall back to bisection.
244  d=m;
245  e=m;
246  }
247  }
248 
249  a=b;
250  fa=fb;
251 
252  if (o2scl::o2abs(d) > tol) {
253  b+=d;
254  } else {
255  b+=(m > 0 ? +tol : -tol);
256  }
257 
258  fb=f(b);
259 
260  // Update the best estimate of the root and bounds on each
261  // iteration
262 
263  root=b;
264 
265  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
266  c=a;
267  }
268  if (b < c) {
269  x_lower=b;
270  x_upper=c;
271  } else {
272  x_lower=c;
273  x_upper=b;
274  }
275 
276  return o2scl::success;
277  }
278 
279  /// Solve \c func in region \f$ x_1<x<x_2 \f$ returning \f$ x_1 \f$.
280  virtual int solve_bkt(fp_t &x1, fp_t x2, func_t &f) {
281 
282  int status, iter=0;
283 
284  int set_ret=set(f,x1,x2);
285  if (set_ret!=0) return set_ret;
286 
287  if (test_form==0) {
288 
289  // Test the bracket size
290 
291  status=gsl_continue;
292  while (status==gsl_continue && iter<this->ntrial) {
293 
294  iter++;
295  iterate(f);
297  this->tol_abs,this->tol_rel);
298 
299  if (this->verbose>0) {
300  fp_t y;
301 
302  y=f(root);
303 
305  this->tol_abs,"root_brent_gsl (abs)");
306  }
307  }
308 
309  } else if (test_form==1) {
310 
311  // Test the residual
312 
313  status=gsl_continue;
314  while (status==gsl_continue && iter<this->ntrial) {
315 
316  iter++;
317  iterate(f);
318 
319  fp_t y=f(root);
320 
321  if (o2scl::o2abs(y)<this->tol_rel) status=o2scl::success;
322 
323  if (this->verbose>0) {
324  this->print_iter(root,y,iter,o2scl::o2abs(y),this->tol_rel,
325  "root_brent_gsl (rel)");
326  }
327  }
328 
329 
330  } else {
331 
332  // Test the bracket size and the residual
333 
334  status=gsl_continue;
335  while (status==gsl_continue && iter<this->ntrial) {
336 
337  iter++;
338  iterate(f);
340  this->tol_abs,this->tol_rel);
341 
342  if (status==o2scl::success) {
343  fp_t y=f(root);
344  if (o2scl::o2abs(y)>=this->tol_rel) status=gsl_continue;
345  if (this->verbose>0) {
346  this->print_iter(root,y,iter,o2scl::o2abs(y),this->tol_rel,
347  "root_brent_gsl (rel)");
348  }
349  } else {
350  if (this->verbose>0) {
351  fp_t y=f(root);
353  this->tol_abs,"root_brent_gsl (abs)");
354  }
355  }
356  }
357 
358  }
359 
360  x1=root;
361 
362  if (iter>=this->ntrial) {
363  O2SCL_CONV2_RET("Function root_brent_gsl::solve_bkt() exceeded ",
364  "maximum number of iterations.",o2scl::exc_emaxiter,
365  this->err_nonconv);
366  }
367 
368  if (status!=o2scl::success) {
369  return status;
370  }
371 
372  return o2scl::success;
373  }
374 
375  /// The type of convergence test applied: 0, 1, or 2 (default 0)
377 
378  /// Get the most recent value of the root
379  fp_t get_root() { return root; }
380 
381  /// Get the lower limit
382  fp_t get_lower() { return x_lower; }
383 
384  /// Get the upper limit
385  fp_t get_upper() { return x_upper; }
386 
387  /** \brief Set the information for the solver
388 
389  This function currently always returns \ref success.
390  */
391  int set(func_t &ff, fp_t lower, fp_t upper) {
392 
393  if (lower > upper) {
394  fp_t tmp=lower;
395  lower=upper;
396  upper=tmp;
397  }
398 
399  x_lower=lower;
400  x_upper=upper;
401  fp_t two=2;
402  root=(lower+upper)/two;
403 
404  fp_t f_lower, f_upper;
405 
406  f_lower=ff(x_lower);
407  f_upper=ff(x_upper);
408 
409  a=x_lower;
410  fa=f_lower;
411 
412  b=x_upper;
413  fb=f_upper;
414 
415  c=x_upper;
416  fc=f_upper;
417 
418  d=x_upper-x_lower;
419  e=x_upper-x_lower;
420 
421  if ((f_lower<0.0 && f_upper<0.0) ||
422  (f_lower>0.0 && f_upper>0.0)) {
423  O2SCL_CONV2_RET("Endpoints don't straddle y=0 in ",
424  "root_brent_gsl::set().",o2scl::exc_einval,
425  this->err_nonconv);
426  }
427 
428  return o2scl::success;
429 
430  }
431 
432 #ifndef DOXYGEN_INTERNAL
433 
434  protected:
435 
436  /// The present solution estimate
437  fp_t root;
438  /// The present lower limit
439  fp_t x_lower;
440  /// The present upper limit
441  fp_t x_upper;
442 
443  /// \name Storage for solver state
444  //@{
445  fp_t a, b, c, d, e;
446  fp_t fa, fb, fc;
447  //@}
448 
449 #endif
450 
451  };
452 
453 #ifndef DOXYGEN_NO_O2NS
454 }
455 #endif
456 
457 #endif
o2scl::root< funct, funct, double >::tol_rel
double tol_rel
The maximum value of the functions for success (default )
Definition: root.h:66
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_brent_gsl::get_upper
fp_t get_upper()
Get the upper limit.
Definition: root_brent_gsl.h:385
o2scl::root_brent_gsl::x_lower
fp_t x_lower
The present lower limit.
Definition: root_brent_gsl.h:439
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
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_brent_gsl::solve_bkt
virtual int solve_bkt(fp_t &x1, fp_t x2, func_t &f)
Solve func in region returning .
Definition: root_brent_gsl.h:280
o2scl::exc_emaxiter
@ exc_emaxiter
exceeded max number of iterations
Definition: err_hnd.h:73
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::root
One-dimensional solver [abstract base].
Definition: root.h:48
o2scl::root_brent_gsl::type
virtual const char * type()
Return the type, "root_brent_gsl".
Definition: root_brent_gsl.h:157
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::root_brent_gsl
One-dimensional root-finding (GSL)
Definition: root_brent_gsl.h:92
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::root_brent_gsl::x_upper
fp_t x_upper
The present upper limit.
Definition: root_brent_gsl.h:441
o2scl::gsl_continue
@ gsl_continue
iteration has not converged
Definition: err_hnd.h:51
o2scl::root_brent_gsl::test_form
int test_form
The type of convergence test applied: 0, 1, or 2 (default 0)
Definition: root_brent_gsl.h:376
o2scl::root_brent_gsl::iterate
int iterate(func_t &f)
Perform an iteration.
Definition: root_brent_gsl.h:163
o2scl::root_brent_gsl::test_interval
int test_interval(fp_t xx_lower, fp_t xx_upper, fp_t epsabs, fp_t epsrel)
Floating point-type agnostic version of gsl_root_test_interval() .
Definition: root_brent_gsl.h:100
o2scl::root_brent_gsl::get_root
fp_t get_root()
Get the most recent value of the root.
Definition: root_brent_gsl.h:379
o2scl::o2abs
float o2abs(const float x)
Absolute value for single precision numbers.
o2scl::root_brent_gsl::root
fp_t root
The present solution estimate.
Definition: root_brent_gsl.h:437
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::exc_ebadtol
@ exc_ebadtol
user specified an invalid tolerance
Definition: err_hnd.h:77
o2scl::root< funct, funct, double >::ntrial
int ntrial
Maximum number of iterations (default 100)
Definition: root.h:77
o2scl::root_brent_gsl::set
int set(func_t &ff, fp_t lower, fp_t upper)
Set the information for the solver.
Definition: root_brent_gsl.h:391
o2scl::root_brent_gsl::get_lower
fp_t get_lower()
Get the lower limit.
Definition: root_brent_gsl.h:382

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