d01eaf computes approximations to the integrals of a vector of similar functions, each defined over the same multidimensional hyper-rectangular region. The routine uses an adaptive subdivision strategy, and also computes absolute error estimates.

2 Specification

Fortran Interface

Subroutine d01eaf (

ndim, a, b, mincls, maxcls, nfun, funsub, absreq, relreq, lenwrk, wrkstr, finest, absest, ifail)

Integer, Intent (In)	::	ndim, maxcls, nfun, lenwrk
Integer, Intent (Inout)	::	mincls, ifail
Real (Kind=nag_wp), Intent (In)	::	a(ndim), b(ndim), absreq, relreq
Real (Kind=nag_wp), Intent (Inout)	::	wrkstr(lenwrk)
Real (Kind=nag_wp), Intent (Out)	::	finest(nfun), absest(nfun)
External	::	funsub

C Header Interface

#include <nag.h>

void

d01eaf_ (const Integer *ndim, const double a[], const double b[], Integer *mincls, const Integer *maxcls, const Integer *nfun,
void (NAG_CALL *funsub)(const Integer *ndim, const double z[], const Integer *nfun, double f[]),
const double *absreq, const double *relreq, const Integer *lenwrk, double wrkstr[], double finest[], double absest[], Integer *ifail)

The routine may be called by the names d01eaf or nagf_quad_md_adapt_multi.

3 Description

d01eaf uses a globally adaptive method based on the algorithm described by van Dooren and de Ridder (1976) and Genz and Malik (1980). It is implemented for integrals in the form:

\int_{a_{1}}^{b_{1}} \int_{a_{2}}^{b_{2}} \dots \int_{a_{n}}^{b_{n}} (f_{1}, f_{2}, \dots, f_{m}) d x_{n} \dots d x_{2} d x_{1},

where

f_{i} = f_{i} (x_{1}, x_{2}, \dots, x_{n})

, for

i = 1, 2, \dots, m

Upon entry, unless mincls has been set to a value less than or equal to

0

, d01eaf divides the integration region into a number of subregions with randomly selected volumes. Inside each subregion the integrals and their errors are estimated. The initial number of subregions is chosen to be as large as possible without using more than mincls calls to funsub. The results are stored in a partially ordered list (a heap). The routine then proceeds in stages. At each stage the subregion with the largest error (measured using the maximum norm) is halved along the coordinate axis where the integrands have largest absolute fourth differences. The basic rule is applied to each half of this subregion and the results are stored in the list. The results from the two halves are used to update the global integral and error estimates (finest and absest) and the routine continues unless

‖ absest ‖ \leq \max (absreq, ‖ finest ‖ \times relreq)

where the norm

‖ . ‖

is the maximum norm, or further subdivision would use more than maxcls calls to funsub. If at some stage there is insufficient working storage to keep the results for the next subdivision, the routine switches to a less efficient mode; only if this mode of operation breaks down is insufficient storage reported.

4 References

Genz A C and Malik A A (1980) An adaptive algorithm for numerical integration over an N-dimensional rectangular region J. Comput. Appl. Math. 6 295–302

van Dooren P and de Ridder L (1976) An adaptive algorithm for numerical integration over an N-dimensional cube J. Comput. Appl. Math. 2 207–217

5 Arguments

1: $ndim$ – Integer Input

On entry:

n

, the number of dimensions of the integrals.

Constraint:

ndim \geq 1

2: $a (ndim)$ – Real (Kind=nag_wp) array Input

On entry: the lower limits of integration,

a_{i}

, for

i = 1, 2, \dots, n

3: $b (ndim)$ – Real (Kind=nag_wp) array Input

On entry: the upper limits of integration,

b_{i}

, for

i = 1, 2, \dots, n

4: $mincls$ – Integer Input/Output

On entry: must be set either to the minimum number of funsub calls to be allowed, in which case

mincls \geq 0

or to a negative value. In this case, the routine continues the calculation started in a previous call with the same integrands and integration limits: no arguments other than mincls, maxcls, absreq, relreq or ifail must be changed between the calls.

On exit: gives the number of funsub calls actually used by d01eaf. For the continuation case (

mincls < 0

on entry) this is the number of new funsub calls on the current call to d01eaf.

5: $maxcls$ – Integer Input

On entry: the maximum number of funsub calls to be allowed. In the continuation case this is the number of new funsub calls to be allowed.

Constraints:

$maxcls \geq mincls$ ;
$maxcls \geq r$ ;
$where r = 2^{n} + 2 n^{2} + 2 n + 1, if n < 11, or r = 1 + n (4 n^{2} - 6 n + 14) / 3, if n \geq 11$ .

6: $nfun$ – Integer Input

On entry:

m

, the number of integrands.

Constraint:

nfun \geq 1

7: $funsub$ – Subroutine, supplied by the user. External Procedure

funsub must evaluate the integrands

f_{i}

at a given point.

The specification of funsub is:

Fortran Interface

Subroutine funsub (

ndim, z, nfun, f)

Integer, Intent (In)	::	ndim, nfun
Real (Kind=nag_wp), Intent (In)	::	z(ndim)
Real (Kind=nag_wp), Intent (Out)	::	f(nfun)

C Header Interface

void	funsub (const Integer ndim, const double z[], const Integer nfun, double f[])

1: $ndim$ – Integer Input: On entry: $n$ , the number of dimensions of the integrals.
2: $z (ndim)$ – Real (Kind=nag_wp) array Input: On entry: the coordinates of the point at which the integrands must be evaluated.
3: $nfun$ – Integer Input: On entry: $m$ , the number of integrands.
4: $f (nfun)$ – Real (Kind=nag_wp) array Output: On exit: the value of the $i$ th integrand at the given point.

funsub must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01eaf is called. Arguments denoted as Input must not be changed by this procedure.

Note: funsub should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01eaf. If your code inadvertently does return any NaNs or infinities, d01eaf is likely to produce unexpected results.

8: $absreq$ – Real (Kind=nag_wp) Input

On entry: the absolute accuracy required by you.

Constraint:

absreq \geq 0.0

9: $relreq$ – Real (Kind=nag_wp) Input

On entry: the relative accuracy required by you.

Constraint:

relreq \geq 0.0

10: $lenwrk$ – Integer Input

On entry: the dimension of the array wrkstr as declared in the (sub)program from which d01eaf is called.

Suggested value:

lenwrk \geq 6 n + 9 m + (n + m + 2) (1 + p / r)

, where

p

is the value of maxcls and

r

is defined under maxcls. If lenwrk is significantly smaller than this, the routine will not work as efficiently and may even fail.

Constraint:

lenwrk \geq 8 \times ndim + 11 \times nfun + 3

11: $wrkstr (lenwrk)$ – Real (Kind=nag_wp) array Input/Output

On entry: if

mincls < 0

, wrkstr must be unchanged from the previous call of d01eaf.

On exit: contains information about the current subdivision which could be used in a continuation call.

12: $finest (nfun)$ – Real (Kind=nag_wp) array Output

On exit:

finest (i)

specifies the best estimate obtained from the

i

th integral, for

i = 1, 2, \dots, m

13: $absest (nfun)$ – Real (Kind=nag_wp) array Output

On exit:

absest (i)

specifies the estimated absolute accuracy of

finest (i)

, for

i = 1, 2, \dots, m

14: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

−1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

Note: in some cases d01eaf may return useful information.

$ifail = 1$: maxcls too small to obtain required accuracy. $maxcls = ⟨ value ⟩$ .

$ifail = 2$: lenwrk too small for the routine to continue. $lenwrk = ⟨ value ⟩$ .

$ifail = 3$: maxcls too small to make any progress. $maxcls = ⟨ value ⟩$ .

$ifail = 4$: On entry, lenwrk is too small. $lenwrk = ⟨ value ⟩$ . Minimum possible dimension: $⟨ value ⟩$ .

On entry, $maxcls = ⟨ value ⟩$ and $mincls = ⟨ value ⟩$ .
Constraint: $maxcls \geq mincls$ .

On entry, $maxcls = ⟨ value ⟩$ and $R = ⟨ value ⟩$ .
Constraint: $maxcls \geq R = 1 + ndim \times (4 \times {ndim}^{2} - 6 \times ndim + 14) / 3$ .

On entry, $maxcls = ⟨ value ⟩$ and $R = ⟨ value ⟩$ .
Constraint: $maxcls \geq R = 2^{ndim} + 2 \times {ndim}^{2} + 2 \times ndim + 1$ .

On entry, $ndim = ⟨ value ⟩$ .
Constraint: $ndim \geq 1$ .

On entry, $nfun = ⟨ value ⟩$ .
Constraint: $nfun \geq 1$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

An absolute error estimate for each integrand is output in the array absest. The routine exits with

ifail = 0

\max_{i} (absest (i)) \leq \max (absreq, relreq \times \max_{i} | finest (i) |) .

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

d01eaf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

Usually the running time for d01eaf will be dominated by the time in funsub, so the maximum time that could be used by d01eaf will be proportional to maxcls multiplied by the cost of a call to funsub.

On a normal call, you should set

mincls = 0

on entry.

For some integrands, particularly those that are poorly behaved in a small part of the integration region, d01eaf may terminate prematurely with values of absest that are significantly smaller than the actual absolute errors. This behaviour should be suspected if the returned value of mincls is small relative to the expected difficulty of the integrals. When this occurs d01eaf should be called again, but with an entry value of

mincls \geq 2 r

, (see specification of maxcls) and the results compared with those from the previous call.

If the routine is called with

mincls \geq 2 r

, the exact values of finest and absest on return will depend (within statistical limits) on the sequence of random numbers generated internally within d01eaf by calls to g05saf. Separate runs will produce identical answers unless the part of the program executed prior to calling d01eaf also calls (directly or indirectly) routines from Chapter G05, and, in addition, the series of such calls differs between runs.

Because of moderate instability in the application of the basic integration rule, approximately the last

1 + \log_{10} (n^{3})

decimal digits may be inaccurate when using d01eaf for large values of

n

10 Example

This example computes

\int_{0}^{1} \int_{0}^{1} \int_{0}^{1} \int_{0}^{1} (f_{1}, f_{2}, \dots, f_{10}) d x_{4} d x_{3} d x_{2} d x_{1},

where

j = 1, 2, \dots, 10

f_{j} = \ln (x_{1} + 2 x_{2} + 3 x_{3} + 4 x_{4}) \sin (j + x_{1} + 2 x_{2} + 3 x_{3} + 4 x_{4})

. The program is intended to show how to exploit the continuation facility provided with d01eaf: the routine exits with

ifail = 1

(printing an explanatory error message) and is re-entered with maxcls reset to a larger value. The program can be used with any values of ndim and nfun, except that the expression for