Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
GSLMultiRootFinder.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 GSLMultiRootFinder
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"
34#include "GSLMultiRootSolver.h"
35#include "Math/Error.h"
36
37#include "gsl/gsl_multiroots.h"
38#include "gsl/gsl_errno.h"
39#include <cmath>
40#include <iomanip>
41
42#include <algorithm>
43#include <functional>
44#include <ctype.h> // need to use c version of tolower defined here
45
46
47namespace ROOT {
48namespace Math {
49
50 // default values
51
52 int gDefaultMaxIter = 100;
53 double gDefaultAbsTolerance = 1.E-6;
54 double gDefaultRelTolerance = 1.E-10;
55
56// impelmentation of static methods
57void GSLMultiRootFinder::SetDefaultTolerance(double abstol, double reltol ) {
58 // set default tolerance
59 gDefaultAbsTolerance = abstol;
60 if (reltol > 0) gDefaultRelTolerance = reltol;
61}
63 // set default max iter
64 gDefaultMaxIter = maxiter;
65}
66
68 fIter(0), fStatus(-1), fPrintLevel(0),
69 fType(type), fUseDerivAlgo(false),
70 fSolver(0)
71{
72 // constructor for non derivative type
73 fFunctions.reserve(2);
74}
75
77 fIter(0), fStatus(-1), fPrintLevel(0),
78 fType(type), fUseDerivAlgo(true),
79 fSolver(0)
80{
81 // constructor for non derivative type
82 fFunctions.reserve(2);
83}
84
86 fIter(0), fStatus(-1), fPrintLevel(0),
87 fType(0), fUseDerivAlgo(false),
88 fSolver(0)
89{
90 // constructor for a string
91 fFunctions.reserve(2);
93}
94
96{
97 // delete function wrapper
99 if (fSolver) delete fSolver;
100}
101
103{
104}
105
107{
108 // dummy operator=
109 if (this == &rhs) return *this; // time saving self-test
110
111 return *this;
112}
113
115 // set type using a string
116 std::pair<bool,int> type = GetType(name);
117 fUseDerivAlgo = type.first;
118 fType = type.second;
119}
120
121
123 // add a new function in the vector
125 if (!f) return 0;
126 fFunctions.push_back(f);
127 return fFunctions.size();
128}
129
131 // clear the function list
132 for (unsigned int i = 0; i < fFunctions.size(); ++i) {
133 if (fFunctions[i] != 0 ) delete fFunctions[i];
134 fFunctions[i] = 0;
135 }
136 fFunctions.clear();
137}
138
140 // clear the function list and the solver
142 if (fSolver) Clear();
143 fSolver = 0;
144}
145
146
147const double * GSLMultiRootFinder::X() const {
148 // return x
149 return (fSolver != 0) ? fSolver->X() : 0;
150}
151const double * GSLMultiRootFinder::Dx() const {
152 // return x
153 return (fSolver != 0) ? fSolver->Dx() : 0;
154}
155const double * GSLMultiRootFinder::FVal() const {
156 // return x
157 return (fSolver != 0) ? fSolver->FVal() : 0;
158}
159const char * GSLMultiRootFinder::Name() const {
160 // get GSL name
161 return (fSolver != 0) ? fSolver->Name().c_str() : "";
162}
163
164// bool GSLMultiRootFinder::AddFunction( const ROOT::Math::IMultiGenFunction & func) {
165// // clone and add function to the list
166// // If using a derivative algorithm the function is checked if it implements
167// // the gradient interface. If this is not the case the type is set to non-derivatibe algo
168// ROOT::Math::IGenMultiFunction * f = func.Clone();
169// if (f != 0) return false;
170// if (fUseDerivAlgo) {
171// bool gradFunc = (dynamic_cast<ROOT::Math::IMultiGradFunction *> (f) != 0 );
172// if (!gradFunc) {
173// MATH_ERROR_MSG("GSLMultiRootFinder::AddFunction","Function does not provide gradient interface");
174// MATH_WARN_MSG("GSLMultiRootFinder::AddFunction","clear the function list");
175// ClearFunctions();
176// return false;
177// }
178// }
179// fFunctions.push_back(f);
180// return true;
181// }
182
183 const gsl_multiroot_fsolver_type * GetGSLType(GSLMultiRootFinder::EType type) {
184 //helper functions to find GSL type
185 switch(type)
186 {
188 return gsl_multiroot_fsolver_hybrids;
190 return gsl_multiroot_fsolver_hybrid;
192 return gsl_multiroot_fsolver_dnewton;
194 return gsl_multiroot_fsolver_broyden;
195 default:
196 return gsl_multiroot_fsolver_hybrids;
197 }
198 return 0;
199}
200
201const gsl_multiroot_fdfsolver_type * GetGSLDerivType(GSLMultiRootFinder::EDerivType type) {
202//helper functions to find GSL deriv type
203 switch(type)
204 {
206 return gsl_multiroot_fdfsolver_hybridsj;
208 return gsl_multiroot_fdfsolver_hybridj;
210 return gsl_multiroot_fdfsolver_newton;
212 return gsl_multiroot_fdfsolver_gnewton;
213 default:
214 return gsl_multiroot_fdfsolver_hybridsj;
215 }
216 return 0; // cannot happen
217}
218
219std::pair<bool,int> GSLMultiRootFinder::GetType(const char * name) {
220 if (name == 0) return std::make_pair<bool,int>(false, -1);
221 std::string aname = name;
222 std::transform(aname.begin(), aname.end(), aname.begin(), (int(*)(int)) tolower );
223
224 if (aname.find("hybridsj") != std::string::npos) return std::make_pair(true, kHybridSJ);
225 if (aname.find("hybridj") != std::string::npos) return std::make_pair(true, kHybridJ);
226 if (aname.find("hybrids") != std::string::npos) return std::make_pair(false, kHybridS);
227 if (aname.find("hybrid") != std::string::npos) return std::make_pair(false, kHybrid);
228 if (aname.find("gnewton") != std::string::npos) return std::make_pair(true, kGNewton);
229 if (aname.find("dnewton") != std::string::npos) return std::make_pair(false, kDNewton);
230 if (aname.find("newton") != std::string::npos) return std::make_pair(true, kNewton);
231 if (aname.find("broyden") != std::string::npos) return std::make_pair(false, kBroyden);
232 MATH_INFO_MSG("GSLMultiRootFinder::GetType","Unknow algorithm - use default one");
233 return std::make_pair(false, -1);
234}
235
236bool GSLMultiRootFinder::Solve (const double * x, int maxIter, double absTol, double relTol)
237{
238 fIter = 0;
239 // create the solvers - delete previous existing solver
240 if (fSolver) delete fSolver;
241 fSolver = 0;
242
243 if (fFunctions.size() == 0) {
244 MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Function list is empty");
245 fStatus = -1;
246 return false;
247 }
248
249 if (fUseDerivAlgo) {
252 }
253 else {
254 EType type = (EType) fType;
256 }
257
258
259 // first set initial values and function
260 assert(fSolver != 0);
261 bool ret = fSolver->InitSolver( fFunctions, x);
262 if (!ret) {
263 MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Error initializing the solver");
264 fStatus = -2;
265 return false;
266 }
267
268 if (maxIter == 0) maxIter = gDefaultMaxIter;
269 if (absTol <= 0) absTol = gDefaultAbsTolerance;
270 if (relTol <= 0) relTol = gDefaultRelTolerance;
271
272 if (fPrintLevel >= 1)
273 std::cout << "GSLMultiRootFinder::Solve:" << Name() << " max iterations " << maxIter << " and tolerance " << absTol << std::endl;
274
275 // find the roots by iterating
276 fStatus = 0;
277 int status = 0;
278 int iter = 0;
279 do {
280 iter++;
281 status = fSolver->Iterate();
282
283 if (fPrintLevel >= 2) {
284 std::cout << "GSLMultiRootFinder::Solve - iteration # " << iter << " status = " << status << std::endl;
285 PrintState();
286 }
287 // act in case of error
288 if (status == GSL_EBADFUNC) {
289 MATH_ERROR_MSG("GSLMultiRootFinder::Solve","The iteration encountered a singolar point due to a bad function value");
290 fStatus = status;
291 break;
292 }
293 if (status == GSL_ENOPROG) {
294 MATH_ERROR_MSG("GSLMultiRootFinder::Solve","The iteration is not making any progress");
295 fStatus = status;
296 break;
297 }
298 if (status != GSL_SUCCESS) {
299 MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Uknown iteration error - exit");
300 fStatus = status;
301 break;
302 }
303
304 // test also residual
305 status = fSolver->TestResidual(absTol);
306
307
308 // should test also the Delta ??
309 int status2 = fSolver->TestDelta(absTol, relTol);
310 if (status2 == GSL_SUCCESS) {
311 MATH_INFO_MSG("GSLMultiRootFinder::Solve","The iteration converged");
312 }
313 }
314 while (status == GSL_CONTINUE && iter < maxIter);
315 if (status == GSL_CONTINUE) {
316 MATH_INFO_MSGVAL("GSLMultiRootFinder::Solve","exceeded max iterations, reached tolerance is not sufficient",absTol);
317 }
318 if (status == GSL_SUCCESS) {
319 if (fPrintLevel>=1) { // print the result
320 MATH_INFO_MSG("GSLMultiRootFinder::Solve","The iteration converged");
321 std::cout << "GSL Algorithm used is : " << fSolver->Name() << std::endl;
322 std::cout << "Number of iterations = " << iter<< std::endl;
323
324 PrintState();
325 }
326 }
327 fIter = iter;
328 fStatus = status;
329 return (fStatus == GSL_SUCCESS);
330
331}
332
333void GSLMultiRootFinder::PrintState(std::ostream & os) {
334 // print current state
335 if (!fSolver) return;
336 double ndigits = std::log10( double( Dim() ) );
337 int wi = int(ndigits)+1;
338 const double * xtmp = fSolver->X();
339 const double * ftmp = fSolver->FVal();
340 os << "Root values = ";
341 for (unsigned int i = 0; i< Dim(); ++i)
342 os << "x[" << std::setw(wi) << i << "] = " << std::setw(12) << xtmp[i] << " ";
343 os << std::endl;
344 os << "Function values = ";
345 for (unsigned int i = 0; i< Dim(); ++i)
346 os << "f[" << std::setw(wi) << i << "] = " << std::setw(12) << ftmp[i] << " ";
347 os << std::endl;
348}
349
350
351
352} // namespace Math
353} // namespace ROOT
#define MATH_INFO_MSG(loc, str)
Pre-processor macro to report messages which can be configured to use ROOT error or simply an std::io...
Definition Error.h:77
#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
char name[80]
Definition TGX11.cxx:110
int type
Definition TGX11.cxx:121
int TestResidual(double absTol) const
test using abs tolerance Sum |f|_i < absTol
const double * FVal() const
return function values
const double * X() const
solution values at the current iteration
virtual const std::string & Name() const =0
return name
const double * Dx() const
return function steps
bool InitSolver(const std::vector< ROOT::Math::IMultiGenFunction * > &funcVec, const double *x)
init the solver with function list and initial values
virtual int Iterate()=0
perform an iteration
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component
GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders using derivatives.
Class for Multidimensional root finding algorithms bassed on GSL.
GSLMultiRootFinder & operator=(const GSLMultiRootFinder &)
unsigned int Dim() const
return the number of sunctions set in the class.
const double * Dx() const
return the last step size
const double * FVal() const
return the function values f(X) solving the system i.e.
void SetType(EType type)
set the type for an algorithm without derivatives
std::vector< ROOT::Math::IMultiGenFunction * > fFunctions
bool Solve(const double *x, int maxIter=0, double absTol=0, double relTol=0)
Find the root starting from the point X; Use the number of iteration and tolerance if given otherwise...
EType
enumeration specifying the types of GSL multi root finders which do not require the derivatives
std::pair< bool, int > GetType(const char *name)
const char * Name() const
Return the algorithm name used for solving Note the name is available only after having called solved...
void PrintState(std::ostream &os=std::cout)
print iteration state
GSLMultiRootBaseSolver * fSolver
EDerivType
enumeration specifying the types of GSL multi root finders requiring the derivatives
GSLMultiRootFinder(EType type)
create a multi-root finder based on an algorithm not requiring function derivative
void Clear()
clear list of functions
static void SetDefaultTolerance(double abstol, double reltol=0)
set tolerance (absolute and relative) relative tolerance is only use to verify the convergence do it ...
int AddFunction(const ROOT::Math::IMultiGenFunction &func)
const double * X() const
return the root X values solving the system
static void SetDefaultMaxIterations(int maxiter)
set maximum number of iterations
GSLMultiRootSolver, internal class for implementing GSL multi-root finders not using derivatives.
Documentation for the abstract class IBaseFunctionMultiDim.
Definition IFunction.h:62
virtual IBaseFunctionMultiDimTempl< T > * Clone() const =0
Clone a function.
Double_t x[n]
Definition legend1.C:17
Namespace for new Math classes and functions.
const gsl_multiroot_fsolver_type * GetGSLType(GSLMultiRootFinder::EType type)
double gDefaultRelTolerance
const gsl_multiroot_fdfsolver_type * GetGSLDerivType(GSLMultiRootFinder::EDerivType type)
double gDefaultAbsTolerance
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...