diff --git a/math/doc/v530/index.html b/math/doc/v530/index.html
index f38eebc0b14382d76d0139098cb37e0778a8082b..7fff654d92d5a1d97964095467bc097d9ffe67ec 100644
--- a/math/doc/v530/index.html
+++ b/math/doc/v530/index.html
@@ -3,6 +3,8 @@
<a name="math"></a>
<h3>Math Libraries</h3>
+<h4>MathCore</h4>
+
<tt>TMath</tt>
<ul>
<li>Add some new functions based on std::numeric_limits:
@@ -22,3 +24,51 @@
</ul>
</li>
</ul>
+
+<p>
+
+<h4>MathMore</h4>
+
+<ul>
+ <li><b>New class <tt>ROOT::Math::GSLMultiRootFinder</tt></b> for finding the root of system of functions.
+ The class is based on the GSL multi-root algorithm
+ (see the GSL <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html"> online
+ manual</A>) and it is used to solve a non-linear system of equations:
+ <pre>
+ f1(x1,....xn) = 0
+ f2(x1,....xn) = 0
+ ..................
+ fn(x1,....xn) = 0
+</pre>
+ The available GSL algorithms require the derivatives of the supplied functions or not (they are
+ computed internally by GSL). In the first case the user needs to provide a list of multidimensional functions implementing the
+ gradient interface (ROOT::Math::IMultiGradFunction) while in the second case it is enough to supply a list of
+ functions implementing the ROOT::Math::IMultiGenFunction interface.
+ The available algorithms requiring derivatives (see also the GSL
+ <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-using-Derivatives.html">documentation</A> )
+ are the followings:
+ <ul>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridSJ</tt> with name <it>"HybridSJ"</it>: modified Powell's hybrid
+ method as implemented in HYBRJ in MINPACK
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridJ</tt> with name <it>"HybridJ"</it>: unscaled version of the
+ previous algorithm</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kNewton</tt> with name <it>"Newton"</it>: Newton method </li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kGNewton</tt> with name <it>"GNewton"</it>: modified Newton method </li>
+ </ul>
+ The algorithms without derivatives (see also the GSL
+ <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-without-Derivatives.html">documentation</A> )
+ are the followings:
+ <ul>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridS</tt> with name <it>"HybridS"</it>: same as HybridSJ but using
+ finate difference approximation for the derivatives</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybrid</tt> with name <it>"Hybrid"</it>: unscaled version of the
+ previous algorithm</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kDNewton</tt> with name <it>"DNewton"</it>: discrete Newton algorithm </li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kBroyden</tt> with name <it>"Broyden"</it>: Broyden algorithm </li>
+ </ul>
+
+ </li>
+
+ </ul>
+
+
\ No newline at end of file
diff --git a/math/mathmore/Module.mk b/math/mathmore/Module.mk
index ddbac8efe47039ed2fb62b54445c8071eb6e5c0e..4fe8eb16aba94ebc6d00087e040704ae2caaded4 100644
--- a/math/mathmore/Module.mk
+++ b/math/mathmore/Module.mk
@@ -46,6 +46,7 @@ MATHMOREDH1 := $(MODDIRI)/Math/DistFuncMathMore.h \
$(MODDIRI)/Math/GSLMinimizer.h \
$(MODDIRI)/Math/GSLNLSMinimizer.h \
$(MODDIRI)/Math/GSLSimAnMinimizer.h \
+ $(MODDIRI)/Math/GSLMultiRootFinder.h \
$(MODDIRI)/Math/Vavilov.h \
$(MODDIRI)/Math/VavilovAccurate.h \
$(MODDIRI)/Math/VavilovAccuratePdf.h \
diff --git a/math/mathmore/inc/Math/GSLMultiRootFinder.h b/math/mathmore/inc/Math/GSLMultiRootFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..ccfe93db2d3c2f4b124e4f5c504e765e8df5fa72
--- /dev/null
+++ b/math/mathmore/inc/Math/GSLMultiRootFinder.h
@@ -0,0 +1,291 @@
+// @(#)root/mathmore:$Id$
+// Author: L. Moneta 03/2011
+
+ /**********************************************************************
+ * *
+ * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Header file for class GSLMultiRootFinder
+//
+
+#ifndef ROOT_Math_GSLMultiRootFinder
+#define ROOT_Math_GSLMultiRootFinder
+
+
+
+#ifndef ROOT_Math_IFunction
+#include "Math/IFunction.h"
+#endif
+
+#ifndef ROOT_Math_WrappedFunction
+#include "Math/WrappedFunction.h"
+#endif
+
+#include <vector>
+
+#include <iostream>
+
+namespace ROOT {
+namespace Math {
+
+
+ class GSLMultiRootBaseSolver;
+
+
+
+//________________________________________________________________________________________________________
+ /**
+ Class for Multidimensional root finding algorithms bassed on GSL. This class is used to solve a
+ non-linear system of equations:
+
+ f1(x1,....xn) = 0
+ f2(x1,....xn) = 0
+ ..................
+ fn(x1,....xn) = 0
+
+ See the GSL <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html"> online manual</A> for
+ information on the GSL MultiRoot finding algorithms
+
+ The available GSL algorithms require the derivatives of the supplied functions or not (they are
+ computed internally by GSL). In the first case the user needs to provide a list of multidimensional functions implementing the
+ gradient interface (ROOT::Math::IMultiGradFunction) while in the second case it is enough to supply a list of
+ functions impelmenting the ROOT::Math::IMultiGenFunction interface.
+ The available algorithms requiring derivatives (see also the GSL
+ <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-using-Derivatives.html">documentation</A> )
+ are the followings:
+ <ul>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridSJ</tt> with name <it>"HybridSJ"</it>: modified Powell's hybrid
+ method as implemented in HYBRJ in MINPACK
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridJ</tt> with name <it>"HybridJ"</it>: unscaled version of the
+ previous algorithm</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kNewton</tt> with name <it>"Newton"</it>: Newton method </li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kGNewton</tt> with name <it>"GNewton"</it>: modified Newton method </li>
+ </ul>
+ The algorithms without derivatives (see also the GSL
+ <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Algorithms-without-Derivatives.html">documentation</A> )
+ are the followings:
+ <ul>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybridS</tt> with name <it>"HybridS"</it>: same as HybridSJ but using
+ finate difference approximation for the derivatives</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kHybrid</tt> with name <it>"Hybrid"</it>: unscaled version of the
+ previous algorithm</li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kDNewton</tt> with name <it>"DNewton"</it>: discrete Newton algorithm </li>
+ <li><tt>ROOT::Math::GSLMultiRootFinder::kBroyden</tt> with name <it>"Broyden"</it>: Broyden algorithm </li>
+ </ul>
+
+ @ingroup MultiRoot
+ */
+
+
+ class GSLMultiRootFinder {
+
+ public:
+
+ /**
+ enumeration specifying the types of GSL multi root finders
+ requiring the derivatives
+ @ingroup MultiRoot
+ */
+ enum EDerivType {
+ kHybridSJ,
+ kHybridJ,
+ kNewton,
+ kGNewton
+ };
+ /**
+ enumeration specifying the types of GSL multi root finders
+ which do not require the derivatives
+ @ingroup MultiRoot
+ */
+ enum EType {
+ kHybridS,
+ kHybrid,
+ kDNewton,
+ kBroyden
+ };
+
+
+
+ /// create a multi-root finder based on an algorithm not requiring function derivative
+ GSLMultiRootFinder(EType type);
+
+ /// create a multi-root finder based on an algorithm requiring function derivative
+ GSLMultiRootFinder(EDerivType type);
+
+ /*
+ create a multi-root finder using a string.
+ The names are those defined in the GSL manuals
+ after having remived the GSL prefix (gsl_multiroot_fsolver).
+ Default algorithm is "hybrids" (without derivative).
+ */
+ GSLMultiRootFinder(const char * name = 0);
+
+ /// destructor
+ virtual ~GSLMultiRootFinder();
+
+ private:
+ // usually copying is non trivial, so we make this unaccessible
+ GSLMultiRootFinder(const GSLMultiRootFinder &);
+ GSLMultiRootFinder & operator = (const GSLMultiRootFinder &);
+
+ public:
+
+ /// set the type of algorithm
+ void SetType(EType type);
+
+ /// set the type of algorithm using derivatives
+ void SetType(EDerivType type);
+
+ /*
+ add the list of functions f1(x1,..xn),...fn(x1,...xn). The list must contain pointers of
+ ROOT::Math::IMultiGenFunctions. The method requires the
+ the begin and end of the list iterator.
+ The list can be any stl container or a simple array of ROOT::Math::IMultiGenFunctions* or
+ whatever implementing an iterator.
+ If using a derivative type algorithm the function pointers must implement the
+ ROOOT::Math::IMultiGradFunction interface
+ */
+ template<class FuncIterator>
+ bool SetFunctionList( FuncIterator begin, FuncIterator end) {
+ bool ret = true;
+ for (FuncIterator itr = begin; itr != end; ++itr) {
+ const ROOT::Math::IMultiGenFunction * f = *itr;
+ ret &= AddFunction( *f);
+ }
+ return ret;
+ }
+
+ /*
+ add (set) a single function fi(x1,...xn) which is part of the system of
+ specifying the begin and end of the iterator.
+ If using a derivative type algorithm the function must implement the
+ ROOOT::Math::IMultiGradFunction interface
+ Return the current number of function in the list and 0 if failed to add the function
+ */
+ int AddFunction( const ROOT::Math::IMultiGenFunction & func);
+
+ /// same method as before but using any function implementing
+ /// the operator(), so can be wrapped in a IMultiGenFunction interface
+ template <class Function>
+ int AddFunction( Function & f, int ndim) {
+ // no need to care about lifetime of wfunc. It will be cloned inside AddFunction
+ WrappedMultiFunction<Function &> wfunc(f, ndim);
+ return AddFunction(wfunc);
+ }
+
+ /**
+ return the number of sunctions set in the class.
+ The number must be equal to the dimension of the functions
+ */
+ unsigned int Dim() const { return fFunctions.size(); }
+
+ /// clear list of functions
+ void Clear();
+
+ /// return the root X values solving the system
+ const double * X() const;
+
+ /// return the function values f(X) solving the system
+ /// i.e. they must be close to zero at the solution
+ const double * FVal() const;
+
+ /// return the last step size
+ const double * Dx() const;
+
+
+ /**
+ Find the root starting from the point X;
+ Use the number of iteration and tolerance if given otherwise use
+ default parameter values which can be defined by
+ the static method SetDefault...
+ */
+ bool Solve(const double * x, int maxIter = 0, double absTol = 0, double relTol = 0);
+
+ /// Return number of iterations
+ int Iterations() const {
+ return fIter;
+ }
+
+ /// Return the status of last root finding
+ int Status() const { return fStatus; }
+
+ /// Return the algorithm name
+ const char * Name() const;
+
+ /*
+ set print level
+ level = 0 quiet (no messages print)
+ = 1 print only the result
+ = 3 max debug. Print result at each iteration
+ */
+ void SetPrintLevel(int level) { fPrintLevel = level; }
+
+ /// return the print level
+ int PrintLevel() const { return fPrintLevel; }
+
+
+ //-- static methods to set configurations
+
+ /// set tolerance (absolute and relative)
+ /// relative tolerance is only use to verify the convergence
+ /// do it is a minor parameter
+ static void SetDefaultTolerance(double abstol, double reltol = 0 );
+
+ /// set maximum number of iterations
+ static void SetDefaultMaxIterations(int maxiter);
+
+ /// print iteration state
+ void PrintState(std::ostream & os = std::cout);
+
+
+ protected:
+
+ // return type given a name
+ std::pair<bool,int> GetType(const char * name);
+ // clear list of functions
+ void ClearFunctions();
+
+
+ private:
+
+ int fIter; // current numer of iterations
+ int fStatus; // current status
+ int fPrintLevel; // print level
+
+ // int fMaxIter; // max number of iterations
+ // double fAbsTolerance; // absolute tolerance
+ // double fRelTolerance; // relative tolerance
+ int fType; // type of algorithm
+ bool fUseDerivAlgo; // algorithm using derivative
+
+ GSLMultiRootBaseSolver * fSolver;
+ std::vector<ROOT::Math::IMultiGenFunction *> fFunctions; //! transient Vector of the functions
+
+
+ };
+
+ // use typedef for most sensible name
+ typedef GSLMultiRootFinder MultiRootFinder;
+
+} // namespace Math
+} // namespace ROOT
+
+
+#endif /* ROOT_Math_GSLMultiRootFinder */
diff --git a/math/mathmore/inc/Math/LinkDef.h b/math/mathmore/inc/Math/LinkDef.h
index 8aadb97202100e70c1ae53484e0adb39930461b1..8abb735b88fa4909da53a269acc4fe15d049a355 100644
--- a/math/mathmore/inc/Math/LinkDef.h
+++ b/math/mathmore/inc/Math/LinkDef.h
@@ -103,6 +103,7 @@
#pragma link C++ class ROOT::Math::VegasParameters+;
#pragma link C++ class ROOT::Math::MiserParameters+;
-
+#pragma link C++ class ROOT::Math::GSLMultiRootFinder+;
+#pragma link C++ typedef ROOT::Math::MultiRootFinder;
#endif //__CINT__
diff --git a/math/mathmore/inc/Math/MultiRootFinder.h b/math/mathmore/inc/Math/MultiRootFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..ac444b2a6f39fe3d0955439e4c59313a1feb88d2
--- /dev/null
+++ b/math/mathmore/inc/Math/MultiRootFinder.h
@@ -0,0 +1,45 @@
+// @(#)root/mathmore:$Id$
+// Author: L. Moneta 03/2011
+
+ /**********************************************************************
+ * *
+ * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Header file for class MultiRootFinder
+//
+#ifndef ROOT_Math_MultiRootFinder
+#define ROOT_Math_MultiRootFinder
+
+
+
+#ifndef ROOT_Math_GSLMultiRootFinder
+#include "Math/GSLMultiRootFinder.h"
+#endif
+
+namespace ROOT {
+namespace Math {
+
+ typedef GSLMultiRootFinder MultiRootFinder;
+
+} // namespace Math
+} // namespace ROOT
+
+
+#endif /* ROOT_Math_MultiRootFinder */
diff --git a/math/mathmore/src/GSLMultiFit.h b/math/mathmore/src/GSLMultiFit.h
index a804fffdf6f60e40fd5a69a2a68a5d44304c7333..3713101770403f1d67737494e0e994ff18706651 100644
--- a/math/mathmore/src/GSLMultiFit.h
+++ b/math/mathmore/src/GSLMultiFit.h
@@ -54,14 +54,16 @@ public:
/**
Default constructor
- No need to specify the type sofar since only one solver exists so far
+ No need to specify the type so far since only one solver exists so far
*/
- GSLMultiFit () :
+ GSLMultiFit (const gsl_multifit_fdfsolver_type * type = 0) :
fSolver(0),
fVec(0),
fCov(0),
- fType(gsl_multifit_fdfsolver_lmsder)
- {}
+ fType(type)
+ {
+ if (fType == 0) fType = gsl_multifit_fdfsolver_lmsder; // default value
+ }
/**
Destructor (no operations)
diff --git a/math/mathmore/src/GSLMultiRootFinder.cxx b/math/mathmore/src/GSLMultiRootFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..efb84febebefce915d29f451caf0a4a290f91d9d
--- /dev/null
+++ b/math/mathmore/src/GSLMultiRootFinder.cxx
@@ -0,0 +1,346 @@
+// @(#)root/mathmore:$Id$
+// Authors: L. Moneta, A. Zsenei 08/2005
+
+ /**********************************************************************
+ * *
+ * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Implementation file for class GSLMultiRootFinder
+//
+// Created by: moneta at Sun Nov 14 11:27:11 2004
+//
+// Last update: Sun Nov 14 11:27:11 2004
+//
+
+#include "Math/IFunction.h"
+#include "Math/GSLMultiRootFinder.h"
+#include "GSLMultiRootSolver.h"
+#include "Math/Error.h"
+
+#include "gsl/gsl_multiroots.h"
+#include "gsl/gsl_errno.h"
+#include <cmath>
+#include <iomanip>
+
+#include <algorithm>
+#include <functional>
+#include <ctype.h> // need to use c version of tolower defined here
+
+
+namespace ROOT {
+namespace Math {
+
+ // default values
+
+ double gDefaultMaxIter = 100;
+ double gDefaultAbsTolerance = 1.E-6;
+ double gDefaultRelTolerance = 1.E-10;
+
+// impelmentation of static methods
+void GSLMultiRootFinder::SetDefaultTolerance(double abstol, double reltol ) {
+ // set default tolerance
+ gDefaultAbsTolerance = abstol;
+ if (reltol > 0) gDefaultRelTolerance = reltol;
+}
+void GSLMultiRootFinder::SetDefaultMaxIterations(int maxiter) {
+ // set default max iter
+ gDefaultMaxIter = maxiter;
+}
+
+GSLMultiRootFinder::GSLMultiRootFinder(EType type) :
+ fIter(0), fStatus(-1), fPrintLevel(0),
+ fType(type), fUseDerivAlgo(false),
+ fSolver(0)
+{
+ // constructor for non derivative type
+ fFunctions.reserve(2);
+}
+
+GSLMultiRootFinder::GSLMultiRootFinder(EDerivType type) :
+ fIter(0), fStatus(-1), fPrintLevel(0),
+ fType(type), fUseDerivAlgo(true),
+ fSolver(0)
+{
+ // constructor for non derivative type
+ fFunctions.reserve(2);
+}
+
+GSLMultiRootFinder::GSLMultiRootFinder(const char * name) :
+ fIter(0), fStatus(-1), fPrintLevel(0),
+ fType(0), fUseDerivAlgo(false),
+ fSolver(0)
+{
+ // constructor for a string
+ fFunctions.reserve(2);
+ std::pair<bool,int> type = GetType(name);
+ fUseDerivAlgo = type.first;
+ fType = type.second;
+}
+
+GSLMultiRootFinder::~GSLMultiRootFinder()
+{
+ // delete function wrapper
+ ClearFunctions();
+ if (fSolver) delete fSolver;
+}
+
+GSLMultiRootFinder::GSLMultiRootFinder(const GSLMultiRootFinder &)
+{
+}
+
+GSLMultiRootFinder & GSLMultiRootFinder::operator = (const GSLMultiRootFinder &rhs)
+{
+ // dummy operator=
+ if (this == &rhs) return *this; // time saving self-test
+
+ return *this;
+}
+
+int GSLMultiRootFinder::AddFunction(const ROOT::Math::IMultiGenFunction & func) {
+ // add a new function in the vector
+ ROOT::Math::IMultiGenFunction * f = func.Clone();
+ if (!f) return 0;
+ fFunctions.push_back(f);
+ return fFunctions.size();
+}
+
+void GSLMultiRootFinder::ClearFunctions() {
+ // clear the function list
+ for (unsigned int i = 0; i < fFunctions.size(); ++i) {
+ if (fFunctions[i] != 0 ) delete fFunctions[i];
+ fFunctions[i] = 0;
+ }
+ fFunctions.clear();
+}
+
+void GSLMultiRootFinder::Clear() {
+ // clear the function list and the solver
+ ClearFunctions();
+ if (fSolver) Clear();
+ fSolver = 0;
+}
+
+
+const double * GSLMultiRootFinder::X() const {
+ // return x
+ return (fSolver != 0) ? fSolver->X() : 0;
+}
+const double * GSLMultiRootFinder::Dx() const {
+ // return x
+ return (fSolver != 0) ? fSolver->Dx() : 0;
+}
+const double * GSLMultiRootFinder::FVal() const {
+ // return x
+ return (fSolver != 0) ? fSolver->FVal() : 0;
+}
+const char * GSLMultiRootFinder::Name() const {
+ // get GSL name
+ return (fSolver != 0) ? fSolver->Name().c_str() : "";
+}
+
+// bool GSLMultiRootFinder::AddFunction( const ROOT::Math::IMultiGenFunction & func) {
+// // clone and add function to the list
+// // If using a derivative algorithm the function is checked if it implements
+// // the gradient interface. If this is not the case the type is set to non-derivatibe algo
+// ROOT::Math::IGenMultiFunction * f = func.Clone();
+// if (f != 0) return false;
+// if (fUseDerivAlgo) {
+// bool gradFunc = (dynamic_cast<ROOT::Math::IMultiGradFunction *> (f) != 0 );
+// if (!gradFunc) {
+// MATH_ERROR_MSG("GSLMultiRootFinder::AddFunction","Function does not provide gradient interface");
+// MATH_WARN_MSG("GSLMultiRootFinder::AddFunction","clear the function list");
+// ClearFunctions();
+// return false;
+// }
+// }
+// fFunctions.push_back(f);
+// return true;
+// }
+
+ const gsl_multiroot_fsolver_type * GetGSLType(GSLMultiRootFinder::EType type) {
+ //helper functions to find GSL type
+ switch(type)
+ {
+ case ROOT::Math::GSLMultiRootFinder::kHybridS:
+ return gsl_multiroot_fsolver_hybrids;
+ case ROOT::Math::GSLMultiRootFinder::kHybrid:
+ return gsl_multiroot_fsolver_hybrid;
+ case ROOT::Math::GSLMultiRootFinder::kDNewton:
+ return gsl_multiroot_fsolver_dnewton;
+ case ROOT::Math::GSLMultiRootFinder::kBroyden:
+ return gsl_multiroot_fsolver_broyden;
+ default:
+ return gsl_multiroot_fsolver_hybrids;
+ }
+ return 0;
+}
+
+const gsl_multiroot_fdfsolver_type * GetGSLDerivType(GSLMultiRootFinder::EDerivType type) {
+//helper functions to find GSL deriv type
+ switch(type)
+ {
+ case ROOT::Math::GSLMultiRootFinder::kHybridSJ :
+ return gsl_multiroot_fdfsolver_hybridsj;
+ case ROOT::Math::GSLMultiRootFinder::kHybridJ :
+ return gsl_multiroot_fdfsolver_hybridj;
+ case ROOT::Math::GSLMultiRootFinder::kNewton :
+ return gsl_multiroot_fdfsolver_newton;
+ case ROOT::Math::GSLMultiRootFinder::kGNewton :
+ return gsl_multiroot_fdfsolver_gnewton;
+ default:
+ return gsl_multiroot_fdfsolver_hybridsj;
+ }
+ return 0; // cannot happen
+}
+
+std::pair<bool,int> GSLMultiRootFinder::GetType(const char * name) {
+ if (name == 0) return std::make_pair<bool,int>(false, -1);
+ std::string aname = name;
+ std::transform(aname.begin(), aname.end(), aname.begin(), (int(*)(int)) tolower );
+
+ if (aname.find("hybridsj") != std::string::npos) return std::make_pair(true, kHybridSJ);
+ if (aname.find("hybridj") != std::string::npos) return std::make_pair(true, kHybridJ);
+ if (aname.find("hybrids") != std::string::npos) return std::make_pair(false, kHybridS);
+ if (aname.find("hybrid") != std::string::npos) return std::make_pair(false, kHybrid);
+ if (aname.find("gnewton") != std::string::npos) return std::make_pair(true, kGNewton);
+ if (aname.find("dnewton") != std::string::npos) return std::make_pair(false, kDNewton);
+ if (aname.find("newton") != std::string::npos) return std::make_pair(true, kNewton);
+ if (aname.find("broyden") != std::string::npos) return std::make_pair(false, kBroyden);
+ MATH_INFO_MSG("GSLMultiRootFinder::GetType","Unknow algorithm - use default one");
+ return std::make_pair(false, -1);
+}
+
+bool GSLMultiRootFinder::Solve (const double * x, int maxIter, double absTol, double relTol)
+{
+ fIter = 0;
+ // create the solvers - delete previous existing solver
+ if (fSolver) delete fSolver;
+ fSolver = 0;
+
+ if (fFunctions.size() == 0) {
+ MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Function list is empty");
+ fStatus = -1;
+ return false;
+ }
+
+ if (fUseDerivAlgo) {
+ EDerivType type = (EDerivType) fType;
+ if (!fSolver) fSolver = new GSLMultiRootDerivSolver( GetGSLDerivType(type), Dim() );
+ }
+ else {
+ EType type = (EType) fType;
+ if (!fSolver) fSolver = new GSLMultiRootSolver( GetGSLType(type), Dim() );
+ }
+
+
+ // first set initial values and function
+ assert(fSolver != 0);
+ bool ret = fSolver->InitSolver( fFunctions, x);
+ if (!ret) {
+ MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Error initializing the solver");
+ fStatus = -2;
+ return false;
+ }
+
+ if (maxIter == 0) maxIter = gDefaultMaxIter;
+ if (absTol <= 0) absTol = gDefaultAbsTolerance;
+ if (relTol <= 0) relTol = gDefaultRelTolerance;
+
+ if (fPrintLevel >= 1)
+ std::cout << "GSLMultiRootFinder::Solve:" << Name() << " max iterations " << maxIter << " and tolerance " << absTol << std::endl;
+
+ // find the roots by iterating
+ fStatus = 0;
+ int status = 0;
+ int iter = 0;
+ do {
+ iter++;
+ status = fSolver->Iterate();
+
+ if (fPrintLevel >= 2) {
+ std::cout << "GSLMultiRootFinder::Solve - iteration # " << iter << " status = " << status << std::endl;
+ PrintState();
+ }
+ // act in case of error
+ if (status == GSL_EBADFUNC) {
+ MATH_ERROR_MSG("GSLMultiRootFinder::Solve","The iteration encountered a singolar point due to a bad function value");
+ fStatus = status;
+ break;
+ }
+ if (status == GSL_ENOPROG) {
+ MATH_ERROR_MSG("GSLMultiRootFinder::Solve","The iteration is not making any progress");
+ fStatus = status;
+ break;
+ }
+ if (status != GSL_SUCCESS) {
+ MATH_ERROR_MSG("GSLMultiRootFinder::Solve","Uknown iteration error - exit");
+ fStatus = status;
+ break;
+ }
+
+ // test also residual
+ status = fSolver->TestResidual(absTol);
+
+
+ // should test also the Delta ??
+ int status2 = fSolver->TestDelta(absTol, relTol);
+ if (status2 == GSL_SUCCESS) {
+ MATH_INFO_MSG("GSLMultiRootFinder::Solve","The iteration converged");
+ }
+ }
+ while (status == GSL_CONTINUE && iter < maxIter);
+ if (status == GSL_CONTINUE) {
+ MATH_INFO_MSGVAL("GSLMultiRootFinder::Solve","exceeded max iterations, reached tolerance is not sufficient",absTol);
+ }
+ if (status == GSL_SUCCESS) {
+ if (fPrintLevel>=1) { // print the result
+ MATH_INFO_MSG("GSLMultiRootFinder::Solve","The iteration converged");
+ std::cout << "GSL Algorithm used is : " << fSolver->Name() << std::endl;
+ std::cout << "Number of iterations = " << iter<< std::endl;
+
+ PrintState();
+ }
+ }
+ fIter = iter;
+ fStatus = status;
+ return (fStatus == GSL_SUCCESS);
+
+}
+
+void GSLMultiRootFinder::PrintState(std::ostream & os) {
+ // print current state
+ if (!fSolver) return;
+ int wi = int(std::log10(Dim() ) )+1;
+ const double * xtmp = fSolver->X();
+ const double * ftmp = fSolver->FVal();
+ os << "Root values = ";
+ for (unsigned int i = 0; i< Dim(); ++i)
+ os << "x[" << std::setw(wi) << i << "] = " << std::setw(12) << xtmp[i] << " ";
+ os << std::endl;
+ os << "Function values = ";
+ for (unsigned int i = 0; i< Dim(); ++i)
+ os << "f[" << std::setw(wi) << i << "] = " << std::setw(12) << ftmp[i] << " ";
+ os << std::endl;
+}
+
+
+
+} // namespace Math
+} // namespace ROOT
diff --git a/math/mathmore/src/GSLMultiRootFunctionAdapter.h b/math/mathmore/src/GSLMultiRootFunctionAdapter.h
new file mode 100644
index 0000000000000000000000000000000000000000..f796182ce09e68ba23c2af6963146812148f55a6
--- /dev/null
+++ b/math/mathmore/src/GSLMultiRootFunctionAdapter.h
@@ -0,0 +1,130 @@
+// @(#)root/mathmore:$Id$
+// Authors: L. Moneta, Mar 2011
+
+ /**********************************************************************
+ * *
+ * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Header file for class GSLMultiMinFunctionAdapter
+//
+// Generic adapter for gsl_multiroot_function signature type
+// usable for any array of function pointers
+// implementing operator()(const double *x) and (if needed)
+// Gradient(const double *x, double * g)
+//
+// The class is very similar to GSLMultiFitFunctionAdapter,
+// but in that case the array is for function references (or value)
+//
+#ifndef ROOT_Math_GSLMultiRootFunctionAdapter
+#define ROOT_Math_GSLMultiRootFunctionAdapter
+
+#include "gsl/gsl_vector.h"
+#include "gsl/gsl_matrix.h"
+
+#include <cassert>
+
+namespace ROOT {
+namespace Math {
+
+
+
+ /**
+ Class for adapting a C++ functor class to C function pointers used by GSL MultiRoot
+ Algorithm
+ The templated C++ function class must implement:
+
+ <em> double operator( const double * x)</em>
+ and if the derivatives are required:
+ <em> void Gradient( const double * x, double * g)</em>
+ and
+ <em> void FdF( const double * x, double &f, double * g)</em>
+
+
+ @ingroup MultiRoot
+ */
+
+
+ // FuncVector must contain a vector of pointers to functions
+ // this same as MultiFit but here need to use pointers where there we used class elements
+
+template<class FuncVector>
+class GSLMultiRootFunctionAdapter {
+
+
+
+public:
+
+ static int F( const gsl_vector * x, void * p, gsl_vector * f ) {
+ // p is a pointer to an iterator of functions
+ unsigned int n = f->size;
+ // need to copy iterator otherwise next time the function is called it wont work
+ FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) );
+ if (n == 0) return -1;
+ for (unsigned int i = 0; i < n ; ++i) {
+ gsl_vector_set(f, i, (*funcVec[i])(x->data) );
+ }
+ return 0;
+ }
+
+
+ static int Df( const gsl_vector * x, void * p, gsl_matrix * h) {
+
+ // p is a pointer to an iterator of functions
+ unsigned int n = h->size1;
+ unsigned int npar = h->size2;
+ if (n == 0) return -1;
+ if (npar == 0) return -2;
+ FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) );
+ for (unsigned int i = 0; i < n ; ++i) {
+ double * g = (h->data)+i*npar; //pointer to start of i-th row
+ assert ( npar == (funcVec[i])->NDim() );
+ (funcVec[i])->Gradient(x->data, g);
+ }
+ return 0;
+ }
+
+ /// evaluate derivative and function at the same time
+ static int FDf( const gsl_vector * x, void * p, gsl_vector * f, gsl_matrix * h) {
+ // should be implemented in the function
+ // p is a pointer to an iterator of functions
+ unsigned int n = h->size1;
+ unsigned int npar = h->size2;
+ if (n == 0) return -1;
+ if (npar == 0) return -2;
+ FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) );
+ assert ( f->size == n);
+ for (unsigned int i = 0; i < n ; ++i) {
+ assert ( npar == (funcVec[i])->NDim() );
+ double fval = 0;
+ double * g = (h->data)+i*npar; //pointer to start of i-th row
+ (funcVec[i])->FdF(x->data, fval, g);
+ gsl_vector_set(f, i, fval );
+ }
+ return 0;
+ }
+
+};
+
+
+} // namespace Math
+} // namespace ROOT
+
+
+#endif /* ROOT_Math_GSLMultiRootFunctionAdapter */
diff --git a/math/mathmore/src/GSLMultiRootFunctionWrapper.h b/math/mathmore/src/GSLMultiRootFunctionWrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..9b531c37479389447558229ddea9e732d4366379
--- /dev/null
+++ b/math/mathmore/src/GSLMultiRootFunctionWrapper.h
@@ -0,0 +1,138 @@
+// @(#)root/mathmore:$Id$
+// Authors: L. Moneta Dec 2006
+
+ /**********************************************************************
+ * *
+ * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Header file for class GSLMultiRootFunctionWrapper
+//
+// Created by: moneta at Sat Nov 13 14:54:41 2004
+//
+// Last update: Sat Nov 13 14:54:41 2004
+//
+#ifndef ROOT_Math_GSLMultiRootFunctionWrapper
+#define ROOT_Math_GSLMultiRootFunctionWrapper
+
+#include "gsl/gsl_multiroots.h"
+
+#include "GSLMultiRootFunctionAdapter.h"
+
+
+#include <cassert>
+
+namespace ROOT {
+namespace Math {
+
+
+
+// can re-use same type for multi-fit
+
+ typedef double ( * GSLMultiRootFPointer ) ( const gsl_vector *, void *, gsl_vector *);
+ typedef void ( * GSLMultiRootDfPointer ) ( const gsl_vector *, void *, gsl_matrix *);
+ typedef void ( * GSLMultiRootFdfPointer ) ( const gsl_vector *, void *, gsl_vector *, gsl_matrix *);
+
+
+/**
+ wrapper to a multi-dim function without derivatives for multi roots
+ algorithm
+
+ @ingroup MultiRoots
+*/
+class GSLMultiRootFunctionWrapper {
+
+public:
+
+ GSLMultiRootFunctionWrapper()
+ {
+ fFunc.f = 0;
+ fFunc.n = 0;
+ fFunc.params = 0;
+ }
+
+
+ /// Fill gsl function structure from a C++ function iterator and size and number of residuals
+ template<class FuncVector>
+ void SetFunctions(const FuncVector & f, unsigned int n ) {
+ const void * p = &f;
+ assert (p != 0);
+ fFunc.f = &GSLMultiRootFunctionAdapter<FuncVector >::F;
+ fFunc.n = n;
+ fFunc.params = const_cast<void *>(p);
+ }
+
+ gsl_multiroot_function * GetFunctions() { return &fFunc; }
+
+
+ private:
+
+ gsl_multiroot_function fFunc;
+
+};
+
+
+/**
+ wrapper to a multi-dim function with derivatives for multi roots
+ algorithm
+
+ @ingroup MultiRoot
+*/
+
+class GSLMultiRootDerivFunctionWrapper {
+
+public:
+
+ GSLMultiRootDerivFunctionWrapper()
+ {
+ fFunc.f = 0;
+ fFunc.df = 0;
+ fFunc.fdf = 0;
+ fFunc.n = 0;
+ fFunc.params = 0;
+ }
+
+
+ /// Fill gsl function structure from a C++ function iterator and size and number of residuals
+ template<class FuncVector>
+ void SetFunctions(const FuncVector & f, unsigned int n ) {
+ const void * p = &f;
+ assert (p != 0);
+ fFunc.f = &GSLMultiRootFunctionAdapter<FuncVector >::F;
+ fFunc.df = &GSLMultiRootFunctionAdapter<FuncVector >::Df;
+ fFunc.fdf = &GSLMultiRootFunctionAdapter<FuncVector >::FDf;
+ fFunc.n = n;
+ fFunc.params = const_cast<void *>(p);
+ }
+
+ gsl_multiroot_function_fdf * GetFunctions() { return &fFunc; }
+
+
+ private:
+
+ gsl_multiroot_function_fdf fFunc;
+
+};
+
+
+
+} // namespace Math
+} // namespace ROOT
+
+#endif /* ROOT_Math_GSLMultiMinFunctionWrapper */
diff --git a/math/mathmore/src/GSLMultiRootSolver.h b/math/mathmore/src/GSLMultiRootSolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..4b557af670fdf0f6482c2f17c3a5008dfb3b2225
--- /dev/null
+++ b/math/mathmore/src/GSLMultiRootSolver.h
@@ -0,0 +1,391 @@
+// @(#)root/mathmore:$Id$
+// Author: L. Moneta Wed Dec 20 17:26:06 2006
+
+/**********************************************************************
+ * *
+ * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT *
+ * *
+ * This library is free software; you can redistribute it and/or *
+ * modify it under the terms of the GNU General Public License *
+ * as published by the Free Software Foundation; either version 2 *
+ * of the License, or (at your option) any later version. *
+ * *
+ * This library is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
+ * General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this library (see file COPYING); if not, write *
+ * to the Free Software Foundation, Inc., 59 Temple Place, Suite *
+ * 330, Boston, MA 02111-1307 USA, or contact the author. *
+ * *
+ **********************************************************************/
+
+// Header file for the class GSLMultiRootBaseSolver,
+// GSLMultiRootSolver and GSLMultiRootDerivSolver
+
+#ifndef ROOT_Math_GSLMultiRootSolver
+#define ROOT_Math_GSLMultiRootSolver
+
+#include "gsl/gsl_vector.h"
+#include "gsl/gsl_matrix.h"
+#include "gsl/gsl_multiroots.h"
+#include "gsl/gsl_blas.h"
+#include "GSLMultiRootFunctionWrapper.h"
+
+#include "Math/IFunction.h"
+#include "Math/Error.h"
+
+#include <vector>
+#include <string>
+#include <cassert>
+
+
+namespace ROOT {
+
+ namespace Math {
+
+
+/**
+ GSLMultiRootBaseSolver, internal class for implementing GSL multi-root finders
+ This is the base class for GSLMultiRootSolver (solver not using derivatives) and
+ GSLMUltiRootDerivSolver (solver using derivatives)
+
+ @ingroup MultiRoot
+*/
+class GSLMultiRootBaseSolver {
+
+public:
+
+ /**
+ virtual Destructor
+ */
+ virtual ~GSLMultiRootBaseSolver () {}
+
+
+public:
+
+
+ /// init the solver with function list and initial values
+ bool InitSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
+ // create a vector of the fit contributions
+ // create function wrapper from an iterator of functions
+ unsigned int n = funcVec.size();
+ if (n == 0) return false;
+
+ unsigned int ndim = funcVec[0]->NDim(); // should also be = n
+
+ if (ndim != n) {
+ MATH_ERROR_MSGVAL("GSLMultiRootSolver::InitSolver","Wrong function dimension",ndim);
+ MATH_ERROR_MSGVAL("GSLMultiRootSolver::InitSolver","Number of functions",n);
+ return false;
+ }
+
+
+ // set function list and initial values in solver
+ int iret = SetSolver(funcVec,x);
+ return (iret == 0);
+ }
+
+ /// return name
+ virtual std::string Name() const = 0;
+
+ /// perform an iteration
+ virtual int Iterate() = 0;
+
+ /// solution values at the current iteration
+ const double * X() const {
+ gsl_vector * x = GetRoot();
+ return x->data;
+ }
+
+ /// return function values
+ const double * FVal() const {
+ gsl_vector * f = GetF();
+ return f->data;
+ }
+
+ /// return function steps
+ const double * Dx() const {
+ gsl_vector * dx = GetDx();
+ return dx->data;
+ }
+
+ /// test using abs and relative tolerance
+ /// |dx| < absTol + relTol*|x| for every component
+ int TestDelta(double absTol, double relTol) const {
+ gsl_vector * x = GetRoot();
+ gsl_vector * dx = GetDx();
+ if (x == 0 || dx == 0) return -1;
+ return gsl_multiroot_test_delta(dx, x, absTol, relTol);
+ }
+
+ /// test using abs tolerance
+ /// Sum |f|_i < absTol
+ int TestResidual(double absTol) const {
+ gsl_vector * f = GetF();
+ if (f == 0) return -1;
+ return gsl_multiroot_test_residual(f, absTol);
+ }
+
+
+private:
+
+ // accessor to be implemented by the derived classes
+
+ virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) = 0;
+
+ virtual gsl_vector * GetRoot() const = 0;
+
+ virtual gsl_vector * GetF() const = 0;
+
+ virtual gsl_vector * GetDx() const = 0;
+
+
+};
+
+
+/**
+ GSLMultiRootSolver, internal class for implementing GSL multi-root finders
+ not using derivatives
+
+ @ingroup MultiRoot
+*/
+class GSLMultiRootSolver : public GSLMultiRootBaseSolver {
+
+public:
+
+ /**
+ Constructor from type and simension of system (number of functions)
+ */
+ GSLMultiRootSolver (const gsl_multiroot_fsolver_type * type, int n ) :
+ fSolver(0),
+ fVec(0)
+ {
+ CreateSolver(type, n);
+ }
+
+ /**
+ Destructor (no operations)
+ */
+ virtual ~GSLMultiRootSolver () {
+ if (fSolver) gsl_multiroot_fsolver_free(fSolver);
+ if (fVec != 0) gsl_vector_free(fVec);
+ }
+
+private:
+ // usually copying is non trivial, so we make this unaccessible
+
+ /**
+ Copy constructor
+ */
+ GSLMultiRootSolver(const GSLMultiRootSolver &) : GSLMultiRootBaseSolver() {}
+
+ /**
+ Assignment operator
+ */
+ GSLMultiRootSolver & operator = (const GSLMultiRootSolver & rhs) {
+ if (this == &rhs) return *this; // time saving self-test
+ return *this;
+ }
+
+
+public:
+
+
+ void CreateSolver(const gsl_multiroot_fsolver_type * type, unsigned int n) {
+
+ /// create the solver from the type and size of number of fitting points and number of parameters
+ if (fSolver) gsl_multiroot_fsolver_free(fSolver);
+ fSolver = gsl_multiroot_fsolver_alloc(type, n);
+
+ }
+
+
+ /// set the solver parameters
+ virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
+ // create a vector of the fit contributions
+ // create function wrapper from an iterator of functions
+ assert(fSolver !=0);
+ unsigned int n = funcVec.size();
+
+ fFunctions.SetFunctions(funcVec, funcVec.size() );
+ // set initial values and create cached vector
+ if (fVec != 0) gsl_vector_free(fVec);
+ fVec = gsl_vector_alloc( n);
+ std::copy(x,x+n, fVec->data);
+ // solver should have been already created at this point
+ assert(fSolver != 0);
+ return gsl_multiroot_fsolver_set(fSolver, fFunctions.GetFunctions(), fVec);
+ }
+
+ virtual std::string Name() const {
+ if (fSolver == 0 ) return "undefined";
+ return std::string(gsl_multiroot_fsolver_name(fSolver) );
+ }
+
+ virtual int Iterate() {
+ if (fSolver == 0) return -1;
+ return gsl_multiroot_fsolver_iterate(fSolver);
+ }
+
+ /// solution values at the current iteration
+ virtual gsl_vector * GetRoot() const {
+ if (fSolver == 0) return 0;
+ return gsl_multiroot_fsolver_root(fSolver);
+ }
+
+ /// return function values
+ virtual gsl_vector * GetF() const {
+ if (fSolver == 0) return 0;
+ return gsl_multiroot_fsolver_f(fSolver);
+ }
+
+ /// return function steps
+ virtual gsl_vector * GetDx() const {
+ if (fSolver == 0) return 0;
+ return gsl_multiroot_fsolver_dx(fSolver);
+ }
+
+
+private:
+
+ GSLMultiRootFunctionWrapper fFunctions;
+ gsl_multiroot_fsolver * fSolver;
+ // cached vector to avoid re-allocating every time a new one
+ mutable gsl_vector * fVec;
+ const gsl_multiroot_fsolver_type * fType;
+
+};
+
+/**
+ GSLMultiRootDerivSolver, internal class for implementing GSL multi-root finders
+ using derivatives
+
+ @ingroup MultiRoot
+*/
+class GSLMultiRootDerivSolver : public GSLMultiRootBaseSolver {
+
+public:
+
+ /**
+ Constructor
+ */
+ GSLMultiRootDerivSolver (const gsl_multiroot_fdfsolver_type * type, int n ) :
+ fDerivSolver(0),
+ fVec(0)
+ {
+ CreateSolver(type, n);
+ }
+
+ /**
+ Destructor (no operations)
+ */
+ virtual ~GSLMultiRootDerivSolver () {
+ if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
+ if (fVec != 0) gsl_vector_free(fVec);
+ }
+
+private:
+ // usually copying is non trivial, so we make this unaccessible
+
+ /**
+ Copy constructor
+ */
+ GSLMultiRootDerivSolver(const GSLMultiRootDerivSolver &) : GSLMultiRootBaseSolver() {}
+
+ /**
+ Assignment operator
+ */
+ GSLMultiRootDerivSolver & operator = (const GSLMultiRootDerivSolver & rhs) {
+ if (this == &rhs) return *this; // time saving self-test
+ return *this;
+ }
+
+
+public:
+
+
+ /// create the solver from the type and size of number of fitting points and number of parameters
+ void CreateSolver(const gsl_multiroot_fdfsolver_type * type, unsigned int n) {
+
+ /// create the solver from the type and size of number of fitting points and number of parameters
+ if (fDerivSolver) gsl_multiroot_fdfsolver_free(fDerivSolver);
+ fDerivSolver = gsl_multiroot_fdfsolver_alloc(type, n);
+ }
+
+
+
+ /// set the solver parameters for the case of derivative
+ virtual int SetSolver(const std::vector<ROOT::Math::IMultiGenFunction*> & funcVec, const double * x) {
+ // create a vector of the fit contributions
+ // need to create a vecctor of gradient functions, convert and store in the class
+ // the new vector pointer
+ assert(fDerivSolver !=0);
+ unsigned int n = funcVec.size();
+ fGradFuncVec.reserve( n );
+ for (unsigned int i = 0; i < n; ++i) {
+ ROOT::Math::IMultiGradFunction * func = dynamic_cast<ROOT::Math::IMultiGradFunction *>(funcVec[i] );
+ if (func == 0) {
+ MATH_ERROR_MSG("GSLMultiRootSolver::SetSolver","Function does not provide gradient interface");
+ return -1;
+ }
+ fGradFuncVec.push_back( func);
+ }
+
+ fDerivFunctions.SetFunctions(fGradFuncVec, funcVec.size() );
+ // set initial values
+ if (fVec != 0) gsl_vector_free(fVec);
+ fVec = gsl_vector_alloc( n);
+ std::copy(x,x+n, fVec->data);
+
+ return gsl_multiroot_fdfsolver_set(fDerivSolver, fDerivFunctions.GetFunctions(), fVec);
+ }
+
+ virtual std::string Name() const {
+ if (fDerivSolver == 0 ) return "undefined";
+ return std::string(gsl_multiroot_fdfsolver_name(fDerivSolver) );
+ }
+
+ virtual int Iterate() {
+ if (fDerivSolver == 0) return -1;
+ return gsl_multiroot_fdfsolver_iterate(fDerivSolver);
+ }
+
+ /// solution values at the current iteration
+ virtual gsl_vector * GetRoot() const {
+ if (fDerivSolver == 0) return 0;
+ return gsl_multiroot_fdfsolver_root(fDerivSolver);
+ }
+
+ /// return function values
+ virtual gsl_vector * GetF() const {
+ if (fDerivSolver == 0) return 0;
+ return gsl_multiroot_fdfsolver_f(fDerivSolver);
+ }
+
+ /// return function steps
+ virtual gsl_vector * GetDx() const {
+ if (fDerivSolver == 0) return 0;
+ return gsl_multiroot_fdfsolver_dx(fDerivSolver);
+ }
+
+
+
+private:
+
+ GSLMultiRootDerivFunctionWrapper fDerivFunctions;
+ gsl_multiroot_fdfsolver * fDerivSolver;
+ // cached vector to avoid re-allocating every time a new one
+ mutable gsl_vector * fVec;
+ std::vector<ROOT::Math::IMultiGradFunction*> fGradFuncVec;
+
+};
+
+ } // end namespace Math
+
+} // end namespace ROOT
+
+
+#endif /* ROOT_Math_GSLMultiRootSolver */
diff --git a/math/mathmore/src/GSLNLSMinimizer.cxx b/math/mathmore/src/GSLNLSMinimizer.cxx
index daa40339d9e7621ae616115c1454b406e692770d..c72dc8ba9cf66c2025f78012762d78aa99e258f6 100644
--- a/math/mathmore/src/GSLNLSMinimizer.cxx
+++ b/math/mathmore/src/GSLNLSMinimizer.cxx
@@ -116,7 +116,7 @@ private:
// GSLNLSMinimizer implementation
-GSLNLSMinimizer::GSLNLSMinimizer( int /* ROOT::Math::EGSLNLSMinimizerType type */ ) :
+GSLNLSMinimizer::GSLNLSMinimizer( int type ) :
fDim(0),
fNFree(0),
fSize(0),
@@ -124,21 +124,31 @@ GSLNLSMinimizer::GSLNLSMinimizer( int /* ROOT::Math::EGSLNLSMinimizerType type *
fMinVal(0)
{
// Constructor implementation : create GSLMultiFit wrapper object
- fGSLMultiFit = new GSLMultiFit( /*type */ );
+ const gsl_multifit_fdfsolver_type * gsl_type = 0; // use default type defined in GSLMultiFit
+ if (type == 1) gsl_type = gsl_multifit_fdfsolver_lmsder; // scaled lmder version
+ if (type == 2) gsl_type = gsl_multifit_fdfsolver_lmder; // unscaled version
+
+ fGSLMultiFit = new GSLMultiFit( gsl_type );
fValues.reserve(10);
fNames.reserve(10);
fSteps.reserve(10);
fEdm = -1;
- fLSTolerance = 0.0001;
- SetMaxIterations(100);
- SetPrintLevel(1);
+
+ // defult tolerance and max iterations
+ int niter = ROOT::Math::MinimizerOptions::DefaultMaxIterations();
+ if (niter <= 0) niter = 100;
+ SetMaxIterations(niter);
+
+ fLSTolerance = ROOT::Math::MinimizerOptions::DefaultTolerance();
+ if (fLSTolerance <=0) fLSTolerance = 0.0001; // default internal value
+
+ SetPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel());
}
GSLNLSMinimizer::~GSLNLSMinimizer () {
assert(fGSLMultiFit != 0);
delete fGSLMultiFit;
-// if (fObjFunc) delete fObjFunc;
}
bool GSLNLSMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) {
@@ -287,7 +297,7 @@ bool GSLNLSMinimizer::Minimize() {
startValues.resize( fNFree );
}
- if (debugLevel >=1 ) std::cout <<"Minimize using GSLNLSMinimizer " << fGSLMultiFit->Name() << std::endl;
+ if (debugLevel >=1 ) std::cout <<"Minimize using GSLNLSMinimizer " << std::endl;
// // use a global step size = min (step vectors)
// double stepSize = 1;
@@ -301,7 +311,7 @@ bool GSLNLSMinimizer::Minimize() {
return false;
}
- if (debugLevel >=1 ) std::cout <<"GSLNLSMinimizer: Start iterating......... " << std::endl;
+ if (debugLevel >=1 ) std::cout <<"GSLNLSMinimizer: " << fGSLMultiFit->Name() << " - start iterating......... " << std::endl;
// start iteration
unsigned int iter = 0;
diff --git a/math/mathmore/test/Makefile b/math/mathmore/test/Makefile
index 670ff7d55ba57ae5ea9adf418a6edabecf803d85..bd8b5c5f398970ae4f8f1fcd27810c1f04104263 100644
--- a/math/mathmore/test/Makefile
+++ b/math/mathmore/test/Makefile
@@ -67,6 +67,10 @@ ROOTFINDEROBJ = testRootFinder.$(ObjSuf)
ROOTFINDERSRC = testRootFinder.$(SrcSuf)
ROOTFINDER = testRootFinder$(ExeSuf)
+MULTIROOTFINDEROBJ = testMultiRootFinder.$(ObjSuf)
+MULTIROOTFINDERSRC = testMultiRootFinder.$(SrcSuf)
+MULTIROOTFINDER = testMultiRootFinder$(ExeSuf)
+
MINIMIZATION1DOBJ = testMinimization1D.$(ObjSuf)
MINIMIZATION1DSRC = testMinimization1D.$(SrcSuf)
MINIMIZATION1D = testMinimization1D$(ExeSuf)
@@ -106,11 +110,11 @@ SIMANTSP = simanTSP$(ExeSuf)
-OBJS = $(CHEBYSHEVOBJ) $(DERIVATIONOBJ) $(INTEGRATIONOBJ) $(INTEGRATIONMDOBJ) $(ROOTFINDEROBJ) \
+OBJS = $(CHEBYSHEVOBJ) $(DERIVATIONOBJ) $(INTEGRATIONOBJ) $(INTEGRATIONMDOBJ) $(ROOTFINDEROBJ) $(MULTIROOTFINDEROBJ)\
$(MINIMIZATION1DOBJ) $(INTERPOLATIONOBJ) $(RANDOMOBJ) $(RANDOMDISTOBJ) \
$(SPECFUNCOBJ) $(STATFUNCOBJ) $(TESTFUNCTOROBJ) $(TESTPERMUTEOBJ) $(TESTVAVILOVOBJ) $(SIMANTSPOBJ)
-PROGRAMS = $(CHEBYSHEV) $(DERIVATION) $(INTEGRATION) $(INTEGRATIONMD) $(ROOTFINDER) \
+PROGRAMS = $(CHEBYSHEV) $(DERIVATION) $(INTEGRATION) $(INTEGRATIONMD) $(ROOTFINDER) $(MULTIROOTFINDER) \
$(MINIMIZATION1D) $(INTERPOLATION) $(RANDOM) $(RANDOMDIST) \
$(SPECFUNC) $(STATFUNC) $(TESTFUNCTOR) $(TESTPERMUTE) $(TESTVAVILOV) $(SIMANTSP)
@@ -140,6 +144,10 @@ $(ROOTFINDER): $(ROOTFINDEROBJ)
$(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@
@echo "$@ done"
+$(MULTIROOTFINDER): $(MULTIROOTFINDEROBJ)
+ $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@
+ @echo "$@ done"
+
$(MINIMIZATION1D): $(MINIMIZATION1DOBJ)
$(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@
@echo "$@ done"
diff --git a/math/mathmore/test/testMultiRootFinder.cxx b/math/mathmore/test/testMultiRootFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ab62dd98abea4919549bce7c5e8f52b680563f1d
--- /dev/null
+++ b/math/mathmore/test/testMultiRootFinder.cxx
@@ -0,0 +1,118 @@
+#include "Math/Functor.h"
+#include "Math/MultiRootFinder.h"
+
+#ifdef HAVE_ROOTLIBS
+#include "TStopwatch.h"
+#else
+struct TStopwatch {
+ void Start(){}
+ void Stop(){}
+ void Reset(){}
+ double RealTime() { return 0; }
+ double CpuTime() { return 0; }
+};
+#endif
+
+#include <iostream>
+
+// solve Roots of rosenbrock function
+// f1(x,y) = a(1-x)
+// f2(x,y) = b(y-x^2)
+// with 1 = 1, b=10
+
+
+// define system of functions to find the roots
+struct FuncSystem {
+
+ double F1(const double *xx) {
+ double x = xx[0];
+ return a * (1. - x );
+ }
+ double F2(const double *xx) {
+ double x = xx[0]; double y = xx[1];
+ return b * (y - x*x );
+ }
+
+ // derivative
+ double DerivF1(const double *, int icoord) {
+ if (icoord == 0) return -a;
+ else return 0;
+ }
+ double DerivF2(const double *xx, int icoord) {
+ double x = xx[0];
+ if (icoord == 0) return -2 * b * x;
+ else return b;
+ }
+
+ double a;
+ double b;
+};
+
+int printlevel = 0;
+
+using namespace ROOT::Math;
+
+int testMultiRootFinder() {
+
+ int status = 0;
+
+ // methods using derivatives
+ FuncSystem f;
+ f.a = 1;
+ f.b = 10;
+
+
+ MultiRootFinder rf(MultiRootFinder::kHybridSJ);
+ GradFunctor g1(&f, &FuncSystem::F1, &FuncSystem::DerivF1,2);
+ GradFunctor g2(&f, &FuncSystem::F2, &FuncSystem::DerivF2,2);
+ rf.AddFunction(g1);
+ rf.AddFunction(g2);
+ rf.SetPrintLevel(printlevel);
+
+ double x0[] = {-1.,-1.};
+
+ bool ret = rf.Solve(x0);
+
+ if (!ret) {
+ std::cerr << "testMultiRootFinder - Error running derivative algorithm " << std::endl;
+ if (printlevel == 0) rf.PrintState(std::cout);
+ status += rf.Status();
+ }
+
+ MultiRootFinder rf2(MultiRootFinder::kHybridS);
+ Functor f1(&f, &FuncSystem::F1, 2);
+ Functor f2(&f, &FuncSystem::F2, 2);
+ std::vector<ROOT::Math::IMultiGenFunction*> funlist;
+ funlist.push_back(&f1);
+ funlist.push_back(&f2);
+ rf2.SetFunctionList(funlist.begin(), funlist.end() );
+
+ rf2.SetPrintLevel(printlevel);
+
+ bool ret2 = rf2.Solve(x0);
+ if (!ret2) {
+ std::cerr << "testMultiRootFinder - Error running non-derivative algorithm " << std::endl;
+ if (printlevel == 0) rf2.PrintState(std::cout);
+ status += 10*rf2.Status();
+ }
+
+ return status;
+
+}
+
+int main (int argc, char **argv) {
+ int status = 0;
+
+ if (argc > 1 )
+ printlevel = atoi(argv[1]);
+
+ status += testMultiRootFinder();
+ if (status == 0) {
+ std::cout << "testMultiRootFinder --- \t" << "OK" << std::endl;
+ }
+ else {
+ std::cerr << "testMultiRootFinder --- \t" << "FAILED ! " << "\t with status = " << status << std::endl;
+ }
+
+ return status;
+}