Logo ROOT   6.13/01
Reference Guide
GSLRootFinder.cxx
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Authors: L. Moneta, A. Zsenei 08/2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
7  * *
8  * This library is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License *
10  * as published by the Free Software Foundation; either version 2 *
11  * of the License, or (at your option) any later version. *
12  * *
13  * This library 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 GNU *
16  * General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this library (see file COPYING); if not, write *
20  * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
21  * 330, Boston, MA 02111-1307 USA, or contact the author. *
22  * *
23  **********************************************************************/
24 
25 // Implementation file for class GSLRootFinder
26 //
27 // Created by: moneta at Sun Nov 14 11:27:11 2004
28 //
29 // Last update: Sun Nov 14 11:27:11 2004
30 //
31 
32 #include "Math/IFunction.h"
33 #include "Math/GSLRootFinder.h"
34 #include "Math/GSLRootHelper.h"
35 #include "Math/Error.h"
36 
37 #include "GSLRootFSolver.h"
38 #include "GSLFunctionWrapper.h"
39 
40 #include "gsl/gsl_roots.h"
41 #include "gsl/gsl_errno.h"
42 #include <cmath>
43 
44 
45 namespace ROOT {
46 namespace Math {
47 
48 
50  fFunction(0), fS(0),
51  fRoot(0), fXlow(0), fXup(0),
52  fIter(0), fStatus(-1),
53  fValidInterval(false)
54 {
55  // create function wrapper
57 }
58 
60 {
61  // delete function wrapper
62  if (fFunction) delete fFunction;
63 }
64 
66 {
67 }
68 
70 {
71  // dummy operator=
72  if (this == &rhs) return *this; // time saving self-test
73 
74  return *this;
75 }
76 
77 bool GSLRootFinder::SetFunction( GSLFuncPointer f, void * p, double xlow, double xup) {
78  // set from GSL function
79  fXlow = xlow;
80  fXup = xup;
82  fFunction->SetParams( p );
83 
84  int status = gsl_root_fsolver_set( fS->Solver(), fFunction->GetFunc(), xlow, xup);
85  if (status == GSL_SUCCESS)
86  fValidInterval = true;
87  else
88  fValidInterval = false;
89 
90  return fValidInterval;
91 }
92 
93 bool GSLRootFinder::SetFunction( const IGenFunction & f, double xlow, double xup) {
94  // set from IGenFunction
95  fStatus = -1; // invalid the status
96  fXlow = xlow;
97  fXup = xup;
98  fFunction->SetFunction( f );
99  int status = gsl_root_fsolver_set( fS->Solver(), fFunction->GetFunc(), xlow, xup);
100  if (status == GSL_SUCCESS)
101  fValidInterval = true;
102  else
103  fValidInterval = false;
104 
105  return fValidInterval;
106 }
107 
109  // set type of solver
110  fS = s;
111 }
112 
114  // free resources
115  if (fS) delete fS;
116 }
117 
119  // iterate
120  int status = 0;
121  if (!fFunction->IsValid() ) {
122  MATH_ERROR_MSG("GSLRootFinder::Iterate"," Function is not valid");
123  status = -1;
124  return status;
125  }
126  if (!fValidInterval ) {
127  MATH_ERROR_MSG("GSLRootFinder::Iterate"," Interval is not valid");
128  status = -2;
129  return status;
130  }
131 
132  status = gsl_root_fsolver_iterate(fS->Solver());
133 
134  // update Root
135  fRoot = gsl_root_fsolver_root(fS->Solver() );
136  // update interval
137  fXlow = gsl_root_fsolver_x_lower(fS->Solver() );
138  fXup = gsl_root_fsolver_x_upper(fS->Solver() );
139 
140  //std::cout << "iterate .." << fRoot << " status " << status << " interval "
141  // << fXlow << " " << fXup << std::endl;
142 
143  return status;
144 }
145 
146 double GSLRootFinder::Root() const {
147  // return cached value
148  return fRoot;
149 }
150 /**
151 double GSLRootFinder::XLower() const {
152  return fXlow;
153 }
154 
155 double GSLRootFinder::XUpper() const {
156  return fXup;
157 }
158 */
159 const char * GSLRootFinder::Name() const {
160  // get GSL name
161  return gsl_root_fsolver_name(fS->Solver() );
162 }
163 
164 bool GSLRootFinder::Solve (int maxIter, double absTol, double relTol)
165 {
166  // find the roots by iterating
167  fStatus = -1;
168  int status = 0;
169  int iter = 0;
170  do {
171  iter++;
172  status = Iterate();
173  //std::cerr << "RF: iteration " << iter << " status = " << status << std::endl;
174  if (status != GSL_SUCCESS) {
175  MATH_ERROR_MSG("GSLRootFinder::Solve","error returned when performing an iteration");
176  fStatus = status;
177  return false;
178  }
179  status = GSLRootHelper::TestInterval(fXlow, fXup, absTol, relTol);
180  if (status == GSL_SUCCESS) {
181  fIter = iter;
182  fStatus = status;
183  return true;
184  }
185  }
186  while (status == GSL_CONTINUE && iter < maxIter);
187  if (status == GSL_CONTINUE) {
188  double tol = std::abs(fXup-fXlow);
189  MATH_INFO_MSGVAL("GSLRootFinder::Solve","exceeded max iterations, reached tolerance is not sufficient",tol);
190  }
191  fStatus = status;
192  return false;
193 }
194 
195 
196 
197 
198 } // namespace Math
199 } // namespace ROOT
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition: IFunction.h:135
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
const char * Name() const
double GSLRootFinder::XLower() const { return fXlow; }
Root-Finder implementation class using GSL.
GSLRootFinder & operator=(const GSLRootFinder &)
Base class for GSL Root-Finding algorithms for one dimensional functions which do not use function de...
Definition: GSLRootFinder.h:73
void SetFunction(const FuncType &f)
fill the GSL C struct from a generic C++ callable object implementing operator()
double(* GSLFuncPointer)(double, void *)
Definition: GSLRootFinder.h:96
Interface for finding function roots of one-dimensional functions.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Find the root.
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
int Iterate()
This method is implemented only by the GSLRootFinder and GSLRootFinderDeriv classes and will return a...
int TestInterval(double xlow, double xup, double epsAbs, double epsRel)
GSLFunctionWrapper * fFunction
bool IsValid()
check if function is valid (has been set)
double Root() const
Returns the previously calculated root.
void SetFuncPointer(GSLFuncPointer f)
set in the GSL C struct the pointer to the function evaluation
Namespace for new Math classes and functions.
void SetSolver(GSLRootFSolver *s)
Wrapper class to the gsl_function C structure.
void SetParams(void *p)
set in the GSL C struct the extra-object pointer
gsl_root_fsolver * Solver() const
bool SetFunction(const IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
#define MATH_INFO_MSGVAL(loc, str, x)
Definition: Error.h:65