Classes
struct	control
	User control parameters for R's integrate. More...

struct	Integral
	Interface to R's adaptive integrate routine. More...

struct	mvIntegral
	Multivariate integral class. More...

struct	mvIntegral0

Functions
template<class S , class T >
double	fmax2 (S x, T y)

template<class S , class T >
double	fmin2 (S x, T y)

template<class S , class T >
int	imin2 (S x, T y)

template<class Integrand >
Integrand::Scalar	integrate (Integrand f, typename Integrand::Scalar a=-INFINITY, typename Integrand::Scalar b=INFINITY, control c=control())
	Integrate function over finite or infinite interval. More...

template<class Integrand >
mvIntegral0< Integrand >	mvIntegrate (Integrand &f, control c=control())
	Multivariate integration. More...

template<class Float , class integr_fn >
void	Rdqagi (integr_fn f, void ex, Float bound, int inf, Float epsabs, Float epsrel, Float result, Float abserr, int neval, int ier, int limit, int lenw, int last, int iwork, Float work)

template<class Float , class integr_fn >
static void	rdqagie (integr_fn f, void ex, Float , int , Float , Float , int , Float , Float , int , int , Float , Float , Float , Float , int , int )

template<class Float , class integr_fn >
void	Rdqags (integr_fn f, void ex, Float a, Float b, Float epsabs, Float epsrel, Float result, Float abserr, int neval, int ier, int limit, int lenw, int last, int iwork, Float work)

template<class Float , class integr_fn >
static void	rdqagse (integr_fn f, void ex, Float , Float , Float , Float , int , Float , Float , int , int , Float , Float , Float , Float , int , int )

template<class Float >
static void	rdqelg (int , Float , Float , Float , Float , int )

template<class Float , class integr_fn >
static void	rdqk15i (integr_fn f, void ex, Float , int , Float , Float , Float , Float , Float , Float *)

template<class Float , class integr_fn >
static void	rdqk21 (integr_fn f, void ex, Float , Float , Float , Float , Float , Float *)

template<class Float >
static void	rdqpsrt (int , int , int , Float , Float , int , int *)

template<class T >
double	value (T x)

Function Documentation

template<class S , class T >

double gauss_kronrod::fmax2	(	S	x,
		T	y
	)

Definition at line 33 of file integrate.hpp.

template<class S , class T >

double gauss_kronrod::fmin2	(	S	x,
		T	y
	)

Definition at line 31 of file integrate.hpp.

Referenced by rdqk15i(), and rdqk21().

template<class S , class T >

int gauss_kronrod::imin2	(	S	x,
		T	y
	)

Definition at line 29 of file integrate.hpp.

template<class Integrand >

Integrand::Scalar gauss_kronrod::integrate	(	Integrand	f,
		typename Integrand::Scalar	a = `-INFINITY`,
		typename Integrand::Scalar	b = `INFINITY`,
		control	c = `control()`
	)

Integrate function over finite or infinite interval.

Parameters

f	Univariate integrand (functor)
a	Lower integration limit. Default is negative infinity.
a	Upper integration limit. Default is positive infinity.
c	Optional control parameters.

Example:

template<class Float>
struct Gauss_t {
  typedef Float Scalar;
  Float a;  // Parameter
  // Evaluate integrand
  Float operator(Float x) () {
    Float ans = exp(- a*x*x);
    return ans;
  }
  // Integrate wrt x
  Float my_integrate() {
    using gauss_kronrod::integrate;
    Float ans = integrate(*this);
    return ans;
  }
};

Definition at line 158 of file integrate.hpp.

template<class Integrand >

mvIntegral0<Integrand> gauss_kronrod::mvIntegrate	(	Integrand &	f,
		control	c = `control()`
	)

Multivariate integration.

Parameters

f	Multivariate integrand (functor)
c	Optional control parameters

Example:

template<class Float>
struct Gauss2D_t {
  typedef Float Scalar;
  Float a, b;  // Parameters
  Float x, y;  // Integration variables
  // Evaluate integrand (u1,u2)
  Float operator() () {
    Float ans = exp(- a*x*x - b*y*y);
    return ans;
  }
  // Integrate wrt (x,y)
  Float my_integrate() {
    using gauss_kronrod::mvIntegrate;
    Float ans = mvIntegrate(*this).
wrt(x, -INFINITY, INFINITY).
wrt(y, -INFINITY, INFINITY) ();
    return ans;
  }
};

Definition at line 254 of file integrate.hpp.

template<class Float , class integr_fn >

void gauss_kronrod::Rdqagi	(	integr_fn	f,
		void *	ex,
		Float *	bound,
		int *	inf,
		Float *	epsabs,
		Float *	epsrel,
		Float *	result,
		Float *	abserr,
		int *	neval,
		int *	ier,
		int *	limit,
		int *	lenw,
		int *	last,
		int *	iwork,
		Float *	work
	)

Definition at line 72 of file integrate.hpp.

Referenced by gauss_kronrod::Integral< gauss_kronrod::mvIntegral::evaluator >::operator()().

template<class Float , class integr_fn >

static void gauss_kronrod::rdqagie	(	integr_fn	f,
		void *	ex,
		Float *	bound,
		int *	inf,
		Float *	epsabs,
		Float *	epsrel,
		int *	limit,
		Float *	result,
		Float *	abserr,
		int *	neval,
		int *	ier,
		Float *	alist,
		Float *	blist,
		Float *	rlist,
		Float *	elist,
		int *	iord,
		int *	last
	)

static

begin prologue dqagie date written 800101 (yymmdd) revision date 830518 (yymmdd) category no. h2a3a1,h2a4a1 keywords automatic integrator, infinite intervals, general-purpose, transformation, extrapolation, globally adaptive author piessens,robert,appl. math & progr. div - k.u.leuven de doncker,elise,appl. math & progr. div - k.u.leuven purpose the routine calculates an approximation result to a given integral i = integral of f over (bound,+infinity) or i = integral of f over (-infinity,bound) or i = integral of f over (-infinity,+infinity), hopefully satisfying following claim for accuracy abs(i-result) <= max(epsabs,epsrel*abs(i)) description

integration over infinite intervals standard fortran subroutine

       f      - double precision
                function subprogram defining the integrand
                function f(x). the actual name for f needs to be
                declared e x t e r n a l in the driver program.

       bound  - double precision
                finite bound of integration range
                (has no meaning if interval is doubly-infinite)

       inf    - double precision
                indicating the kind of integration range involved
                inf = 1 corresponds to  (bound,+infinity),
                inf = -1            to  (-infinity,bound),
                inf = 2             to (-infinity,+infinity).

       epsabs - double precision
                absolute accuracy requested
       epsrel - double precision
                relative accuracy requested
                if  epsabs <= 0
                and epsrel < max(50*rel.mach.acc.,0.5d-28),
                the routine will end with ier = 6.

       limit  - int
                gives an upper bound on the number of subintervals
                in the partition of (a,b), limit >= 1

    on return
       result - double precision
                approximation to the integral

       abserr - double precision
                estimate of the modulus of the absolute error,
                which should equal or exceed abs(i-result)

       neval  - int
                number of integrand evaluations

       ier    - int
                ier = 0 normal and reliable termination of the
                        routine. it is assumed that the requested
                        accuracy has been achieved.
              - ier > 0 abnormal termination of the routine. the
                        estimates for result and error are less
                        reliable. it is assumed that the requested
                        accuracy has not been achieved.
       error messages
                ier = 1 maximum number of subdivisions allowed
                        has been achieved. one can allow more
                        subdivisions by increasing the value of
                        limit (and taking the according dimension
                        adjustments into account). however,if
                        this yields no improvement it is advised
                        to analyze the integrand in order to
                        determine the integration difficulties.
                        if the position of a local difficulty can
                        be determined (e.g. singularity,
                        discontinuity within the interval) one
                        will probably gain from splitting up the
                        interval at this point and calling the
                        integrator on the subranges. if possible,
                        an appropriate special-purpose integrator
                        should be used, which is designed for
                        handling the type of difficulty involved.
                    = 2 the occurrence of roundoff error is
                        detected, which prevents the requested
                        tolerance from being achieved.
                        the error may be under-estimated.
                    = 3 extremely bad integrand behaviour occurs
                        at some points of the integration
                        interval.
                    = 4 the algorithm does not converge.
                        roundoff error is detected in the
                        extrapolation table.
                        it is assumed that the requested tolerance
                        cannot be achieved, and that the returned
                        result is the best which can be obtained.
                    = 5 the integral is probably divergent, or
                        slowly convergent. it must be noted that
                        divergence can occur with any other value
                        of ier.
                    = 6 the input is invalid, because
                        (epsabs <= 0 and
                         epsrel < max(50*rel.mach.acc.,0.5d-28),
                        result, abserr, neval, last, rlist(1),
                        elist(1) and iord(1) are set to zero.
                        alist(1) and blist(1) are set to 0
                        and 1 respectively.

       alist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the left
                end points of the subintervals in the partition
                of the transformed integration range (0,1).

       blist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the right
                end points of the subintervals in the partition
                of the transformed integration range (0,1).

       rlist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the integral
                approximations on the subintervals

       elist  - double precision
                vector of dimension at least limit,  the first
                last elements of which are the moduli of the
                absolute error estimates on the subintervals

       iord   - int
                vector of dimension limit, the first k
                elements of which are pointers to the
                error estimates over the subintervals,
                such that elist(iord(1)), ..., elist(iord(k))
                form a decreasing sequence, with k = last
                if last <= (limit/2+2), and k = limit+1-last
                otherwise

       last   - int
                number of subintervals actually produced
                in the subdivision process

routines called dqelg,dqk15i,dqpsrt end prologue dqagie

       the dimension of rlist2 is determined by the value of
       limexp in subroutine dqelg.

       list of major variables
       -----------------------

      alist     - list of left end points of all subintervals
                  considered up to now
      blist     - list of right end points of all subintervals
                  considered up to now
      rlist(i)  - approximation to the integral over
                  (alist(i),blist(i))
      rlist2    - array of dimension at least (limexp+2),
                  containing the part of the epsilon table
                  wich is still needed for further computations
      elist(i)  - error estimate applying to rlist(i)
      maxerr    - pointer to the interval with largest error
                  estimate
      errmax    - elist(maxerr)
      erlast    - error on the interval currently subdivided
                  (before that subdivision has taken place)
      area      - sum of the integrals over the subintervals
      errsum    - sum of the errors over the subintervals
      errbnd    - requested accuracy max(epsabs,epsrel*
                  abs(result))
           1    - variable for the left subinterval
           2    - variable for the right subinterval
      last      - index for subdivision
      nres      - number of calls to the extrapolation routine
      numrl2    - number of elements currently in rlist2. if an
                  appropriate approximation to the compounded
                  integral has been obtained, it is put in
                  rlist2(numrl2) after numrl2 has been increased
                  by one.
      small     - length of the smallest interval considered up
                  to now, multiplied by 1.5
      erlarg    - sum of the errors over the intervals larger
                  than the smallest interval considered up to now
      extrap    - logical variable denoting that the routine
                  is attempting to perform extrapolation. i.e.
                  before subdividing the smallest interval we
                  try to decrease the value of erlarg.
      noext     - logical variable denoting that extrapolation
                  is no longer allowed (true-value)

       machine dependent constants
       ---------------------------

      epmach is the largest relative spacing.
      uflow is the smallest positive magnitude.
      oflow is the largest positive magnitude.

begin prologue dqagie date written 800101 (yymmdd) revision date 830518 (yymmdd) category no. h2a3a1,h2a4a1 keywords automatic integrator, infinite intervals, general-purpose, transformation, extrapolation, globally adaptive author piessens,robert,appl. math & progr. div - k.u.leuven de doncker,elise,appl. math & progr. div - k.u.leuven purpose the routine calculates an approximation result to a given integral i = integral of f over (bound,+infinity) or i = integral of f over (-infinity,bound) or i = integral of f over (-infinity,+infinity), hopefully satisfying following claim for accuracy abs(i-result) <= max(epsabs,epsrel*abs(i)) description

integration over infinite intervals standard fortran subroutine

       f      - double precision
                function subprogram defining the integrand
                function f(x). the actual name for f needs to be
                declared e x t e r n a l in the driver program.

       bound  - double precision
                finite bound of integration range
                (has no meaning if interval is doubly-infinite)

       inf    - double precision
                indicating the kind of integration range involved
                inf = 1 corresponds to  (bound,+infinity),
                inf = -1            to  (-infinity,bound),
                inf = 2             to (-infinity,+infinity).

       epsabs - double precision
                absolute accuracy requested
       epsrel - double precision
                relative accuracy requested
                if  epsabs <= 0
                and epsrel < max(50*rel.mach.acc.,0.5d-28),
                the routine will end with ier = 6.

       limit  - int
                gives an upper bound on the number of subintervals
                in the partition of (a,b), limit >= 1

    on return
       result - double precision
                approximation to the integral

       abserr - double precision
                estimate of the modulus of the absolute error,
                which should equal or exceed abs(i-result)

       neval  - int
                number of integrand evaluations

       ier    - int
                ier = 0 normal and reliable termination of the
                        routine. it is assumed that the requested
                        accuracy has been achieved.
              - ier > 0 abnormal termination of the routine. the
                        estimates for result and error are less
                        reliable. it is assumed that the requested
                        accuracy has not been achieved.
       error messages
                ier = 1 maximum number of subdivisions allowed
                        has been achieved. one can allow more
                        subdivisions by increasing the value of
                        limit (and taking the according dimension
                        adjustments into account). however,if
                        this yields no improvement it is advised
                        to analyze the integrand in order to
                        determine the integration difficulties.
                        if the position of a local difficulty can
                        be determined (e.g. singularity,
                        discontinuity within the interval) one
                        will probably gain from splitting up the
                        interval at this point and calling the
                        integrator on the subranges. if possible,
                        an appropriate special-purpose integrator
                        should be used, which is designed for
                        handling the type of difficulty involved.
                    = 2 the occurrence of roundoff error is
                        detected, which prevents the requested
                        tolerance from being achieved.
                        the error may be under-estimated.
                    = 3 extremely bad integrand behaviour occurs
                        at some points of the integration
                        interval.
                    = 4 the algorithm does not converge.
                        roundoff error is detected in the
                        extrapolation table.
                        it is assumed that the requested tolerance
                        cannot be achieved, and that the returned
                        result is the best which can be obtained.
                    = 5 the integral is probably divergent, or
                        slowly convergent. it must be noted that
                        divergence can occur with any other value
                        of ier.
                    = 6 the input is invalid, because
                        (epsabs <= 0 and
                         epsrel < max(50*rel.mach.acc.,0.5d-28),
                        result, abserr, neval, last, rlist(1),
                        elist(1) and iord(1) are set to zero.
                        alist(1) and blist(1) are set to 0
                        and 1 respectively.

       alist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the left
                end points of the subintervals in the partition
                of the transformed integration range (0,1).

       blist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the right
                end points of the subintervals in the partition
                of the transformed integration range (0,1).

       rlist  - double precision
                vector of dimension at least limit, the first
                 last  elements of which are the integral
                approximations on the subintervals

       elist  - double precision
                vector of dimension at least limit,  the first
                last elements of which are the moduli of the
                absolute error estimates on the subintervals

       iord   - int
                vector of dimension limit, the first k
                elements of which are pointers to the
                error estimates over the subintervals,
                such that elist(iord(1)), ..., elist(iord(k))
                form a decreasing sequence, with k = last
                if last <= (limit/2+2), and k = limit+1-last
                otherwise

       last   - int
                number of subintervals actually produced
                in the subdivision process

routines called dqelg,dqk15i,dqpsrt end prologue dqagie

       the dimension of rlist2 is determined by the value of
       limexp in subroutine dqelg.

       list of major variables
       -----------------------

      alist     - list of left end points of all subintervals
                  considered up to now
      blist     - list of right end points of all subintervals
                  considered up to now
      rlist(i)  - approximation to the integral over
                  (alist(i),blist(i))
      rlist2    - array of dimension at least (limexp+2),
                  containing the part of the epsilon table
                  wich is still needed for further computations
      elist(i)  - error estimate applying to rlist(i)
      maxerr    - pointer to the interval with largest error
                  estimate
      errmax    - elist(maxerr)
      erlast    - error on the interval currently subdivided
                  (before that subdivision has taken place)
      area      - sum of the integrals over the subintervals
      errsum    - sum of the errors over the subintervals
      errbnd    - requested accuracy max(epsabs,epsrel*
                  abs(result))
           1    - variable for the left subinterval
           2    - variable for the right subinterval
      last      - index for subdivision
      nres      - number of calls to the extrapolation routine
      numrl2    - number of elements currently in rlist2. if an
                  appropriate approximation to the compounded
                  integral has been obtained, it is put in
                  rlist2(numrl2) after numrl2 has been increased
                  by one.
      small     - length of the smallest interval considered up
                  to now, multiplied by 1.5
      erlarg    - sum of the errors over the intervals larger
                  than the smallest interval considered up to now
      extrap    - logical variable denoting that the routine
                  is attempting to perform extrapolation. i.e.
                  before subdividing the smallest interval we
                  try to decrease the value of erlarg.
      noext     - logical variable denoting that extrapolation
                  is no longer allowed (true-value)

       machine dependent constants
       ---------------------------

      epmach is the largest relative spacing.
      uflow is the smallest positive magnitude.
      oflow is the largest positive magnitude.

Definition at line 255 of file integrate.hpp.

template<class Float , class integr_fn >

void gauss_kronrod::Rdqags	(	integr_fn	f,
		void *	ex,
		Float *	a,
		Float *	b,
		Float *	epsabs,
		Float *	epsrel,
		Float *	result,
		Float *	abserr,
		int *	neval,
		int *	ier,
		int *	limit,
		int *	lenw,
		int *	last,
		int *	iwork,
		Float *	work
	)

Definition at line 806 of file integrate.hpp.

Referenced by gauss_kronrod::Integral< gauss_kronrod::mvIntegral::evaluator >::operator()().

template<class Float , class integr_fn >

static void gauss_kronrod::rdqagse	(	integr_fn	f,
		void *	ex,
		Float *	a,
		Float *	b,
		Float *	epsabs,
		Float *	epsrel,
		int *	limit,
		Float *	result,
		Float *	abserr,
		int *	neval,
		int *	ier,
		Float *	alist,
		Float *	blist,
		Float *	rlist,
		Float *	elist,
		int *	iord,
		int *	last
	)

static

Definition at line 987 of file integrate.hpp.

template<class Float >

static void gauss_kronrod::rdqelg	(	int *	n,
		Float *	epstab,
		Float *	result,
		Float *	abserr,
		Float *	res3la,
		int *	nres
	)

static

Definition at line 1723 of file integrate.hpp.

template<class Float , class integr_fn >

static void gauss_kronrod::rdqk15i	(	integr_fn	f,
		void *	ex,
		Float *	boun,
		int *	inf,
		Float *	a,
		Float *	b,
		Float *	result,
		Float *	abserr,
		Float *	resabs,
		Float *	resasc
	)

static

Definition at line 1493 of file integrate.hpp.

template<class Float , class integr_fn >

static void gauss_kronrod::rdqk21	(	integr_fn	f,
		void *	ex,
		Float *	a,
		Float *	b,
		Float *	result,
		Float *	abserr,
		Float *	resabs,
		Float *	resasc
	)

static

Definition at line 1924 of file integrate.hpp.

template<class Float >

static void gauss_kronrod::rdqpsrt	(	int *	limit,
		int *	last,
		int *	maxerr,
		Float *	ermax,
		Float *	elist,
		int *	iord,
		int *	nrmax
	)

static

Definition at line 2133 of file integrate.hpp.

template<class T >

double gauss_kronrod::value ( T x )

Definition at line 27 of file integrate.hpp.

Referenced by fmax2(), and fmin2().

Classes

Functions

Function Documentation