ROOT  6.06/09
Reference Guide
GSLMultiRootFinder.h
Go to the documentation of this file.
1 // @(#)root/mathmore:$Id$
2 // Author: L. Moneta 03/2011
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 // Header file for class GSLMultiRootFinder
26 //
27 
28 #ifndef ROOT_Math_GSLMultiRootFinder
29 #define ROOT_Math_GSLMultiRootFinder
30 
31 
32 
33 #ifndef ROOT_Math_IFunction
34 #include "Math/IFunction.h"
35 #endif
36 
37 #ifndef ROOT_Math_WrappedFunction
38 #include "Math/WrappedFunction.h"
39 #endif
40 
41 #include <vector>
42 
43 #include <iostream>
44 
45 namespace ROOT {
46 namespace Math {
47 
48 
49  class GSLMultiRootBaseSolver;
50 
51  /** @defgroup MultiRoot Multidimensional ROOT finding
52  Classes for finding the roots of a multi-dimensional system.
53  @ingroup NumAlgo
54  */
55 
56  /**
57  Class for Multidimensional root finding algorithms bassed on GSL. This class is used to solve a
58  non-linear system of equations:
59 
60  f1(x1,....xn) = 0
61  f2(x1,....xn) = 0
62  ..................
63  fn(x1,....xn) = 0
64 
65  See the GSL <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html"> online manual</A> for
66  information on the GSL MultiRoot finding algorithms
67 
68  The available GSL algorithms require the derivatives of the supplied functions or not (they are
69  computed internally by GSL). In the first case the user needs to provide a list of multidimensional functions implementing the
70  gradient interface (ROOT::Math::IMultiGradFunction) while in the second case it is enough to supply a list of
71  functions impelmenting the ROOT::Math::IMultiGenFunction interface.
72  The available algorithms requiring derivatives (see also the GSL
73  <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-using-Derivatives.html">documentation</A> )
74  are the followings:
75  <ul>
76  <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridSJ</tt> with name <it>"HybridSJ"</it>: modified Powell's hybrid
77  method as implemented in HYBRJ in MINPACK
78  <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridJ</tt> with name <it>"HybridJ"</it>: unscaled version of the
79  previous algorithm</li>
80  <li><tt>ROOT::Math::GSLMultiRootFinder::kNewton</tt> with name <it>"Newton"</it>: Newton method </li>
81  <li><tt>ROOT::Math::GSLMultiRootFinder::kGNewton</tt> with name <it>"GNewton"</it>: modified Newton method </li>
82  </ul>
83  The algorithms without derivatives (see also the GSL
84  <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-without-Derivatives.html">documentation</A> )
85  are the followings:
86  <ul>
87  <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridS</tt> with name <it>"HybridS"</it>: same as HybridSJ but using
88  finate difference approximation for the derivatives</li>
89  <li><tt>ROOT::Math::GSLMultiRootFinder::kHybrid</tt> with name <it>"Hybrid"</it>: unscaled version of the
90  previous algorithm</li>
91  <li><tt>ROOT::Math::GSLMultiRootFinder::kDNewton</tt> with name <it>"DNewton"</it>: discrete Newton algorithm </li>
92  <li><tt>ROOT::Math::GSLMultiRootFinder::kBroyden</tt> with name <it>"Broyden"</it>: Broyden algorithm </li>
93  </ul>
94 
95  @ingroup MultiRoot
96  */
97 
98 
100 
101  public:
102 
103  /**
104  enumeration specifying the types of GSL multi root finders
105  requiring the derivatives
106 
107  */
108  enum EDerivType {
113  };
114  /**
115  enumeration specifying the types of GSL multi root finders
116  which do not require the derivatives
117 
118  */
119  enum EType {
124  };
125 
126 
127 
128  /// create a multi-root finder based on an algorithm not requiring function derivative
130 
131  /// create a multi-root finder based on an algorithm requiring function derivative
133 
134  /*
135  create a multi-root finder using a string.
136  The names are those defined in the GSL manuals
137  after having remived the GSL prefix (gsl_multiroot_fsolver).
138  Default algorithm is "hybrids" (without derivative).
139  */
140  GSLMultiRootFinder(const char * name = 0);
141 
142  /// destructor
143  virtual ~GSLMultiRootFinder();
144 
145  private:
146  // usually copying is non trivial, so we make this unaccessible
149 
150  public:
151 
152  /// set the type for an algorithm without derivatives
154  fType = type; fUseDerivAlgo = false;
155  }
156 
157  /// set the type of algorithm using derivatives
159  fType = type; fUseDerivAlgo = true;
160  }
161 
162  /// set the type using a string
163  void SetType(const char * name);
164 
165  /*
166  add the list of functions f1(x1,..xn),...fn(x1,...xn). The list must contain pointers of
167  ROOT::Math::IMultiGenFunctions. The method requires the
168  the begin and end of the list iterator.
169  The list can be any stl container or a simple array of ROOT::Math::IMultiGenFunctions* or
170  whatever implementing an iterator.
171  If using a derivative type algorithm the function pointers must implement the
172  ROOT::Math::IMultiGradFunction interface
173  */
174  template<class FuncIterator>
175  bool SetFunctionList( FuncIterator begin, FuncIterator end) {
176  bool ret = true;
177  for (FuncIterator itr = begin; itr != end; ++itr) {
178  const ROOT::Math::IMultiGenFunction * f = *itr;
179  ret &= AddFunction( *f);
180  }
181  return ret;
182  }
183 
184  /*
185  add (set) a single function fi(x1,...xn) which is part of the system of
186  specifying the begin and end of the iterator.
187  If using a derivative type algorithm the function must implement the
188  ROOT::Math::IMultiGradFunction interface
189  Return the current number of function in the list and 0 if failed to add the function
190  */
192 
193  /// same method as before but using any function implementing
194  /// the operator(), so can be wrapped in a IMultiGenFunction interface
195  template <class Function>
196  int AddFunction( Function & f, int ndim) {
197  // no need to care about lifetime of wfunc. It will be cloned inside AddFunction
198  WrappedMultiFunction<Function &> wfunc(f, ndim);
199  return AddFunction(wfunc);
200  }
201 
202  /**
203  return the number of sunctions set in the class.
204  The number must be equal to the dimension of the functions
205  */
206  unsigned int Dim() const { return fFunctions.size(); }
207 
208  /// clear list of functions
209  void Clear();
210 
211  /// return the root X values solving the system
212  const double * X() const;
213 
214  /// return the function values f(X) solving the system
215  /// i.e. they must be close to zero at the solution
216  const double * FVal() const;
217 
218  /// return the last step size
219  const double * Dx() const;
220 
221 
222  /**
223  Find the root starting from the point X;
224  Use the number of iteration and tolerance if given otherwise use
225  default parameter values which can be defined by
226  the static method SetDefault...
227  */
228  bool Solve(const double * x, int maxIter = 0, double absTol = 0, double relTol = 0);
229 
230  /// Return number of iterations
231  int Iterations() const {
232  return fIter;
233  }
234 
235  /// Return the status of last root finding
236  int Status() const { return fStatus; }
237 
238  /// Return the algorithm name
239  const char * Name() const;
240 
241  /*
242  set print level
243  level = 0 quiet (no messages print)
244  = 1 print only the result
245  = 3 max debug. Print result at each iteration
246  */
247  void SetPrintLevel(int level) { fPrintLevel = level; }
248 
249  /// return the print level
250  int PrintLevel() const { return fPrintLevel; }
251 
252 
253  //-- static methods to set configurations
254 
255  /// set tolerance (absolute and relative)
256  /// relative tolerance is only use to verify the convergence
257  /// do it is a minor parameter
258  static void SetDefaultTolerance(double abstol, double reltol = 0 );
259 
260  /// set maximum number of iterations
261  static void SetDefaultMaxIterations(int maxiter);
262 
263  /// print iteration state
264  void PrintState(std::ostream & os = std::cout);
265 
266 
267  protected:
268 
269  // return type given a name
270  std::pair<bool,int> GetType(const char * name);
271  // clear list of functions
272  void ClearFunctions();
273 
274 
275  private:
276 
277  int fIter; // current numer of iterations
278  int fStatus; // current status
279  int fPrintLevel; // print level
280 
281  // int fMaxIter; // max number of iterations
282  // double fAbsTolerance; // absolute tolerance
283  // double fRelTolerance; // relative tolerance
284  int fType; // type of algorithm
285  bool fUseDerivAlgo; // algorithm using derivative
286 
288  std::vector<ROOT::Math::IMultiGenFunction *> fFunctions; //! transient Vector of the functions
289 
290 
291  };
292 
293  // use typedef for most sensible name
295 
296 } // namespace Math
297 } // namespace ROOT
298 
299 
300 #endif /* ROOT_Math_GSLMultiRootFinder */
Class for Multidimensional root finding algorithms bassed on GSL.
std::vector< ROOT::Math::IMultiGenFunction * > fFunctions
const double absTol
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
const double * FVal() const
return the function values f(X) solving the system i.e.
GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders This is the base class...
void Clear()
clear list of functions
GSLMultiRootBaseSolver * fSolver
EType
enumeration specifying the types of GSL multi root finders which do not require the derivatives ...
GSLMultiRootFinder & operator=(const GSLMultiRootFinder &)
int Iterations() const
Return number of iterations.
int PrintLevel() const
return the print level
Double_t x[n]
Definition: legend1.C:17
bool SetFunctionList(FuncIterator begin, FuncIterator end)
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 ...
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 Status() const
Return the status of last root finding.
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...
unsigned int Dim() const
return the number of sunctions set in the class.
int AddFunction(const ROOT::Math::IMultiGenFunction &func)
EDerivType
enumeration specifying the types of GSL multi root finders requiring the derivatives ...
double f(double x)
GSLMultiRootFinder MultiRootFinder
int type
Definition: TGX11.cxx:120
const double * Dx() const
return the last step size
double func(double *x, double *p)
Definition: stressTF1.cxx:213
int AddFunction(Function &f, int ndim)
same method as before but using any function implementing the operator(), so can be wrapped in a IMul...
Template class to wrap any C++ callable object implementing operator() (const double * x) in a multi-...
Namespace for new Math classes and functions.
#define name(a, b)
Definition: linkTestLib0.cpp:5
void SetType(EType type)
set the type for an algorithm without derivatives
const char * Name() const
Return the algorithm name.
std::pair< bool, int > GetType(const char *name)
void SetType(EDerivType type)
set the type of algorithm using derivatives
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:63
void PrintState(std::ostream &os=std::cout)
print iteration state
const double * X() const
return the root X values solving the system
virtual ~GSLMultiRootFinder()
destructor