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