root_brent_gsl.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 /* 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=funct11> class root_brent_gsl :
93  public root_bkt<func_t> {
94 
95  public:
96 
97  root_brent_gsl() {
98  test_form=0;
99  }
100 
101  /// Return the type, \c "root_brent_gsl".
102  virtual const char *type() { return "root_brent_gsl"; }
103 
104  /** \brief Perform an iteration
105 
106  This function currently always returns \ref success.
107  */
108  int iterate(func_t &f) {
109 
110  double tol, m;
111 
112  int ac_equal=0;
113 
114  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
115  ac_equal=1;
116  c=a;
117  fc=fa;
118  d=b-a;
119  e=b-a;
120  }
121 
122  if (fabs(fc) < fabs(fb)) {
123  ac_equal=1;
124  a=b;
125  b=c;
126  c=a;
127  fa=fb;
128  fb=fc;
129  fc=fa;
130  }
131 
132  tol=0.5*fabs(b)*std::numeric_limits<double>::epsilon();
133  m=0.5*(c-b);
134 
135  if (fb == 0) {
136  root=b;
137  x_lower=b;
138  x_upper=b;
139 
140  return o2scl::success;
141  }
142  if (fabs(m) <= tol) {
143  root=b;
144 
145  if (b < c) {
146  x_lower=b;
147  x_upper=c;
148  } else {
149  x_lower=c;
150  x_upper=b;
151  }
152 
153  return o2scl::success;
154  }
155 
156  if (fabs(e) < tol || fabs(fa) <= fabs(fb)) {
157  // [GSL] Use bisection
158  d=m;
159  e=m;
160  } else {
161 
162  // [GSL] Use inverse cubic interpolation
163  double p, q, r;
164  double s=fb/fa;
165 
166  if (ac_equal) {
167  p=2*m*s;
168  q=1-s;
169  } else {
170  q=fa/fc;
171  r=fb/fc;
172  p=s*(2*m*q*(q-r)-(b-a)*(r-1));
173  q=(q-1)*(r-1)*(s-1);
174  }
175 
176  if (p > 0) {
177  q=-q;
178  } else {
179  p=-p;
180  }
181  double dtmp;
182  if (3*m*q-fabs(tol*q)<fabs(e*q)) dtmp=3*m*q-fabs(tol*q);
183  else dtmp=fabs(e*q);
184  if (2*p<dtmp) {
185  e=d;
186  d=p/q;
187  } else {
188  // [GSL] Interpolation failed, fall back to bisection.
189  d=m;
190  e=m;
191  }
192  }
193 
194  a=b;
195  fa=fb;
196 
197  if (fabs(d) > tol) {
198  b+=d;
199  } else {
200  b+=(m > 0 ? +tol : -tol);
201  }
202 
203  fb=f(b);
204 
205  // Update the best estimate of the root and bounds on each
206  // iteration
207 
208  root=b;
209 
210  if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0)) {
211  c=a;
212  }
213  if (b < c) {
214  x_lower=b;
215  x_upper=c;
216  } else {
217  x_lower=c;
218  x_upper=b;
219  }
220 
221  return o2scl::success;
222  }
223 
224  /// Solve \c func in region \f$ x_1<x<x_2 \f$ returning \f$ x_1 \f$.
225  virtual int solve_bkt(double &x1, double x2, func_t &f) {
226 
227  int status, iter=0;
228 
229  if (set(f,x1,x2)!=0) {
230  O2SCL_ERR2("Function set() failed in",
231  "root_brent_gsl::solve_bkt().",exc_einval);
232  }
233 
234  if (test_form==0) {
235 
236  // Test the bracket size
237 
238  status=gsl_continue;
239  while (status==gsl_continue && iter<this->ntrial) {
240 
241  iter++;
242  iterate(f);
243  status=gsl_root_test_interval(x_lower,x_upper,
244  this->tol_abs,this->tol_rel);
245 
246  if (this->verbose>0) {
247  double y;
248 
249  y=f(root);
250 
251  this->print_iter(root,y,iter,fabs(x_upper-x_lower),this->tol_abs,
252  "root_brent_gsl (abs)");
253  }
254  }
255 
256  } else if (test_form==1) {
257 
258  // Test the residual
259 
260  status=gsl_continue;
261  while (status==gsl_continue && iter<this->ntrial) {
262 
263  iter++;
264  iterate(f);
265 
266  double y=f(root);
267 
268  if (fabs(y)<this->tol_rel) status=o2scl::success;
269 
270  if (this->verbose>0) {
271  this->print_iter(root,y,iter,fabs(y),this->tol_rel,
272  "root_brent_gsl (rel)");
273  }
274  }
275 
276 
277  } else {
278 
279  // Test the bracket size and the residual
280 
281  status=gsl_continue;
282  while (status==gsl_continue && iter<this->ntrial) {
283 
284  iter++;
285  iterate(f);
286  status=gsl_root_test_interval(x_lower,x_upper,
287  this->tol_abs,this->tol_rel);
288 
289  if (status==o2scl::success) {
290  double y=f(root);
291  if (fabs(y)>=this->tol_rel) status=gsl_continue;
292  if (this->verbose>0) {
293  this->print_iter(root,y,iter,fabs(y),this->tol_rel,
294  "root_brent_gsl (rel)");
295  }
296  } else {
297  if (this->verbose>0) {
298  double y=f(root);
299  this->print_iter(root,y,iter,fabs(x_upper-x_lower),this->tol_abs,
300  "root_brent_gsl (abs)");
301  }
302  }
303  }
304 
305  }
306 
307  x1=root;
308 
309  if (iter>=this->ntrial) {
310  O2SCL_CONV2_RET("Function root_brent_gsl::solve_bkt() exceeded ",
311  "maximum number of iterations.",o2scl::exc_emaxiter,
312  this->err_nonconv);
313  }
314 
315  if (status!=o2scl::success) {
316  return status;
317  }
318 
319  return o2scl::success;
320  }
321 
322  /// The type of convergence test applied: 0, 1, or 2 (default 0)
324 
325  /// Get the most recent value of the root
326  double get_root() { return root; }
327 
328  /// Get the lower limit
329  double get_lower() { return x_lower; }
330 
331  /// Get the upper limit
332  double get_upper() { return x_upper; }
333 
334  /** \brief Set the information for the solver
335 
336  This function currently always returns \ref success.
337  */
338  int set(func_t &ff, double lower, double upper) {
339 
340  if (lower > upper) {
341  double tmp=lower;
342  lower=upper;
343  upper=tmp;
344  }
345 
346  x_lower=lower;
347  x_upper=upper;
348  root=0.5*(lower+upper);
349 
350  double f_lower, f_upper;
351 
352  f_lower=ff(x_lower);
353  f_upper=ff(x_upper);
354 
355  a=x_lower;
356  fa=f_lower;
357 
358  b=x_upper;
359  fb=f_upper;
360 
361  c=x_upper;
362  fc=f_upper;
363 
364  d=x_upper-x_lower;
365  e=x_upper-x_lower;
366 
367  if ((f_lower<0.0 && f_upper<0.0) ||
368  (f_lower>0.0 && f_upper>0.0)) {
369  O2SCL_ERR2("Endpoints don't straddle y=0 in ",
370  "root_brent_gsl::set().",o2scl::exc_einval);
371  }
372 
373  return o2scl::success;
374 
375  }
376 
377 #ifndef DOXYGEN_INTERNAL
378 
379  protected:
380 
381  /// The present solution estimate
382  double root;
383  /// The present lower limit
384  double x_lower;
385  /// The present upper limit
386  double x_upper;
387 
388  /// \name Storage for solver state
389  //@{
390  double a, b, c, d, e;
391  double fa, fb, fc;
392  //@}
393 
394 #endif
395 
396  };
397 
398 #ifndef DOXYGEN_NO_O2NS
399 }
400 #endif
401 
402 #endif
double get_lower()
Get the lower limit.
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
double x_upper
The present upper limit.
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
double root
The present solution estimate.
iteration has not converged
Definition: err_hnd.h:51
double get_upper()
Get the upper limit.
double get_root()
Get the most recent value of the root.
One-dimensional bracketing solver [abstract base].
Definition: root.h:145
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_brent_gsl".
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
int iterate(func_t &f)
Perform an iteration.
#define O2SCL_ERR2(d, d2, n)
Set an error, two-string version.
Definition: err_hnd.h:281
virtual int solve_bkt(double &x1, double x2, func_t &f)
Solve func in region returning .
int test_form
The type of convergence test applied: 0, 1, or 2 (default 0)
One-dimensional root-finding (GSL)
double x_lower
The present lower limit.
static const double x2[5]
Definition: inte_qng_gsl.h:66
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).