Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:10:16

0001 // @(#)root/mathmore:$Id$
0002 // Authors: L. Moneta, A. Zsenei   08/2005
0003 
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2004 ROOT Foundation,  CERN/PH-SFT                   *
0007   *                                                                    *
0008   * This library is free software; you can redistribute it and/or      *
0009   * modify it under the terms of the GNU General Public License        *
0010   * as published by the Free Software Foundation; either version 2     *
0011   * of the License, or (at your option) any later version.             *
0012   *                                                                    *
0013   * This library is distributed in the hope that it will be useful,    *
0014   * but WITHOUT ANY WARRANTY; without even the implied warranty of     *
0015   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU   *
0016   * General Public License for more details.                           *
0017   *                                                                    *
0018   * You should have received a copy of the GNU General Public License  *
0019   * along with this library (see file COPYING); if not, write          *
0020   * to the Free Software Foundation, Inc., 59 Temple Place, Suite      *
0021   * 330, Boston, MA 02111-1307 USA, or contact the author.             *
0022   *                                                                    *
0023   **********************************************************************/
0024 
0025 // Header file for class GSLIntegrator
0026 //
0027 // Created by: Lorenzo Moneta  at Thu Nov 11 14:22:32 2004
0028 //
0029 // Last update: Thu Nov 11 14:22:32 2004
0030 //
0031 #ifndef ROOT_Math_GSLIntegrator
0032 #define ROOT_Math_GSLIntegrator
0033 
0034 
0035 #include "Math/VirtualIntegrator.h"
0036 
0037 #include "Math/IntegrationTypes.h"
0038 
0039 #include "Math/IFunctionfwd.h"
0040 
0041 
0042 
0043 
0044 #include "Math/GSLFunctionAdapter.h"
0045 
0046 #include <vector>
0047 
0048 
0049 
0050 namespace ROOT {
0051 namespace Math {
0052 
0053 
0054 
0055    class GSLIntegrationWorkspace;
0056    class GSLFunctionWrapper;
0057 
0058    //_________________________________________________________________________
0059    /**
0060 
0061    Class for performing numerical integration of a function in one dimension.
0062    It uses the numerical integration algorithms of GSL, which reimplements the
0063    algorithms used in the QUADPACK, a numerical integration package written in Fortran.
0064 
0065    Various types of adaptive and non-adaptive integration are supported. These include
0066    integration over infinite and semi-infinite ranges and singular integrals.
0067 
0068    The integration type is selected using the Integration::type enumeration
0069    in the class constructor.
0070    The default type is adaptive integration with singularity
0071    (ADAPTIVESINGULAR or QAGS in the QUADPACK convention) applying a Gauss-Kronrod 21-point integration rule.
0072    In the case of ADAPTIVE type, the integration rule can also be specified via the
0073    Integration::GKRule. The default rule is 31 points.
0074 
0075    In the case of integration over infinite and semi-infinite ranges, the type used is always
0076    ADAPTIVESINGULAR applying a transformation from the original interval into (0,1).
0077 
0078    The ADAPTIVESINGULAR type is the most sophicticated type. When performances are
0079    important, it is then recommended to use the NONADAPTIVE type in case of smooth functions or
0080    ADAPTIVE with a lower Gauss-Kronrod rule.
0081 
0082    For detailed description on GSL integration algorithms see the
0083    <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html">GSL Manual</A>.
0084 
0085 
0086    @ingroup Integration
0087    */
0088 
0089 
0090    class GSLIntegrator : public VirtualIntegratorOneDim  {
0091 
0092    public:
0093 
0094 
0095 
0096       // constructors
0097 
0098 
0099       /** Default constructor of GSL Integrator for Adaptive Singular integration
0100 
0101       @param absTol desired absolute Error
0102       @param relTol desired relative Error
0103       @param size maximum number of sub-intervals
0104       */
0105 
0106       GSLIntegrator(double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
0107 
0108 
0109 
0110 
0111       /** constructor of GSL Integrator. In the case of Adaptive integration the Gauss-Krond rule of 31 points is used
0112 
0113          @param type type of integration. The possible types are defined in the Integration::Type enumeration
0114          @param absTol desired absolute Error
0115          @param relTol desired relative Error
0116          @param size maximum number of sub-intervals
0117          */
0118 
0119 
0120       GSLIntegrator(const Integration::Type type, double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
0121 
0122 
0123       /**
0124          generic constructor for GSL Integrator
0125 
0126        @param type type of integration. The possible types are defined in the Integration::Type enumeration
0127        @param rule Gauss-Kronrod rule. It is used only for ADAPTIVE::Integration types. The possible rules are defined in the Integration::GKRule enumeration
0128        @param absTol desired absolute Error
0129        @param relTol desired relative Error
0130        @param size maximum number of sub-intervals
0131 
0132        */
0133 
0134       GSLIntegrator(const Integration::Type type, const Integration::GKRule rule, double absTol = 1.E-9, double relTol = 1E-6, size_t size = 1000);
0135 
0136 
0137       /** constructor of GSL Integrator. In the case of Adaptive integration the Gauss-Krond rule of 31 points is used
0138           This is used by the plug-in manager (need a char * instead of enumerations)
0139 
0140          @param type type of integration. The possible types are defined in the Integration::Type enumeration
0141          @param rule Gauss-Kronrod rule (from 1 to 6)
0142          @param absTol desired absolute Error
0143          @param relTol desired relative Error
0144          @param size maximum number of sub-intervals
0145          */
0146       GSLIntegrator(const char *  type, int rule, double absTol, double relTol, size_t size );
0147 
0148       ~GSLIntegrator() override;
0149       //~GSLIntegrator();
0150 
0151       // disable copy ctrs
0152    private:
0153 
0154       GSLIntegrator(const GSLIntegrator &);
0155       GSLIntegrator & operator=(const GSLIntegrator &);
0156 
0157    public:
0158 
0159 
0160          // template methods for generic functors
0161 
0162          /**
0163          method to set the a generic integration function
0164 
0165           @param f integration function. The function type must implement the assignment operator, <em>  double  operator() (  double  x ) </em>
0166 
0167           */
0168 
0169 
0170       void SetFunction(const IGenFunction &f) override;
0171 
0172       /**
0173          Set function from a GSL pointer function type
0174        */
0175       void SetFunction( GSLFuncPointer f, void * p = nullptr);
0176 
0177       // methods using IGenFunction
0178 
0179       /**
0180          evaluate the Integral of a function f over the defined interval (a,b)
0181        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0182        @param a lower value of the integration interval
0183        @param b upper value of the integration interval
0184        */
0185 
0186       double Integral(const IGenFunction & f, double a, double b);
0187 
0188 
0189       /**
0190          evaluate the Integral of a function f over the infinite interval (-inf,+inf)
0191        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0192        */
0193 
0194       double Integral(const IGenFunction & f);
0195 
0196       /**
0197        evaluate the Cauchy principal value of the integral of  a previously defined function f over
0198         the defined interval (a,b) with a singularity at c
0199         @param a lower interval value
0200         @param b lower interval value
0201         @param c singular value of f
0202         */
0203       double IntegralCauchy(double a, double b, double c) override;
0204 
0205       /**
0206        evaluate the Cauchy principal value of the integral of  a function f over the defined interval (a,b)
0207         with a singularity at c
0208         @param f integration function. The function type must implement the mathlib::IGenFunction interface
0209         @param a lower interval value
0210         @param b lower interval value
0211         @param c singular value of f
0212       */
0213       double IntegralCauchy(const IGenFunction & f, double a, double b, double c);
0214 
0215       /**
0216          evaluate the Integral of a function f over the semi-infinite interval (a,+inf)
0217        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0218        @param a lower value of the integration interval
0219 
0220        */
0221       double IntegralUp(const IGenFunction & f, double a );
0222 
0223       /**
0224          evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b)
0225        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0226        @param b upper value of the integration interval
0227        */
0228       double IntegralLow(const IGenFunction & f, double b );
0229 
0230       /**
0231          evaluate the Integral of a function f with known singular points over the defined Integral (a,b)
0232        @param f integration function. The function type must implement the mathlib::IGenFunction interface
0233        @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
0234 
0235        */
0236       double Integral(const IGenFunction & f, const std::vector<double> & pts );
0237 
0238       // evaluate using cached function
0239 
0240       /**
0241          evaluate the Integral over the defined interval (a,b) using the function previously set with GSLIntegrator::SetFunction method
0242        @param a lower value of the integration interval
0243        @param b upper value of the integration interval
0244        */
0245 
0246       double Integral(double a, double b) override;
0247 
0248 
0249       /**
0250          evaluate the Integral over the infinite interval (-inf,+inf) using the function previously set with GSLIntegrator::SetFunction method.
0251        */
0252       double Integral( ) override;
0253 
0254       /**
0255          evaluate the Integral of a function f over the semi-infinite interval (a,+inf) using the function previously set with GSLIntegrator::SetFunction method.
0256        @param a lower value of the integration interval
0257        */
0258       double IntegralUp(double a ) override;
0259 
0260       /**
0261          evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b) using the function previously set with GSLIntegrator::SetFunction method.
0262        @param b upper value of the integration interval
0263        */
0264       double IntegralLow( double b ) override;
0265 
0266       /**
0267          evaluate the Integral over the defined interval (a,b) using the function previously set with GSLIntegrator::SetFunction method. The function has known singular points.
0268        @param pts vector containing both the function singular points and the lower/upper edges of the interval. The vector must have as first element the lower edge of the integration Integral ( \a a) and last element the upper value.
0269 
0270        */
0271       double Integral( const std::vector<double> & pts) override;
0272 
0273       // evaluate using free function pointer (same GSL signature)
0274 
0275       /**
0276          signature for function pointers used by GSL
0277        */
0278       //typedef double ( * GSLFuncPointer ) ( double, void * );
0279 
0280       /**
0281          evaluate the Integral of  of a function f over the defined interval (a,b) passing a free function pointer
0282        The integration function must be a free function and have a signature consistent with GSL functions:
0283 
0284        <em>double my_function ( double x, void * p ) { ...... } </em>
0285 
0286        This method is the most efficient since no internal adapter to GSL function is created.
0287        @param f pointer to the integration function
0288        @param p pointer to the Parameters of the function
0289        @param a lower value of the integration interval
0290        @param b upper value of the integration interval
0291 
0292        */
0293       double Integral(GSLFuncPointer f, void * p, double a, double b);
0294 
0295       /**
0296          evaluate the Integral  of a function f over the infinite interval (-inf,+inf) passing a free function pointer
0297        */
0298       double Integral(GSLFuncPointer f, void * p);
0299 
0300       /**
0301          evaluate the Integral of a function f over the semi-infinite interval (a,+inf) passing a free function pointer
0302        */
0303       double IntegralUp(GSLFuncPointer f, void * p, double a);
0304 
0305       /**
0306          evaluate the Integral of a function f over the over the semi-infinite interval (-inf,b) passing a free function pointer
0307        */
0308       double IntegralLow(GSLFuncPointer f, void * p, double b);
0309 
0310       /**
0311          evaluate the Integral of a function f with knows singular points over the over a defined interval passing a free function pointer
0312        */
0313       double Integral(GSLFuncPointer f, void * p, const std::vector<double> & pts);
0314 
0315       /**
0316          return  the Result of the last Integral calculation
0317        */
0318       double Result() const override;
0319 
0320       /**
0321          return the estimate of the absolute Error of the last Integral calculation
0322        */
0323       double Error() const override;
0324 
0325       /**
0326          return the Error Status of the last Integral calculation
0327        */
0328       int Status() const override;
0329 
0330       /**
0331           return number of function evaluations in calculating the integral
0332       */
0333       int NEval() const override { return fNEval; }
0334 
0335       // setter for control Parameters  (getters are not needed so far )
0336 
0337       /**
0338          set the desired relative Error
0339        */
0340       void SetRelTolerance(double relTolerance) override;
0341 
0342 
0343       /**
0344          set the desired absolute Error
0345        */
0346       void SetAbsTolerance(double absTolerance) override;
0347 
0348       /**
0349          set the integration rule (Gauss-Kronrod rule).
0350        The possible rules are defined in the Integration::GKRule enumeration.
0351        The integration rule can be modified only for ADAPTIVE type integrations
0352        */
0353       void SetIntegrationRule(Integration::GKRule );
0354 
0355       /// set the options
0356       void SetOptions(const ROOT::Math::IntegratorOneDimOptions & opt) override;
0357 
0358       ///  get the option used for the integration
0359       ROOT::Math::IntegratorOneDimOptions Options() const override;
0360 
0361       /// get type name
0362       IntegrationOneDim::Type GetType() const { return fType; }
0363 
0364       /**
0365           return the name
0366       */
0367       const char * GetTypeName() const;
0368 
0369 
0370    protected:
0371 
0372       // internal method to check validity of GSL function pointer
0373       bool CheckFunction();
0374 
0375 
0376    private:
0377 
0378       Integration::Type fType;
0379       Integration::GKRule fRule;
0380       double fAbsTol;
0381       double fRelTol;
0382       size_t fSize;
0383       size_t fMaxIntervals;
0384 
0385       // cache Error, Result and Status of integration
0386 
0387       double fResult;
0388       double fError;
0389       int fStatus;
0390       int fNEval;
0391 
0392       // GSLIntegrationAlgorithm * fAlgorithm;
0393 
0394       GSLFunctionWrapper  *     fFunction;
0395       GSLIntegrationWorkspace * fWorkspace;
0396 
0397    };
0398 
0399 
0400 
0401 
0402 
0403 } // namespace Math
0404 } // namespace ROOT
0405 
0406 
0407 #endif /* ROOT_Math_GSLIntegrator */