Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
45namespace ROOT {
46namespace 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
77bool 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
93bool 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;
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
146double GSLRootFinder::Root() const {
147 // return cached value
148 return fRoot;
149}
150/**
151double GSLRootFinder::XLower() const {
152 return fXlow;
153}
154
155double GSLRootFinder::XUpper() const {
156 return fXup;
157}
158*/
159const char * GSLRootFinder::Name() const {
160 // get GSL name
161 return gsl_root_fsolver_name(fS->Solver() );
162}
163
164bool 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
#define MATH_INFO_MSGVAL(loc, txt, x)
Definition Error.h:101
#define MATH_ERROR_MSG(loc, str)
Definition Error.h:83
#define f(i)
Definition RSha256.hxx:104
Wrapper class to the gsl_function C structure.
void SetFunction(const FuncType &f)
fill the GSL C struct from a generic C++ callable object implementing operator()
void SetFuncPointer(GSLFuncPointer f)
set in the GSL C struct the pointer to the function evaluation
void SetParams(void *p)
set in the GSL C struct the extra-object pointer
bool IsValid()
check if function is valid (has been set)
Root-Finder implementation class using GSL.
gsl_root_fsolver * Solver() const
Base class for GSL Root-Finding algorithms for one dimensional functions which do not use function de...
GSLFunctionWrapper * fFunction
bool SetFunction(const IGenFunction &f, double xlow, double xup)
Sets the function for the rest of the algorithms.
void SetSolver(GSLRootFSolver *s)
GSLRootFinder & operator=(const GSLRootFinder &)
const char * Name() const
double GSLRootFinder::XLower() const { return fXlow; }
int Iterate()
This method is implemented only by the GSLRootFinder and GSLRootFinderDeriv classes and will return a...
double Root() const
Returns the previously calculated root.
bool Solve(int maxIter=100, double absTol=1E-8, double relTol=1E-10)
Find the root.
Interface (abstract class) for generic functions objects of one-dimension Provides a method to evalua...
Definition IFunction.h:135
Interface for finding function roots of one-dimensional functions.
Namespace for new Math classes and functions.
int TestInterval(double xlow, double xup, double epsAbs, double epsRel)
double(* GSLFuncPointer)(double, void *)
Function pointer corresponding to gsl_function signature.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...