Logo ROOT   6.10/09
Reference Guide
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 
47 namespace ROOT {
48 namespace 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
57 void 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);
92  SetType(name);
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 
114 void GSLMultiRootFinder::SetType(const char * name) {
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
141  ClearFunctions();
142  if (fSolver) Clear();
143  fSolver = 0;
144 }
145 
146 
147 const double * GSLMultiRootFinder::X() const {
148  // return x
149  return (fSolver != 0) ? fSolver->X() : 0;
150 }
151 const double * GSLMultiRootFinder::Dx() const {
152  // return x
153  return (fSolver != 0) ? fSolver->Dx() : 0;
154 }
155 const double * GSLMultiRootFinder::FVal() const {
156  // return x
157  return (fSolver != 0) ? fSolver->FVal() : 0;
158 }
159 const 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 
201 const 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 
219 std::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 
236 bool 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;
255  if (!fSolver) fSolver = new GSLMultiRootSolver( GetGSLType(type), Dim() );
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 
333 void 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
Class for Multidimensional root finding algorithms bassed on GSL.
std::vector< ROOT::Math::IMultiGenFunction * > fFunctions
double gDefaultAbsTolerance
const double absTol
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
EType
enumeration specifying the types of GSL multi root finders which do not require the derivatives ...
void Clear()
clear list of functions
GSLMultiRootBaseSolver * fSolver
virtual IBaseFunctionMultiDimTempl< T > * Clone() const =0
Clone a function.
const char * Name() const
Return the algorithm name used for solving Note the name is available only after having called solved...
GSLMultiRootFinder & operator=(const GSLMultiRootFinder &)
int TestResidual(double absTol) const
test using abs tolerance Sum |f|_i < absTol
unsigned int Dim() const
return the number of sunctions set in the class.
bool InitSolver(const std::vector< ROOT::Math::IMultiGenFunction *> &funcVec, const double *x)
init the solver with function list and initial values
int TestDelta(double absTol, double relTol) const
test using abs and relative tolerance |dx| < absTol + relTol*|x| for every component ...
Double_t x[n]
Definition: legend1.C:17
double log10(double)
const double * FVal() const
return function values
virtual int Iterate()=0
perform an iteration
const double * X() const
return the root X values solving the system
double gDefaultRelTolerance
#define MATH_ERROR_MSG(loc, str)
Definition: Error.h:50
virtual const std::string & Name() const =0
return name
static void SetDefaultMaxIterations(int maxiter)
set maximum number of iterations
GSLMultiRootFinder(EType type)
create a multi-root finder based on an algorithm not requiring function derivative ...
const double * Dx() const
return function steps
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
static void SetDefaultTolerance(double abstol, double reltol=0)
set tolerance (absolute and relative) relative tolerance is only use to verify the convergence do it ...
#define MATH_INFO_MSG(loc, str)
Definition: Error.h:44
GSLMultiRootSolver, internal class for implementing GSL multi-root finders not using derivatives...
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...
const gsl_multiroot_fdfsolver_type * GetGSLDerivType(GSLMultiRootFinder::EDerivType type)
const double * FVal() const
return the function values f(X) solving the system i.e.
EDerivType
enumeration specifying the types of GSL multi root finders requiring the derivatives ...
int AddFunction(const ROOT::Math::IMultiGenFunction &func)
const double * Dx() const
return the last step size
double f(double x)
int type
Definition: TGX11.cxx:120
double func(double *x, double *p)
Definition: stressTF1.cxx:213
Namespace for new Math classes and functions.
void SetType(EType type)
set the type for an algorithm without derivatives
const gsl_multiroot_fsolver_type * GetGSLType(GSLMultiRootFinder::EType type)
GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders using derivatives...
std::pair< bool, int > GetType(const char *name)
const double * X() const
solution values at the current iteration
void PrintState(std::ostream &os=std::cout)
print iteration state
virtual ~GSLMultiRootFinder()
destructor
#define MATH_INFO_MSGVAL(loc, str, x)
Definition: Error.h:65