Gaudi Framework, version v21r11
Home   Generated: 30 Sep 2010

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral Class Reference

This class allows the numerical evaluation of the following functions:. More...

#include <NumericalDefiniteIntegral.h>

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalDefiniteIntegral:
[legend]

List of all members.

Public Types

typedef std::vector< double > Points
 typedef for vector of singular points

Public Member Functions

 NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const double b, const GaudiMath::Integration::Type type=GaudiMath::Integration::Adaptive, const GaudiMath::Integration::KronrodRule rule=GaudiMath::Integration::Default, const double epsabs=1.e-10, const double epsrel=1.e-7, const size_t size=1000)
 From CLHEP/GenericFunctions.
 NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const double b, const Points &points, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
 Standard constructor The function created with this constructor compute the following integral:.
 NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const GaudiMath::Integration::Inf b=GaudiMath::Integration::Infinity, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
 Standard constructor The function created with this constructor compute the following integral:.
 NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const GaudiMath::Integration::Inf a, const double b, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
 Standard constructor The function created with this constructor compute the following integral:.
 NumericalDefiniteIntegral (const AbsFunction &function, const size_t index, const float epsabs=1e-9, const float epsrel=1.e-6, const size_t size=1000)
 Standard constructor The function created with this constructor compute the following integral:.
 NumericalDefiniteIntegral (const NumericalDefiniteIntegral &)
 copy constructor
virtual ~NumericalDefiniteIntegral ()
 destructor
virtual unsigned int dimensionality () const
 dimensionality of the problem
virtual double operator() (double argument) const
 Function value.
virtual double operator() (const Argument &argument) const
 Function value.
virtual bool hasAnalyticDerivative () const
 Does this function have an analytic derivative?
virtual Genfun::Derivative partial (unsigned int index) const
 Derivatives.
const AbsFunction & function () const
 accessor to the function itself
double a () const
 integration limit
double b () const
const Points & points () const
 known singularities
double epsabs () const
 absolute precision
double epsrel () const
 relatiove precision
double result () const
 previous result
double error () const
 evaluate of previous error
size_t size () const
GaudiMath::Integration::Type type () const
 integration type
GaudiMath::Integration::Category category () const
 integration category
GaudiMath::Integration::KronrodRule rule () const
 integration rule

Protected Member Functions

double QAGI (_Function *fun) const
double QAGP (_Function *fun) const
double QNG (_Function *fun) const
double QAG (_Function *fun) const
double QAGS (_Function *fun) const
_Workspace * allocate () const
 allocate the integration workspace
_Workspace * ws () const
StatusCode Exception (const std::string &message, const StatusCode &sc=StatusCode::FAILURE) const

Private Member Functions

 NumericalDefiniteIntegral ()
NumericalDefiniteIntegral & operator= (const NumericalDefiniteIntegral &)

Private Attributes

const AbsFunction * m_function
size_t m_DIM
size_t m_index
double m_a
double m_b
bool m_ia
bool m_ib
GaudiMath::Integration::Type m_type
GaudiMath::Integration::Category m_category
GaudiMath::Integration::KronrodRule m_rule
Points m_points
double * m_pdata
double m_epsabs
double m_epsrel
double m_result
double m_error
size_t m_size
_Workspace * m_ws
Argument m_argument
Argument m_argF

Classes

struct  _Function
struct  _Workspace

Detailed Description

This class allows the numerical evaluation of the following functions:.

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

Author:
Vanya BELYAEV Ivan.Belyaev@itep.ru
Date:
2003-08-31

Definition at line 72 of file NumericalDefiniteIntegral.h.

Member Typedef Documentation

typedef std::vector<double> Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Points

typedef for vector of singular points

Definition at line 76 of file NumericalDefiniteIntegral.h.

Constructor & Destructor Documentation

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const AbsFunction &  function,
const size_t  index,
const double  a,
const double  b,
const GaudiMath::Integration::Type  type = GaudiMath::Integration::Adaptive,
const GaudiMath::Integration::KronrodRule  rule = GaudiMath::Integration::Default,
const double  epsabs = 1.e-10,
const double  epsrel = 1.e-7,
const size_t  size = 1000 
)

From CLHEP/GenericFunctions.

from CLHEP/GenericFunctions

Standard constructor The function created with this constructor compute the following integral:

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

If function contains singularities, the type = Type::AdaptiveSingular need to be used

For faster integration of smooth function non-adaptive integration can be used: type = Type::NonAdaptive need to be used

For adaptive integration type = Type>>Adaptive one can specify the order of Gauss-Kronrad integration rule rule = KronrodRule::Gauss15 The higher-order rules give better accuracy for smooth functions, while lower-order rules save the time when the function contains local difficulties, such as discontinuites.

  • The GSL routine gsl_integration_qng is used for type = Type:NonAdaptive :
  • The GSL routine gsl_integration_qag is used for type = Type:Adaptive :
  • The GSL routine gsl_integration_qags is used for type = Type:AdaptiveSingular :

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b high integration limit
type integration type
rule integration rule (for adaptive integration)
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration
Standard constructor The function created with this constructor compute the following integral:

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

If function contains singularities, the type = Type::AdaptiveSingular need to be used

For faster integration of smooth function non-adaptive integration can be used: type = Type::NonAdaptive need to be used

For adaptive integration type = Type>>Adaptive one can specify the order of Gauss-Kronrad integration rule rule = KronrodRule::Gauss15 The higher-order rules give better accuracy for smooth functions, while lower-order rules save the time when the function contains local difficulties, such as discontinuites.

  • The GSL routine gsl_integration_qng is used for type = Type:NonAdaptive :
  • The GSL routine gsl_integration_qag is used for type = Type:Adaptive :
  • The GSL routine gsl_integration_qags is used for type = Type:AdaptiveSingular :

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b high integration limit
type integration type
rule integration rule (for adaptive integration)
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

Definition at line 104 of file NumericalDefiniteIntegral.cpp.

00113       : AbsFunction () 
00114       , m_function  ( function.clone          ()    ) 
00115       , m_DIM       ( 0                             )
00116       , m_index     ( index                         )
00117       , m_a         ( a      ) 
00118       , m_b         ( b      )
00119       , m_ia        ( false  ) 
00120       , m_ib        ( false  )
00121       , m_type      ( type   ) 
00122       , m_category  ( GaudiMath::Integration::Finite )
00123       , m_rule      ( rule       ) 
00124       , m_points    (            ) 
00125       , m_pdata     ( 0          ) 
00126       , m_epsabs    ( epsabs     )
00127       , m_epsrel    ( epsrel     )
00128       , m_result    ( GSL_NEGINF ) 
00129       , m_error     ( GSL_POSINF ) 
00130       , m_size      ( size       ) 
00131       , m_ws        ()
00132       , m_argument  ()
00133       , m_argF      () 
00134     {
00135       if ( GaudiMath::Integration::Fixed == m_rule ) 
00136         { m_rule = GaudiMath::Integration::Default ; }
00137       if ( function.dimensionality() < 2          ) 
00138         { Exception("::constructor: invalid dimensionality ") ; }        
00139       if ( m_index >= function.dimensionality()   ) 
00140         { Exception("::constructor: invalid variable index") ; }
00141 
00142       m_DIM = function.dimensionality() - 1 ;
00143       m_argument = Argument( m_DIM      ) ;
00144       m_argF     = Argument( m_DIM + 1  ) ;
00145 
00146     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const AbsFunction &  function,
const size_t  index,
const double  a,
const double  b,
const Points &  points,
const double  epsabs = 1e-9,
const double  epsrel = 1.e-6,
const size_t  size = 1000 
)

Standard constructor The function created with this constructor compute the following integral:.

standard constructor

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

where function have known singulatities

  • The GSL routine gsl_integration_qagp is used for integration

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b high integration limit
pointns vector of know singular points
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration
function the base function
index the variable index
a integration limit
b integration limit
points list of known function singularities
epsabs absolute precision for integration
epsrel relative precision for integration

Definition at line 159 of file NumericalDefiniteIntegral.cpp.

00167       : AbsFunction () 
00168       , m_function  ( function.clone()             )
00169       , m_DIM       ( 0                            ) 
00170       , m_index     ( index                        ) 
00171       , m_a         ( a                            )
00172       , m_b         ( b                            )
00173       , m_ia        ( false  ) 
00174       , m_ib        ( false  )
00175       , m_type      ( GaudiMath::Integration::    Other )
00176       , m_category  ( GaudiMath::Integration:: Singular )
00177       , m_rule      ( GaudiMath::Integration::    Fixed )
00178       , m_points    ( points  ) 
00179       , m_pdata     ( 0       ) 
00180       , m_epsabs    ( epsabs  )
00181       , m_epsrel    ( epsrel  ) 
00182       //
00183       , m_result    ( GSL_NEGINF                          ) 
00184       , m_error     ( GSL_POSINF                          ) 
00185       //
00186       , m_size      ( size                                ) 
00187       , m_ws        ( 0                                   ) 
00188       , m_argument  () 
00189       , m_argF      () 
00190     {
00191       if ( function.dimensionality() < 2          ) 
00192         { Exception("::constructor: invalid dimensionality ") ; }        
00193       if ( m_index >= function.dimensionality()   ) 
00194         { Exception("::constructor: invalid variable index") ; }
00195 
00196       m_DIM = function.dimensionality() - 1 ;
00197       m_argument = Argument( m_DIM      ) ;
00198       m_argF     = Argument( m_DIM + 1  ) ;
00199 
00200       const double l1 = std::min ( a , b ) ;
00201       const double l2 = std::max ( a , b ) ;
00202       m_points.push_back ( l1 ) ;
00203       m_points.push_back ( l2 ) ;
00204       std::sort ( m_points.begin() , m_points.end() ) ;
00205       m_points.erase( std::unique( m_points.begin () , 
00206                                    m_points.end   () ) , 
00207                       m_points.end() );
00208       
00209       Points::iterator lower = 
00210         std::lower_bound ( m_points.begin () , m_points.end () , l1 ) ;
00211       m_points.erase     ( m_points.begin () , lower                ) ;
00212       Points::iterator upper = 
00213         std::upper_bound ( m_points.begin () , m_points.end () , l2 ) ;
00214       m_points.erase     ( upper             , m_points.end ()      ) ;
00215       
00216       m_pdata = new double[ m_points.size() ] ;
00217       std::copy( m_points.begin() , m_points.end() , m_pdata );
00218     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const AbsFunction &  function,
const size_t  index,
const double  a,
const GaudiMath::Integration::Inf  b = GaudiMath::Integration::Infinity,
const double  epsabs = 1e-9,
const double  epsrel = 1.e-6,
const size_t  size = 1000 
)

Standard constructor The function created with this constructor compute the following integral:.

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagiu is used for integration

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b indication for infinity
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagiu is used for integration

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b indication for infinity
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

Definition at line 246 of file NumericalDefiniteIntegral.cpp.

00253       : AbsFunction() 
00254       , m_function  ( function.clone()             )
00255       , m_DIM       ( 0                            ) 
00256       , m_index     ( index                        ) 
00257       , m_a         ( a                            )
00258       , m_b         ( GSL_POSINF                   )
00259       , m_ia        ( false  ) 
00260       , m_ib        ( true   )
00261       , m_type      ( GaudiMath::Integration::    Other )
00262       , m_category  ( GaudiMath::Integration:: Infinite )
00263       , m_rule      ( GaudiMath::Integration::    Fixed )
00264       , m_points    (         ) 
00265       , m_pdata     ( 0       ) 
00266       , m_epsabs    ( epsabs  )
00267       , m_epsrel    ( epsrel  ) 
00268       //
00269       , m_result    ( GSL_NEGINF                          ) 
00270       , m_error     ( GSL_POSINF                          ) 
00271       //
00272       , m_size      ( size                                ) 
00273       , m_ws        ( 0                                   ) 
00274       , m_argument  () 
00275       , m_argF      () 
00276     {
00277       if ( function.dimensionality() < 2          ) 
00278         { Exception("::constructor: invalid dimensionality ") ; }        
00279       if ( m_index >= function.dimensionality()   ) 
00280         { Exception("::constructor: invalid variable index") ; }
00281 
00282       m_DIM = function.dimensionality() - 1 ;
00283       m_argument = Argument( m_DIM      ) ;
00284       m_argF     = Argument( m_DIM + 1  ) ;
00285 
00286     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const AbsFunction &  function,
const size_t  index,
const GaudiMath::Integration::Inf  a,
const double  b,
const double  epsabs = 1e-9,
const double  epsrel = 1.e-6,
const size_t  size = 1000 
)

Standard constructor The function created with this constructor compute the following integral:.

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagiu is used for integration

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b indication for infinity
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{b} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagiu is used for integration

Parameters:
function funtion to be integrated
index variable index
a lower intgeration limit
b indication for infinity
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

Definition at line 315 of file NumericalDefiniteIntegral.cpp.

00322       : AbsFunction() 
00323       , m_function  ( function.clone()             )
00324       , m_DIM       ( 0                            ) 
00325       , m_index     ( index                        ) 
00326       , m_a         ( GSL_NEGINF                   )
00327       , m_b         ( b                            )
00328       , m_ia        ( true   ) 
00329       , m_ib        ( false  )
00330       , m_type      ( GaudiMath::Integration::    Other )
00331       , m_category  ( GaudiMath::Integration:: Infinite )
00332       , m_rule      ( GaudiMath::Integration::    Fixed )
00333       , m_points    (         ) 
00334       , m_pdata     ( 0       ) 
00335       , m_epsabs    ( epsabs  )
00336       , m_epsrel    ( epsrel  ) 
00337       //
00338       , m_result    ( GSL_NEGINF                          ) 
00339       , m_error     ( GSL_POSINF                          ) 
00340       //
00341       , m_size      ( size                                ) 
00342       , m_ws        ( 0                                   ) 
00343       , m_argument  () 
00344       , m_argF      () 
00345     {
00346       if ( function.dimensionality() < 2          ) 
00347         { Exception("::constructor: invalid dimensionality ") ; }        
00348       if ( m_index >= function.dimensionality()   ) 
00349         { Exception("::constructor: invalid variable index") ; }
00350 
00351       m_DIM = function.dimensionality() - 1 ;
00352       m_argument = Argument( m_DIM      ) ;
00353       m_argF     = Argument( m_DIM + 1  ) ;
00354     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const AbsFunction &  function,
const size_t  index,
const float  epsabs = 1e-9,
const float  epsrel = 1.e-6,
const size_t  size = 1000 
)

Standard constructor The function created with this constructor compute the following integral:.

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagi is used for integration

Parameters:
function funtion to be integrated
index variable index
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

\[ {\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{+\infty} f \left(x_1, \dots , x_{i-1}, \hat{x}_i , x_{i+1}, \dots , x_n \right) \, {\mathrm{d}} \hat{x}_i \]

  • The GSL routine gsl_integration_qagi is used for integration

Parameters:
function funtion to be integrated
index variable index
epsabs required absolute precision
epsrel required relative precision
size maximal number of bisections for adaptive integration

Definition at line 380 of file NumericalDefiniteIntegral.cpp.

00385       : AbsFunction() 
00386       , m_function  ( function.clone()             )
00387       , m_DIM       ( 0      ) 
00388       , m_index     ( index                        ) 
00389       , m_a         ( GSL_NEGINF                   )
00390       , m_b         ( GSL_POSINF                   )
00391       , m_ia        ( true   ) 
00392       , m_ib        ( true   )
00393       , m_type      ( GaudiMath::Integration::    Other )
00394       , m_category  ( GaudiMath::Integration:: Infinite )
00395       , m_rule      ( GaudiMath::Integration::    Fixed )
00396       , m_points    (         ) 
00397       , m_pdata     ( 0       ) 
00398       , m_epsabs    ( epsabs  )
00399       , m_epsrel    ( epsrel  ) 
00400       //
00401       , m_result    ( GSL_NEGINF                          ) 
00402       , m_error     ( GSL_POSINF                          ) 
00403       //
00404       , m_size      ( size                                ) 
00405       , m_ws        ( 0                                   ) 
00406       , m_argument  () 
00407       , m_argF      () 
00408     {
00409       if ( function.dimensionality() < 2          ) 
00410         { Exception("::constructor: invalid dimensionality ") ; }        
00411       if ( m_index >= function.dimensionality()   ) 
00412         { Exception("::constructor: invalid variable index") ; }
00413       
00414       m_DIM = function.dimensionality() - 1 ;
00415       m_argument = Argument( m_DIM      ) ;
00416       m_argF     = Argument( m_DIM + 1  ) ;
00417       
00418     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral ( const NumericalDefiniteIntegral &  right  ) 

copy constructor

Definition at line 422 of file NumericalDefiniteIntegral.cpp.

00423       : AbsFunction () 
00424       , m_function  ( right.m_function->clone() ) 
00425       , m_DIM       ( right.m_DIM      )
00426       , m_index     ( right.m_index    )
00427       , m_a         ( right.m_a        ) 
00428       , m_b         ( right.m_b        )
00429       , m_ia        ( right.m_ia       ) 
00430       , m_ib        ( right.m_ib       )
00431       , m_type      ( right.m_type     ) 
00432       , m_category  ( right.m_category )
00433       , m_rule      ( right.m_rule     ) 
00434       , m_points    ( right.m_points   ) 
00435       , m_pdata     ( 0                ) 
00436       , m_epsabs    ( right.m_epsabs   )
00437       , m_epsrel    ( right.m_epsrel   )
00438       , m_result    ( GSL_NEGINF       ) 
00439       , m_error     ( GSL_POSINF       ) 
00440       , m_size      ( right.m_size     ) 
00441       , m_ws        ( 0                )
00442       , m_argument  ( right.m_argument )
00443       , m_argF      ( right.m_argF     )
00444     {
00445       if( 0 != right.m_pdata ) 
00446         {
00447           m_pdata = new double[m_points.size()] ;
00448           std::copy( m_points.begin() , m_points.end() , m_pdata );
00449         } 
00450     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::~NumericalDefiniteIntegral (  )  [virtual]

destructor

Definition at line 452 of file NumericalDefiniteIntegral.cpp.

00453     {
00454       if( 0 != m_function ) { delete m_function ; m_function = 0 ; }
00455       if( 0 != m_pdata    ) { delete m_pdata    ; m_pdata    = 0 ; }      
00456     }

Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::NumericalDefiniteIntegral (  )  [private]

Member Function Documentation

virtual unsigned int Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::dimensionality (  )  const [inline, virtual]

dimensionality of the problem

Definition at line 304 of file NumericalDefiniteIntegral.h.

00304 { return m_DIM ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator() ( double  argument  )  const [virtual]

Function value.

evaluate the function

Definition at line 475 of file NumericalDefiniteIntegral.cpp.

00476     {
00477       // reset the result and the error  
00478       m_result = GSL_NEGINF ;
00479       m_error  = GSL_POSINF ;
00480       
00481       // check the argument 
00482       if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
00483       
00484       m_argument[0] = argument ;
00485       return (*this) ( m_argument );
00486     }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator() ( const Argument &  argument  )  const [virtual]

Function value.

evaluate the function

Definition at line 508 of file NumericalDefiniteIntegral.cpp.

00509     {
00510       // reset the result and the error  
00511       m_result = GSL_NEGINF ;
00512       m_error  = GSL_POSINF ;
00513       
00514       // check the argument 
00515       if( argument.dimension() != m_DIM ) 
00516         { Exception ( "operator(): invalid argument size " ) ; };
00517       
00518       // copy the argument 
00519       
00520       for( size_t i  = 0 ; i < m_DIM ; ++i )
00521         { 
00522           m_argument [i] = argument [i] ;
00523           const size_t j =  i < m_index ? i : i + 1 ;
00524           m_argF     [j] = argument [i] ;
00525         };
00526       
00527       // create the helper object 
00528       GSL_Helper helper( *m_function , m_argF , m_index );
00529       
00530       // use GSL to evaluate the numerical derivative 
00531       gsl_function F ;
00532       F.function = &GSL_Adaptor ;
00533       F.params   = &helper                 ;
00534       _Function F1    ;
00535       F1.fn      = &F ;
00536       
00537       if        (  GaudiMath::Integration::Infinite         == category () )
00538         { return   QAGI ( &F1 ) ; }                                // RETURN
00539       
00540       if ( m_a == m_b ) 
00541         { 
00542           m_result = 0    ; m_error  = 0    ;                      // EXACT 
00543           return m_result ;                                       // RETURN 
00544         }
00545       
00546       if        (  GaudiMath::Integration::Singular         == category () )
00547          { return   QAGP ( &F1 ) ; }                                // RETURN
00548       else if   (  GaudiMath::Integration::Finite           == category () )
00549         if      (  GaudiMath::Integration::NonAdaptive      == type     () ) 
00550           { return QNG  ( &F1 ) ; }                                // RETURN
00551         else if (  GaudiMath::Integration::Adaptive         == type     () ) 
00552           { return QAG  ( &F1 ) ; }                                // RETURN
00553         else if (  GaudiMath::Integration::AdaptiveSingular == type     () ) 
00554           { return QAGS ( &F1 ) ; }                                // RETURN
00555         else 
00556           { Exception ( "::operator(): invalid type "  ); }
00557       else 
00558         { Exception ( "::operator(): invalid category "  ); }
00559       
00560       return 0 ;
00561     }

virtual bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::hasAnalyticDerivative (  )  const [inline, virtual]

Does this function have an analytic derivative?

Definition at line 312 of file NumericalDefiniteIntegral.h.

00312 { return true ;}

Genfun::Derivative Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::partial ( unsigned int  index  )  const [virtual]

Derivatives.

Definition at line 493 of file NumericalDefiniteIntegral.cpp.

00494     {
00495       if      ( idx >= m_DIM   )
00496         { Exception ( "::partial(i): invalid variable index " ) ; };
00497       //
00498       const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
00499       return Genfun::FunctionNoop( &aux ) ;
00500     }

const AbsFunction& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::function (  )  const [inline]

accessor to the function itself

Definition at line 320 of file NumericalDefiniteIntegral.h.

00320 { return *m_function ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::a (  )  const [inline]

integration limit

Definition at line 322 of file NumericalDefiniteIntegral.h.

00322 { return         m_a ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::b (  )  const [inline]

Definition at line 323 of file NumericalDefiniteIntegral.h.

00323 { return         m_b ; }

const Points& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::points (  )  const [inline]

known singularities

Definition at line 325 of file NumericalDefiniteIntegral.h.

00325 { return    m_points ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::epsabs (  )  const [inline]

absolute precision

Definition at line 327 of file NumericalDefiniteIntegral.h.

00327 { return    m_epsabs ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::epsrel (  )  const [inline]

relatiove precision

Definition at line 329 of file NumericalDefiniteIntegral.h.

00329 { return    m_epsrel ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::result (  )  const [inline]

previous result

Definition at line 332 of file NumericalDefiniteIntegral.h.

00332 { return    m_result ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::error (  )  const [inline]

evaluate of previous error

Definition at line 334 of file NumericalDefiniteIntegral.h.

00334 { return     m_error ; }

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::size (  )  const [inline]

Definition at line 337 of file NumericalDefiniteIntegral.h.

00337 { return      m_size ; }

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::type (  )  const [inline]

integration type

Definition at line 341 of file NumericalDefiniteIntegral.h.

00341 { return      m_type ; }

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::category (  )  const [inline]

integration category

Definition at line 344 of file NumericalDefiniteIntegral.h.

00344 { return  m_category ; }

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::rule (  )  const [inline]

integration rule

Definition at line 347 of file NumericalDefiniteIntegral.h.

00347 { return      m_rule ; }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGI ( _Function *  fun  )  const [protected]

Definition at line 583 of file NumericalDefiniteIntegral.cpp.

00584     {
00585       // check the argument 
00586       if ( 0 == F    ) { Exception("::QAGI: invalid function"); }
00587       
00588       // allocate workspace 
00589       if ( 0 == ws() ) { allocate() ; }
00590       
00591       int ierror = 0 ;
00592       
00593       if ( m_ia && m_ib ) 
00594         {
00595           ierror = gsl_integration_qagi  ( F->fn                , 
00596                                            m_epsabs  , m_epsrel , 
00597                                            size ()   , ws()->ws , 
00598                                            &m_result , &m_error ) ;
00599         }
00600       else if ( m_ia )
00601         {
00602           ierror = gsl_integration_qagil ( F->fn     , m_b      , 
00603                                            m_epsabs  , m_epsrel , 
00604                                            size ()   , ws()->ws , 
00605                                            &m_result , &m_error ) ;
00606         }
00607       else if ( m_ib )
00608         {
00609           ierror = gsl_integration_qagiu ( F->fn     , m_a      , 
00610                                            m_epsabs  , m_epsrel , 
00611                                            size ()   , ws()->ws , 
00612                                            &m_result , &m_error ) ;
00613         }
00614       else 
00615         { Exception ( "::QAGI: invalid mode" ) ; };
00616       
00617       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI" , 
00618                                 __FILE__ , __LINE__ , ierror ) ;}
00619       
00620       return m_result ;
00621     }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGP ( _Function *  fun  )  const [protected]

Definition at line 627 of file NumericalDefiniteIntegral.cpp.

00628     {
00629       if( 0 == F ) { Exception("QAGP::invalid function") ; }
00630       
00631       // no known singular points ?
00632       if( points().empty() || 0 == m_pdata ) { return QAGS( F ) ; }
00633       
00634       const size_t npts = points().size();
00635       
00636       // use GSL 
00637       int ierror =
00638         gsl_integration_qagp ( F->fn                , 
00639                                m_pdata   , npts     , 
00640                                m_epsabs  , m_epsrel ,
00641                                size ()   , ws()->ws , 
00642                                &m_result , &m_error ) ;
00643       
00644       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGI " , 
00645                                 __FILE__ , __LINE__ , ierror ) ; }
00646       
00647       // the sign
00648       if ( m_a > m_b ) { m_result *= -1 ; }
00649       
00650       return m_result ;
00651     }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QNG ( _Function *  fun  )  const [protected]

Definition at line 657 of file NumericalDefiniteIntegral.cpp.

00658     {
00659       if( 0 == F ) { Exception("QNG::invalid function") ; }
00660       
00661       // integration limits 
00662       const double low  = std::min ( m_a , m_b ) ;
00663       const double high = std::max ( m_a , m_b ) ;
00664       
00665       size_t neval = 0 ;
00666       int ierror = 
00667         gsl_integration_qng ( F->fn                 ,  
00668                               low       ,      high , 
00669                               m_epsabs  ,  m_epsrel , 
00670                               &m_result , &m_error  , &neval  ) ;
00671       
00672       if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " , 
00673                                 __FILE__ , __LINE__ , ierror ) ; }
00674       
00675       // sign
00676       if ( m_a > m_b ) { m_result *= -1 ; }
00677       
00678       return m_result ;
00679     }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAG ( _Function *  fun  )  const [protected]

Definition at line 686 of file NumericalDefiniteIntegral.cpp.

00687     {
00688       if( 0 == F ) { Exception("QAG::invalid function") ; }
00689       
00690       // allocate workspace 
00691       if( 0 == ws () ) { allocate () ; }
00692       
00693       // integration limits 
00694       const double low  = std::min ( m_a , m_b ) ;
00695       const double high = std::max ( m_a , m_b ) ;
00696       
00697       int ierror = 
00698         gsl_integration_qag ( F->fn                    , 
00699                               low       ,         high , 
00700                               m_epsabs  ,     m_epsrel , 
00701                               size   () , (int) rule() , ws ()->ws , 
00702                               &m_result ,     &m_error             );
00703       
00704       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAG " , 
00705                                 __FILE__ , __LINE__ , ierror ) ; }
00706       
00707       // sign
00708       if ( m_a > m_b ) { m_result *= -1 ; }
00709       
00710       return m_result ;
00711     }

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::QAGS ( _Function *  fun  )  const [protected]

Definition at line 717 of file NumericalDefiniteIntegral.cpp.

00718     {
00719       if( 0 == F ) { Exception("QAG::invalid function") ; }
00720       
00721       if ( m_a == m_b ) 
00722         { 
00723           m_result = 0    ;
00724           m_error  = 0    ;   // EXACT !
00725           return m_result ; 
00726         }
00727       
00728       // allocate workspace 
00729       if( 0 == ws () ) { allocate () ; }
00730       
00731       // integration limits 
00732       const double low  = std::min ( m_a , m_b ) ;
00733       const double high = std::max ( m_a , m_b ) ;
00734       
00735       int ierror = 
00736         gsl_integration_qags ( F->fn                , 
00737                                low       , high     , 
00738                                m_epsabs  , m_epsrel , 
00739                                size   () , ws()->ws , 
00740                                &m_result , &m_error );
00741       
00742       if( ierror ) { gsl_error( "NumericalDefiniteIntegral::QAGS " , 
00743                                 __FILE__ , __LINE__ , ierror ) ; }
00744       
00745       // sign
00746       if ( m_a > m_b ) { m_result *= -1 ; }
00747       
00748       return m_result ;
00749     }

NumericalDefiniteIntegral::_Workspace * Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::allocate (  )  const [protected]

allocate the integration workspace

Definition at line 568 of file NumericalDefiniteIntegral.cpp.

00569     {
00570       if ( 0 != m_ws ) { return m_ws; }
00571       gsl_integration_workspace* aux = 
00572         gsl_integration_workspace_alloc( size () );
00573       if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
00574       m_ws = new _Workspace() ;
00575       m_ws->ws = aux ;
00576       return m_ws ;
00577     }

_Workspace* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::ws (  )  const [inline, protected]

Definition at line 365 of file NumericalDefiniteIntegral.h.

00366       { return m_ws ; };

StatusCode Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::Exception ( const std::string &  message,
const StatusCode &  sc = StatusCode::FAILURE 
) const [protected]

Definition at line 462 of file NumericalDefiniteIntegral.cpp.

00464     {
00465       throw GaudiException( "NumericalDefiniteIntegral::" + message , 
00466                             "*GaudiMath*" , sc ) ;
00467       return sc ;
00468     }

NumericalDefiniteIntegral& Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::operator= ( const NumericalDefiniteIntegral &   )  [private]

Member Data Documentation

const AbsFunction* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_function [private]

Definition at line 382 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_DIM [private]

Definition at line 383 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_index [private]

Definition at line 384 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_a [private]

Definition at line 386 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_b [private]

Definition at line 387 of file NumericalDefiniteIntegral.h.

bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ia [private]

Definition at line 388 of file NumericalDefiniteIntegral.h.

bool Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ib [private]

Definition at line 389 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_type [private]

Definition at line 391 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_category [private]

Definition at line 392 of file NumericalDefiniteIntegral.h.

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_rule [private]

Definition at line 393 of file NumericalDefiniteIntegral.h.

Points Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_points [private]

Definition at line 395 of file NumericalDefiniteIntegral.h.

double* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_pdata [private]

Definition at line 396 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_epsabs [private]

Definition at line 398 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_epsrel [private]

Definition at line 399 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_result [mutable, private]

Definition at line 401 of file NumericalDefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_error [mutable, private]

Definition at line 402 of file NumericalDefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_size [private]

Definition at line 404 of file NumericalDefiniteIntegral.h.

_Workspace* Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_ws [mutable, private]

Definition at line 405 of file NumericalDefiniteIntegral.h.

Argument Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_argument [mutable, private]

Definition at line 407 of file NumericalDefiniteIntegral.h.

Argument Genfun::GaudiMathImplementation::NumericalDefiniteIntegral::m_argF [mutable, private]

Definition at line 408 of file NumericalDefiniteIntegral.h.

The documentation for this class was generated from the following files:
Generated at Thu Sep 30 09:59:04 2010 for Gaudi Framework, version v21r11 by Doxygen version 1.5.6 written by Dimitri van Heesch, © 1997-2004