NAG Library Routine Document
d01eaf (md_adapt_multi)
1
Purpose
d01eaf computes approximations to the integrals of a vector of similar functions, each defined over the same multidimensional hyperrectangular 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 nagmk26.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) 

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:
where
${f}_{i}={f}_{i}\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$, 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
$\Vert {\mathbf{absest}}\Vert \le \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{absreq}},\Vert {\mathbf{finest}}\Vert \times {\mathbf{relreq}}\right)$ where the norm
$\Vert .\Vert $ 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 Ndimensional 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 Ndimensional cube J. Comput. Appl. Math. 2 207–217
5
Arguments
 1: $\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integrals.
Constraint:
${\mathbf{ndim}}\ge 1$.
 2: $\mathbf{a}\left({\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the lower limits of integration,
${a}_{i}$, for $\mathit{i}=1,2,\dots ,n$.
 3: $\mathbf{b}\left({\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the upper limits of integration,
${b}_{i}$, for $\mathit{i}=1,2,\dots ,n$.
 4: $\mathbf{mincls}$ – IntegerInput/Output

On entry: must be set either to the minimum number of
funsub calls to be allowed, in which case
${\mathbf{mincls}}\ge 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 (
${\mathbf{mincls}}<0$ on entry) this is the number of new
funsub calls on the current call to
d01eaf.
 5: $\mathbf{maxcls}$ – IntegerInput

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:
 ${\mathbf{maxcls}}\ge {\mathbf{mincls}}$;
 ${\mathbf{maxcls}}\ge r$;
 $\text{where}r={2}^{n}+2{n}^{2}+2n+1\text{, if}n<11\text{, or}r=1+n\left(4{n}^{2}6n+14\right)/3\text{, if}n\ge 11$.
 6: $\mathbf{nfun}$ – IntegerInput

On entry: $m$, the number of integrands.
Constraint:
${\mathbf{nfun}}\ge 1$.
 7: $\mathbf{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
Integer, Intent (In)  ::  ndim, nfun  Real (Kind=nag_wp), Intent (In)  ::  z(ndim)  Real (Kind=nag_wp), Intent (Out)  ::  f(nfun) 

C Header Interface
#include nagmk26.h
void 
funsub (const Integer *ndim, const double z[], const Integer *nfun, double f[]) 

 1: $\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integrals.
 2: $\mathbf{z}\left({\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the coordinates of the point at which the integrands must be evaluated.
 3: $\mathbf{nfun}$ – IntegerInput

On entry: $m$, the number of integrands.
 4: $\mathbf{f}\left({\mathbf{nfun}}\right)$ – Real (Kind=nag_wp) arrayOutput

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 floatingpoint 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: $\mathbf{absreq}$ – Real (Kind=nag_wp)Input

On entry: the absolute accuracy required by you.
Constraint:
${\mathbf{absreq}}\ge 0.0$.
 9: $\mathbf{relreq}$ – Real (Kind=nag_wp)Input

On entry: the relative accuracy required by you.
Constraint:
${\mathbf{relreq}}\ge 0.0$.
 10: $\mathbf{lenwrk}$ – IntegerInput

On entry: the dimension of the array
wrkstr as declared in the (sub)program from which
d01eaf is called.
Suggested value:
${\mathbf{lenwrk}}\ge 6n+9m+\left(n+m+2\right)\left(1+p/r\right)$, 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:
${\mathbf{lenwrk}}\ge 8\times {\mathbf{ndim}}+11\times {\mathbf{nfun}}+3$.
 11: $\mathbf{wrkstr}\left({\mathbf{lenwrk}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{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: $\mathbf{finest}\left({\mathbf{nfun}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: ${\mathbf{finest}}\left(\mathit{i}\right)$ specifies the best estimate obtained from the $\mathit{i}$th integral, for $\mathit{i}=1,2,\dots ,m$.
 13: $\mathbf{absest}\left({\mathbf{nfun}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: ${\mathbf{absest}}\left(\mathit{i}\right)$ specifies the estimated absolute accuracy of ${\mathbf{finest}}\left(\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,m$.
 14: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if
${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit, the recommended value is
$1$.
When the value $\mathbf{1}\text{ or}1$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

maxcls was too small for
d01eaf to obtain the required accuracy. The arrays
finest and
absest respectively contain current estimates for the integrals and errors.
 ${\mathbf{ifail}}=2$

lenwrk is too small for the routine to continue. The arrays
finest and
absest respectively contain current estimates for the integrals and errors.
 ${\mathbf{ifail}}=3$

On a continuation call,
maxcls was set too small to make any progress. Increase
maxcls before calling
d01eaf again.
 ${\mathbf{ifail}}=4$

On entry,  ${\mathbf{ndim}}<1$, 
or  ${\mathbf{nfun}}<1$, 
or  ${\mathbf{maxcls}}<{\mathbf{mincls}}$, 
or  ${\mathbf{maxcls}}<r$ (see maxcls), 
or  ${\mathbf{absreq}}<0.0$, 
or  ${\mathbf{relreq}}<0.0$, 
or  ${\mathbf{lenwrk}}<8\times {\mathbf{ndim}}+11\times {\mathbf{nfun}}+3$. 
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
An absolute error estimate for each integrand is output in the array
absest. The routine exits with
${\mathbf{ifail}}={\mathbf{0}}$ if
8
Parallelism and Performance
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 implementationspecific information.
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 ${\mathbf{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
${\mathbf{mincls}}\ge 2r$, (see specification of
maxcls) and the results compared with those from the previous call.
If the routine is called with
${\mathbf{mincls}}\ge 2r$, 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+{\mathrm{log}}_{10}\left({n}^{3}\right)$ decimal digits may be inaccurate when using d01eaf for large values of $n$.
10
Example
This example computes
where
$j=1,2,\dots ,10$,
${f}_{j}=\mathrm{ln}\left({x}_{1}+2{x}_{2}+3{x}_{3}+4{x}_{4}\right)\mathrm{sin}\left(j+{x}_{1}+2{x}_{2}+3{x}_{3}+4{x}_{4}\right)$. The program is intended to show how to exploit the continuation facility provided with
d01eaf: the routine exits with
${\mathbf{ifail}}={\mathbf{1}}$ (printing an explanatory error message) and is reentered with
maxcls reset to a larger value. The program can be used with any values of
ndim and
nfun, except that the expression for
$r$ must be changed if
${\mathbf{ndim}}>10$ (see specification of
maxcls).
10.1
Program Text
Program Text (d01eafe.f90)
10.2
Program Data
None.
10.3
Program Results
Program Results (d01eafe.r)