## ▸▿ Contents

Settings help

FL Name Style:

FL Specification Language:

## 1Purpose

d01fcf attempts to evaluate a multidimensional integral (up to $15$ dimensions), with constant and finite limits, to a specified relative accuracy, using an adaptive subdivision strategy.

## 2Specification

Fortran Interface
 Subroutine d01fcf ( ndim, a, b, f, eps, acc,
 Integer, Intent (In) :: ndim, maxpts, lenwrk Integer, Intent (Inout) :: minpts, ifail Real (Kind=nag_wp), External :: f Real (Kind=nag_wp), Intent (In) :: a(ndim), b(ndim), eps Real (Kind=nag_wp), Intent (Out) :: acc, wrkstr(lenwrk), finval
#include <nag.h>
 void d01fcf_ (const Integer *ndim, const double a[], const double b[], Integer *minpts, const Integer *maxpts, double (NAG_CALL *f)(const Integer *ndim, const double z[]),const double *eps, double *acc, const Integer *lenwrk, double wrkstr[], double *finval, Integer *ifail)

## 3Description

d01fcf returns an estimate of a multidimensional integral over a hyper-rectangle (i.e., with constant limits), and also an estimate of the relative error. You set the relative accuracy required, return values for the integrand via a routine argument f, and also set the minimum and maximum acceptable number of calls to f (in minpts and maxpts).
The routine operates by repeated subdivision of the hyper-rectangular region into smaller hyper-rectangles. In each subregion, the integral is estimated using a seventh-degree rule, and an error estimate is obtained by comparison with a fifth-degree rule which uses a subset of the same points. The fourth differences of the integrand along each coordinate axis are evaluated, and the subregion is marked for possible future subdivision in half along that coordinate axis which has the largest absolute fourth difference.
If the estimated errors, totalled over the subregions, exceed the requested relative error (or if fewer than minpts calls to f have been made), further subdivision is necessary, and is performed on the subregion with the largest estimated error, that subregion being halved along the appropriate coordinate axis.
The routine will fail if the requested relative error level has not been attained by the time maxpts calls to f have been made; or, if the amount lenwrk of working storage is insufficient. A formula for the recommended value of lenwrk is given in Section 5. If a smaller value is used, and is exhausted in the course of execution, the routine switches to a less efficient mode of operation; only if this mode also breaks down is insufficient storage reported.
d01fcf is based on the HALF subroutine developed by van Dooren and de Ridder (1976). It uses a different basic rule, described in Genz and Malik (1980).

## 4References

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

## 5Arguments

1: $\mathbf{ndim}$Integer Input
On entry: $n$, the number of dimensions of the integral.
Constraint: $2\le {\mathbf{ndim}}\le 15$.
2: $\mathbf{a}\left({\mathbf{ndim}}\right)$Real (Kind=nag_wp) array Input
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) array Input
On entry: the upper limits of integration, ${b}_{i}$, for $\mathit{i}=1,2,\dots ,n$.
4: $\mathbf{minpts}$Integer Input/Output
On entry: must be set to the minimum number of integrand evaluations to be allowed.
On exit: contains the actual number of integrand evaluations used by d01fcf.
5: $\mathbf{maxpts}$Integer Input
On entry: the maximum number of integrand evaluations to be allowed.
Constraints:
• ${\mathbf{maxpts}}\ge {\mathbf{minpts}}$;
• ${\mathbf{maxpts}}\ge \alpha$, where $\alpha ={2}^{{\mathbf{ndim}}}+2×{{\mathbf{ndim}}}^{2}+2×{\mathbf{ndim}}+1$.
6: $\mathbf{f}$real (Kind=nag_wp) Function, supplied by the user. External Procedure
f must return the value of the integrand at a given point.
The specification of f is:
Fortran Interface
 Function f ( ndim, z)
 Real (Kind=nag_wp) :: f Integer, Intent (In) :: ndim Real (Kind=nag_wp), Intent (In) :: z(ndim)
 double f (const Integer *ndim, const double z[])
1: $\mathbf{ndim}$Integer Input
On entry: $n$, the number of dimensions of the integral.
2: $\mathbf{z}\left({\mathbf{ndim}}\right)$Real (Kind=nag_wp) array Input
On entry: the coordinates of the point at which the integrand $f$ must be evaluated.
f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d01fcf is called. Arguments denoted as Input must not be changed by this procedure.
Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01fcf. If your code inadvertently does return any NaNs or infinities, d01fcf is likely to produce unexpected results.
7: $\mathbf{eps}$Real (Kind=nag_wp) Input
On entry: the relative error acceptable to you. When the solution is zero or very small relative accuracy may not be achievable but you may still set eps to a reasonable value and check for the error exit ${\mathbf{ifail}}={\mathbf{2}}$.
Constraint: ${\mathbf{eps}}>0.0$.
8: $\mathbf{acc}$Real (Kind=nag_wp) Output
On exit: the estimated relative error in finval.
9: $\mathbf{lenwrk}$Integer Input
On entry: the dimension of the array wrkstr as declared in the (sub)program from which d01fcf is called.
Suggested value: for maximum efficiency, ${\mathbf{lenwrk}}\ge \left({\mathbf{ndim}}+2\right)×\left(1+{\mathbf{maxpts}}/\alpha \right)$ (see argument maxpts for $\alpha$).
If lenwrk is less than this, d01fcf will usually run less efficiently and may fail.
Constraint: ${\mathbf{lenwrk}}\ge 2×{\mathbf{ndim}}+4$.
10: $\mathbf{wrkstr}\left({\mathbf{lenwrk}}\right)$Real (Kind=nag_wp) array Workspace
11: $\mathbf{finval}$Real (Kind=nag_wp) Output
On exit: the best estimate obtained for the integral.
12: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $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$ or $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 ${\mathbf{ifail}}\ne {\mathbf{0}}$ on exit. When the value $-\mathbf{1}$ or $\mathbf{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).

## 6Error 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:
Note: in some cases d01fcf may return useful information.
${\mathbf{ifail}}=1$
On entry, ${\mathbf{eps}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{eps}}>0.0$.
On entry, lenwrk is too small. ${\mathbf{lenwrk}}=⟨\mathit{\text{value}}⟩$. Minimum possible dimension: $⟨\mathit{\text{value}}⟩$.
On entry, maxpts is too small. ${\mathbf{maxpts}}=⟨\mathit{\text{value}}⟩$. Minimum possible dimension: $⟨\mathit{\text{value}}⟩$.
On entry, ${\mathbf{minpts}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{maxpts}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{minpts}}\le {\mathbf{maxpts}}$.
On entry, ${\mathbf{ndim}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{ndim}}\le 15$.
On entry, ${\mathbf{ndim}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{ndim}}\ge 2$.
${\mathbf{ifail}}=2$
maxpts too small to obtain requested accuracy eps: ${\mathbf{maxpts}}=⟨\mathit{\text{value}}⟩$ and ${\mathbf{eps}}=⟨\mathit{\text{value}}⟩$.
${\mathbf{ifail}}=3$
lenwrk was too small to complete computation. finval and acc contain estimates of integral and relative error, but acc is greater than eps.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{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.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

A relative error estimate is output through the argument acc.

## 8Parallelism and Performance

d01fcf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d01fcf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
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.

Execution time will usually be dominated by the time taken to evaluate f, and hence the maximum time that could be taken will be proportional to maxpts.

## 10Example

This example estimates the integral
 $∫01 ∫01 ∫01 ∫01 4 z1 z32 exp(2z1z3) (1+z2+z4) 2 dz4 dz3 dz2 dz1 = 0.575364 .$
The accuracy requested is one part in $10000$.

### 10.1Program Text

Program Text (d01fcfe.f90)

None.

### 10.3Program Results

Program Results (d01fcfe.r)