The simple class for numerical integrations. More...

#include <NumericalIndefiniteIntegral.h>

Inheritance diagram for Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral:

Collaboration diagram for Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral:

Classes
struct	_Function

struct	_Workspace

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

Public Member Functions
	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const GaudiMath::Integration::Type type=GaudiMath::Integration::Adaptive, const GaudiMath::Integration::KronrodRule rule=GaudiMath::Integration::Default, const double epsabs=1e-10, const double epsrel=1.e-7, const size_t size=1000)
	From CLHEP/GenericFunctions.

	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const double a, const Points &points, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	standard constructor

	NumericalIndefiniteIntegral (const AbsFunction &function, const size_t index, const GaudiMath::Integration::Limit limit=GaudiMath::Integration::VariableHighLimit, const double epsabs=1e-9, const double epsrel=1.e-6, const size_t size=1000)
	standard constructor The function, created with this constructor evaluates following indefinite integral:

	NumericalIndefiniteIntegral (const NumericalIndefiniteIntegral &)
	copy constructor

virtual	~NumericalIndefiniteIntegral ()
	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

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::Limit	limit () const
	integration limit

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
	NumericalIndefiniteIntegral ()

NumericalIndefiniteIntegral &	operator= (const NumericalIndefiniteIntegral &)

Private Attributes
const AbsFunction *	m_function

size_t	m_DIM

size_t	m_index

double	m_a

GaudiMath::Integration::Limit	m_limit

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

Detailed Description

The simple class for numerical integrations.

It allows to evaluate following indefinite integrals:

${\mathcal{F}}_i \left(x_1, \dots , x_{i-1}, x_i , x_{i+1}, \dots , x_n \right) = \int\limits_{a}^{x_i} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{x_i}^{a} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{-\infty}^{x_i} 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 , x_{i+1}, \dots , x_n \right) = \int\limits_{x_i}^{+\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..nosp@m.Bely.nosp@m.aev@i.nosp@m.tep..nosp@m.ru

Date: 2003-08-31

Definition at line 76 of file NumericalIndefiniteIntegral.h.

Member Typedef Documentation

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

typedef for vector of singular points

Definition at line 80 of file NumericalIndefiniteIntegral.h.

Constructor & Destructor Documentation

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const GaudiMath::Integration::Limit	limit = `GaudiMath::Integration::VariableHighLimit`,
		const GaudiMath::Integration::Type	type = `GaudiMath::Integration::Adaptive`,
		const GaudiMath::Integration::KronrodRule	rule = `GaudiMath::Integration::Default`,
		const double	epsabs = `1e-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 evaluates following indefinite integral:

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

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 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	the base function
index	the variable index
a	integration limit
limit	flag to distinguisch low variable limit from high variable limit
type	the integration type (adaptive, non-adaptive or adaptive for singular functions
key	Gauss-Kronrad integration rule
epsabs	absolute precision for integration
epsrel	relative precision for integration
lim	bisection limit

Standard constructor

Parameters

function	the base function
index	the variable index
a	integration limit
limit	flag to distinguisch low variable limit from high variable limit
type	the integration type (adaptive, non-adaptive or adaptive for singular functions
key	Gauss-Kronrad integration rule
epsabs	absolute precision for integration
epsrel	relative precision for integration
lim	bisection limit

Definition at line 83 of file NumericalIndefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( function.clone()                    )
      , m_DIM       ( function.dimensionality()           )
      , m_index     ( index                               )
      , m_a         ( a                                   )
      , m_limit     ( limit                               )
      , m_type      ( type                                )
      , m_category  ( GaudiMath::Integration::Finite      )
      , m_rule      ( rule                                )
      //
      , m_points    (                                     )
      , m_pdata     ( 0                                   )
      //
      , m_epsabs    ( epsabs                              )
      , m_epsrel    ( epsrel                              )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ( function.dimensionality()           )
    {
      if ( GaudiMath::Integration::Fixed == m_rule )
        { m_rule = GaudiMath::Integration::Default ; }
      if ( m_index >= m_DIM  )
        { Exception("::constructor: invalid variable index") ; }
    }

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const double	a,
		const Points &	points,
		const GaudiMath::Integration::Limit	limit = `GaudiMath::Integration::VariableHighLimit`,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

standard constructor

The function, created with this constructor evaluates following indefinite integral:

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

The integrand is assumed to have a known discontinuities

The GSL routine gsl_integration_qagp is used for integration

Parameters

function	the base function
index	the variable index
a	integration limit
limit	flag to distinguisch low variable limit from high variable limit
points	list of known function singularities
epsabs	absolute precision for integration
epsrel	relative precision for integration
function	the base function
index	the variable index
a	integration limit
limit	flag to distinguisch low variable limit from high variable limit
points	list of known function singularities
epsabs	absolute precision for integration
epsrel	relative precision for integration

Definition at line 133 of file NumericalIndefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( function.clone()          )
      , m_DIM       ( function.dimensionality() )
      , m_index     ( index                     )
      , m_a         ( a                         )
      , m_limit     ( limit                     )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Singular )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    ( points  )
      , m_pdata     ( 0       )
      , m_epsabs    ( epsabs  )
      , m_epsrel    ( epsrel  )
      //
      , m_result    ( GSL_NEGINF                          )
      , m_error     ( GSL_POSINF                          )
      //
      , m_size      ( size                                )
      , m_ws        ( 0                                   )
      , m_argument  ( function.dimensionality()           )
    {
      if ( m_index >= m_DIM )
        { Exception("::constructor: invalid variable index") ; }
      m_pdata = new double[ 2 + m_points.size() ] ;
      m_points.push_back( a ) ;
      std::sort( m_points.begin() , m_points.end() ) ;
      m_points.erase ( std::unique( m_points.begin () ,
                                    m_points.end   () ) , m_points.end() );
    }

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral	(	const AbsFunction &	function,
		const size_t	index,
		const GaudiMath::Integration::Limit	limit = `GaudiMath::Integration::VariableHighLimit`,
		const double	epsabs = `1e-9`,
		const double	epsrel = `1.e-6`,
		const size_t	size = `1000`
	)

standard constructor The function, created with this constructor evaluates following indefinite integral:

standard constructor the integral limt is assumed to be infinity

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

for value of limit = VariableHighLimit and the integral

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

for value of limit = VariableLowLimit

The GSL routines gsl_integration_qagil and gsl_integration_qagiu are used for adapive integration

Parameters

function	the base function
index	the variable index
limit	flag to distinguisch low variable limit from high variable limit
singularities	list of known function singularities
function	the base function
index	the variable index
limit	flag to distinguisch low variable limit from high variable limit

Definition at line 181 of file NumericalIndefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( function.clone() )
      , m_DIM       ( function.dimensionality() )
      , m_index     ( index            )
      , m_a         ( GSL_NEGINF       ) // should not be used!
      , m_limit     ( limit            )
      , m_type      ( GaudiMath::Integration::    Other )
      , m_category  ( GaudiMath::Integration:: Infinite )
      , m_rule      ( GaudiMath::Integration::    Fixed )
      , m_points    (            )
      , m_pdata     ( 0          )
      , m_epsabs    ( epsabs     )
      , m_epsrel    ( epsrel     )
      , m_result    ( GSL_NEGINF )
      , m_error     ( GSL_POSINF )
      , m_size      ( size       )
      , m_ws        ( 0          )
      , m_argument  ( function.dimensionality()           )
    {
      if ( m_index >= m_DIM )
        { Exception("::constructor: invalid variable index") ; }
    }

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

copy constructor

Definition at line 216 of file NumericalIndefiniteIntegral.cpp.

      : AbsFunction ()
      , m_function  ( right.m_function->clone() )
      , m_DIM       ( right.m_DIM      )
      , m_index     ( right.m_index    )
      , m_a         ( right.m_a        )
      , m_limit     ( right.m_limit    )
      , m_type      ( right.m_type     )
      , m_category  ( right.m_category )
      , m_rule      ( right.m_rule     )
      , m_points    ( right.m_points   )
      , m_pdata     ( 0 )           // attention
      , m_epsabs    ( right.m_epsabs   )
      , m_epsrel    ( right.m_epsrel   )
      , m_result    ( GSL_NEGINF       )
      , m_error     ( GSL_POSINF       )
      , m_size      ( right.m_size     )
      , m_ws        ( 0                )
      , m_argument  ( right.m_argument )
    {
      m_pdata = new double[ 2 + m_points.size() ] ; // attention!
    }

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::~NumericalIndefiniteIntegral ( )

virtual

destructor

Definition at line 244 of file NumericalIndefiniteIntegral.cpp.

    {
      if( 0 != m_ws )
        {
          gsl_integration_workspace_free ( m_ws->ws ) ;
          delete m_ws ;
          m_ws = 0 ;
        }
      if ( 0 != m_pdata    ) { delete m_pdata    ; m_pdata    = 0 ; }
      if ( 0 != m_function ) { delete m_function ; m_function = 0 ; }
    }

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::NumericalIndefiniteIntegral ( )

private

Member Function Documentation

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::a ( ) const

inline

integration limit

Definition at line 296 of file NumericalIndefiniteIntegral.h.

{ return m_a ; }

NumericalIndefiniteIntegral::_Workspace * Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::allocate ( ) const

protected

allocate the integration workspace

Definition at line 361 of file NumericalIndefiniteIntegral.cpp.

    {
      if ( 0 != m_ws ) { return m_ws; }
      gsl_integration_workspace* aux =
        gsl_integration_workspace_alloc( size () );
      if ( 0 == aux ) { Exception ( "allocate()::invalid workspace" ) ; };
      m_ws = new _Workspace() ;
      m_ws->ws = aux ;
      return m_ws ;
    }

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::category ( ) const

inline

integration category

Definition at line 320 of file NumericalIndefiniteIntegral.h.

{ return m_category ; }

virtual unsigned int Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::dimensionality ( ) const

inlinevirtual

dimensionality of the problem

Definition at line 278 of file NumericalIndefiniteIntegral.h.

{ return m_DIM ; }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::epsabs ( ) const

inline

absolute precision

Definition at line 300 of file NumericalIndefiniteIntegral.h.

{ return m_epsabs ; }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::epsrel ( ) const

inline

relatiove precision

Definition at line 302 of file NumericalIndefiniteIntegral.h.

{ return m_epsrel ; }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::error ( ) const

inline

evaluate of previous error

Definition at line 307 of file NumericalIndefiniteIntegral.h.

{ return m_error ; }

StatusCode Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::Exception	(	const std::string &	message,
		const StatusCode &	sc = `StatusCode::FAILURE`
	)		const

protected

Definition at line 261 of file NumericalIndefiniteIntegral.cpp.

    {
      throw GaudiException( "NumericalIndefiniteIntegral::" + message ,
                            "*GaudiMath*" , sc ) ;
      return sc ;
    }

const AbsFunction& Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::function ( ) const

inline

accessor to the function itself

Definition at line 294 of file NumericalIndefiniteIntegral.h.

{ return *m_function ; }

virtual bool Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::hasAnalyticDerivative ( ) const

inlinevirtual

Does this function have an analytic derivative?

Definition at line 286 of file NumericalIndefiniteIntegral.h.

{ return true ;}

GaudiMath::Integration::Limit Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::limit ( ) const

inline

integration limit

Definition at line 314 of file NumericalIndefiniteIntegral.h.

{ return m_limit ; }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::operator() ( double argument ) const

virtual

Function value.

evaluate the function

Definition at line 274 of file NumericalIndefiniteIntegral.cpp.

    {
      // reset the result and the error
      m_result = GSL_NEGINF ;
      m_error  = GSL_POSINF ;
      // check the argument
      if( 1 != m_DIM ) { Exception ( "operator(): invalid argument size " ) ; };
      m_argument[0] = argument ;
      return (*this) ( m_argument );
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::operator() ( const Argument & argument ) const

virtual

Function value.

evaluate the function

Definition at line 314 of file NumericalIndefiniteIntegral.cpp.

    {
      // reset the result and the error
      m_result = GSL_NEGINF ;
      m_error  = GSL_POSINF ;
      // check the argument
      if( argument.dimension() != m_DIM )
        { Exception ( "operator(): invalid argument size " ) ; };
      // copy the argument
      {for( size_t i  = 0 ; i < m_DIM ; ++i ){ m_argument[i] = argument[i];}}
      // create the helper object
      GSL_Helper helper( *m_function , m_argument , m_index );
      // use GSL to evaluate the numerical derivative
      gsl_function F ;
      F.function = &GSL_Adaptor ;
      F.params   = &helper                 ;
      _Function F1    ;
      F1.fn      = &F ;
      if        (  GaudiMath::Integration::Infinite         == category () )
        { return   QAGI ( &F1 ) ; }                                // RETURN
      else if   (  GaudiMath::Integration::Singular         == category () )
        { return   QAGP ( &F1 ) ; }                                // RETURN
      else if   (  GaudiMath::Integration::Finite           == category () )
        if      (  GaudiMath::Integration::NonAdaptive      == type     () )
          { return QNG  ( &F1 ) ; }                                // RETURN
        else if (  GaudiMath::Integration::Adaptive         == type     () )
          { return QAG  ( &F1 ) ; }                                // RETURN
        else if (  GaudiMath::Integration::AdaptiveSingular == type     () )
          { return QAGS ( &F1 ) ; }                                // RETURN
        else
          { Exception ( "::operator(): invalid type "  ); }
      else
        { Exception ( "::operator(): invalid category "  ); }
      return 0 ;
    }

NumericalIndefiniteIntegral& Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::operator= ( const NumericalIndefiniteIntegral & )

private

Genfun::Derivative Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::partial ( unsigned int index ) const

virtual

Derivatives.

Definition at line 292 of file NumericalIndefiniteIntegral.cpp.

    {
      if      ( idx >= m_DIM   )
        { Exception ( "::partial(i): invalid variable index " ) ; };
      if      ( idx != m_index )
        {
          const AbsFunction& aux = NumericalDerivative( *this , idx ) ;
          return Genfun::FunctionNoop( &aux ) ;
        }
      else if ( GaudiMath::Integration::VariableLowLimit == limit () )
        {
          const AbsFunction& aux = -1 * function() ;
          return Genfun::FunctionNoop( &aux ) ;
        }
      const AbsFunction& aux = function() ;
      return Genfun::FunctionNoop( &aux ) ;
    }

const Points& Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::points ( ) const

inline

known singularities

Definition at line 298 of file NumericalIndefiniteIntegral.h.

{ return m_points ; }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::QAG ( _Function * fun ) const

protected

Definition at line 507 of file NumericalIndefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAG::invalid function") ; }
      const double x = m_argument[m_index] ;
      if ( m_a == x  )
        {
          m_result = 0    ;
          m_error  = 0    ;   // EXACT !
          return m_result ;
        }
      // allocate workspace
      if( 0 == ws () ) { allocate () ; }
      // integration limits
      const double a = std::min ( m_a , x ) ;
      const double b = std::max ( m_a , x ) ;
      int ierror =
        gsl_integration_qag ( F->fn                    ,
                              a         ,            b ,
                              m_epsabs  ,     m_epsrel ,
                              size   () , (int) rule() , ws ()->ws ,
                              &m_result ,     &m_error             );
      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAG " ,
                                __FILE__ , __LINE__ , ierror ) ; }
      // sign
      if      ( GaudiMath::Integration::VariableHighLimit == limit()
                &&  x < m_a  ) { m_result *= -1 ; }
      else if ( GaudiMath::Integration::VariableLowLimit  == limit()
                &&  x > m_a  ) { m_result *= -1 ; }
      return m_result ;
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::QAGI ( _Function * fun ) const

protected

Definition at line 376 of file NumericalIndefiniteIntegral.cpp.

    {
      // check the argument
      if( 0 == F ) { Exception("::QAGI: invalid function"); }
      const double x = m_argument[m_index] ;
      // allocate workspace
      if( 0 == ws() ) { allocate() ; }
      int ierror = 0 ;
      switch ( limit() )
        {
        case GaudiMath::Integration::VariableLowLimit  :
          ierror = gsl_integration_qagiu ( F->fn     , x        ,
                                           m_epsabs  , m_epsrel ,
                                           size ()   , ws()->ws ,
                                           &m_result , &m_error ) ; break ;
        case GaudiMath::Integration::VariableHighLimit :
          ierror = gsl_integration_qagil ( F->fn     , x        ,
                                           m_epsabs  , m_epsrel ,
                                           size ()   , ws()->ws ,
                                           &m_result , &m_error ) ; break ;
        default :
          Exception ( "::QAGI: invalid mode" ) ;
        };
      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI" ,
                                __FILE__ , __LINE__ , ierror ) ;}
      return m_result ;
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::QAGP ( _Function * fun ) const

protected

Definition at line 413 of file NumericalIndefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAGP::invalid function") ; }
      const double x = m_argument[m_index] ;
      if ( m_a == x  )
        {
          m_result = 0    ;
          m_error  = 0    ;   // EXACT !
          return m_result ;
        }
      // no known singular points ?
      if( points().empty() ) { return QAGS( F ) ; }
      // integration limits
      const double a = std::min ( m_a , x ) ;
      const double b = std::max ( m_a , x ) ;
      // "active" singular points
      Points::const_iterator lower =
        std::lower_bound ( points().begin() , points().end() , a ) ;
      Points::const_iterator upper =
        std::upper_bound ( points().begin() , points().end() , b ) ;
      Points pnts ( upper - lower ) ;
      std::copy( lower , upper , pnts.begin() );
      if ( *lower != a       ) { pnts.insert( pnts.begin () , a ) ; }
      if ( *upper != b       ) { pnts.insert( pnts.end   () , b ) ; }
      std::copy( pnts.begin() , pnts.end() , m_pdata );
      const size_t npts = pnts.size() ;
      // use GSL
      int ierror =
        gsl_integration_qagp ( F->fn                ,
                               m_pdata   , npts     ,
                               m_epsabs  , m_epsrel ,
                               size ()   , ws()->ws ,
                               &m_result , &m_error ) ;
      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGI " ,
                                __FILE__ , __LINE__ , ierror ) ; }
      // sign
      if      ( GaudiMath::Integration::VariableHighLimit == limit()
                &&  x < m_a  ) { m_result *= -1 ; }
      else if ( GaudiMath::Integration::VariableLowLimit  == limit()
                &&  x > m_a  ) { m_result *= -1 ; }
      return m_result ;
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::QAGS ( _Function * fun ) const

protected

Definition at line 549 of file NumericalIndefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QAG::invalid function") ; }
      const double x = m_argument[m_index] ;
      if ( m_a == x  )
        {
          m_result = 0    ;
          m_error  = 0    ;   // EXACT !
          return m_result ;
        }
      // allocate workspace
      if( 0 == ws () ) { allocate () ; }
      // integration limits
      const double a = std::min ( m_a , x ) ;
      const double b = std::max ( m_a , x ) ;
      int ierror =
        gsl_integration_qags ( F->fn                ,
                               a         , b        ,
                               m_epsabs  , m_epsrel ,
                               size   () , ws()->ws ,
                               &m_result , &m_error );
      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QAGS " ,
                                __FILE__ , __LINE__ , ierror ) ; }
      // sign
      if      ( GaudiMath::Integration::VariableHighLimit == limit()
                &&  x < m_a  ) { m_result *= -1 ; }
      else if ( GaudiMath::Integration::VariableLowLimit  == limit()
                &&  x > m_a  ) { m_result *= -1 ; }
      return m_result ;
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::QNG ( _Function * fun ) const

protected

Definition at line 469 of file NumericalIndefiniteIntegral.cpp.

    {
      if( 0 == F ) { Exception("QNG::invalid function") ; }
      const double x = m_argument[m_index] ;
      if ( m_a == x  )
        {
          m_result = 0    ;
          m_error  = 0    ;   // EXACT !
          return m_result ;
        }
      // integration limits
      const double a = std::min ( m_a , x ) ;
      const double b = std::max ( m_a , x ) ;
      size_t neval = 0 ;
      int ierror =
        gsl_integration_qng ( F->fn                 ,
                              a         ,         b ,
                              m_epsabs  ,  m_epsrel ,
                              &m_result , &m_error  , &neval  ) ;
      if( ierror ) { gsl_error( "NumericalIndefiniteIntegral::QNG " ,
                                __FILE__ , __LINE__ , ierror ) ; }
      // sign
      if      ( GaudiMath::Integration::VariableHighLimit == limit()
                &&  x < m_a  ) { m_result *= -1 ; }
      else if ( GaudiMath::Integration::VariableLowLimit  == limit()
                &&  x > m_a  ) { m_result *= -1 ; }
      return m_result ;
    }

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::result ( ) const

inline

previous result

Definition at line 305 of file NumericalIndefiniteIntegral.h.

{ return m_result ; }

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::rule ( ) const

inline

integration rule

Definition at line 323 of file NumericalIndefiniteIntegral.h.

{ return m_rule ; }

size_t Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::size ( ) const

inline

Definition at line 310 of file NumericalIndefiniteIntegral.h.

{ return m_size ; }

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::type ( ) const

inline

integration type

Definition at line 317 of file NumericalIndefiniteIntegral.h.

{ return m_type ; }

_Workspace* Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::ws ( ) const

inlineprotected

Definition at line 341 of file NumericalIndefiniteIntegral.h.

{ return m_ws ; };

Member Data Documentation

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_a

private

Definition at line 363 of file NumericalIndefiniteIntegral.h.

Argument Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_argument

mutableprivate

Definition at line 382 of file NumericalIndefiniteIntegral.h.

GaudiMath::Integration::Category Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_category

private

Definition at line 367 of file NumericalIndefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_DIM

private

Definition at line 360 of file NumericalIndefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_epsabs

private

Definition at line 373 of file NumericalIndefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_epsrel

private

Definition at line 374 of file NumericalIndefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_error

mutableprivate

Definition at line 377 of file NumericalIndefiniteIntegral.h.

const AbsFunction* Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_function

private

Definition at line 359 of file NumericalIndefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_index

private

Definition at line 361 of file NumericalIndefiniteIntegral.h.

GaudiMath::Integration::Limit Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_limit

private

Definition at line 365 of file NumericalIndefiniteIntegral.h.

double* Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_pdata

private

Definition at line 371 of file NumericalIndefiniteIntegral.h.

Points Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_points

private

Definition at line 370 of file NumericalIndefiniteIntegral.h.

double Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_result

mutableprivate

Definition at line 376 of file NumericalIndefiniteIntegral.h.

GaudiMath::Integration::KronrodRule Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_rule

private

Definition at line 368 of file NumericalIndefiniteIntegral.h.

size_t Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_size

private

Definition at line 379 of file NumericalIndefiniteIntegral.h.

GaudiMath::Integration::Type Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_type

private

Definition at line 366 of file NumericalIndefiniteIntegral.h.

_Workspace* Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral::m_ws

mutableprivate

Definition at line 380 of file NumericalIndefiniteIntegral.h.

The documentation for this class was generated from the following files:

/afs/cern.ch/sw/Gaudi/releases/GAUDI/GAUDI_v23r8/GaudiGSL/GaudiMath/NumericalIndefiniteIntegral.h
/afs/cern.ch/sw/Gaudi/releases/GAUDI/GAUDI_v23r8/GaudiGSL/src/Lib/NumericalIndefiniteIntegral.cpp

Gaudi Framework, version v23r8

Classes

Public Types

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation