The simple class for numerical integrations. More...

#include <NumericalIndefiniteIntegral.h>

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.Belyaev@itep.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 [inline, virtual]

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 [inline, virtual]

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() ( 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 ;
    }

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 );
    }

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 ( void ) 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 [inline, protected]

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 [mutable, private]

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 [mutable, private]

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 [mutable, private]

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 [mutable, private]

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_v23r0/GaudiGSL/GaudiMath/NumericalIndefiniteIntegral.h
/afs/cern.ch/sw/Gaudi/releases/GAUDI/GAUDI_v23r0/GaudiGSL/src/Lib/NumericalIndefiniteIntegral.cpp

Gaudi Framework, version v23r0

Genfun::GaudiMathImplementation::NumericalIndefiniteIntegral Class Reference

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