NAG CL Interface
d01rac (dim1_gen_vec_multi_rcomm)
Note: this function uses optional parameters to define choices in the problem specification and in the details of the algorithm. If you wish to use default
settings for all of the optional parameters, you need only read Sections 1 to 10 of this document. If, however, you wish to reset some or all of the settings please refer to Section 11 for a detailed description of the specification of the optional parameters.
1
Purpose
d01rac is a general purpose adaptive integrator which calculates an approximation to a vector of definite integrals
$\mathbf{F}$ over a finite range
$\left[a,b\right]$, given the vector of integrands
$\mathbf{f}\left(x\right)$.
If the same subdivisions of the range are equally good for functions ${f}_{1}\left(x\right)$ and ${f}_{2}\left(x\right)$, because ${f}_{1}\left(x\right)$ and ${f}_{2}\left(x\right)$ have common areas of the range where they vary slowly and where they vary quickly, then we say that ${f}_{1}\left(x\right)$ and ${f}_{2}\left(x\right)$ are ‘similar’. d01rac is particularly effective for the integration of a vector of similar functions.
2
Specification
void 
d01rac (Integer *irevcm,
Integer ni,
double a,
double b,
Integer *sid,
Integer needi[],
double x[],
Integer lenx,
Integer *nx,
const double fm[],
Integer ldfm,
double dinest[],
double errest[],
const Integer iopts[],
const double opts[],
Integer icom[],
Integer licom,
double com[],
Integer lcom,
NagError *fail) 

The function may be called by the names: d01rac, nag_quad_dim1_gen_vec_multi_rcomm or nag_quad_1d_gen_vec_multi_rcomm.
3
Description
d01rac is an extension to various QUADPACK routines, including QAG, QAGS and QAGP. The extensions made allow multiple integrands to be evaluated simultaneously, using a vectorized interface and reverse communication.
The quadrature scheme employed by
d01rac can be chosen by you. Six Gauss–Kronrod schemes are available. The algorithm incorporates a global acceptance criterion (as defined by
Malcolm and Simpson (1976)), optionally together with the εalgorithm (see
Wynn (1956)) to perform extrapolation. The local error estimation is described in
Piessens et al. (1983).
d01rac is the integration function in the suite of functions
d01rac and
d01rcc. It also uses optional parameters, which can be set and queried using the functions
d01zkc and
d01zlc respectively. The options available are described in
Section 11.
First, the option arrays
iopts and
opts must be initialized using
d01zkc. Thereafter any required options must be set before calling
d01rac, or the function
d01rcc.
A typical usage of this suite of functions is (in pseudocode for clarity),
Setup phase
liopts = 100
lopts = 100
allocate(iopts(liopts),opts(lopts))
d01zkc('initialize = d01rac',iopts,liopts,opts,lopts,fail)
d01zkc('option = value',iopts,liopts,opts,lopts,fail)
...
d01rcc(ni,lenxrq,ldfmrq,sdfmrq,licmin,licmax,lcmin,lcmax,
iopts,opts,fail)
lenx = lenxrq
ldfm = ldfmrq
sdfm = sdfmrq
licom = licmax
lcom = lcmax
allocate(icom(licom),com(lcom),x(lenx),fm(ldfm,sdfm),needi(ni),
dinest(ni),errest(ni))
Solve phase
irevcm = 1
while irevcm≠0
d01rac(irevcm,ni,a,b,sid,needi,x,lenx,nx,fm,ldfm,
dinest,errest,iopts,opts,icom,licom,com,
lcom,fail)
select case(irevcm)
case(11)
Initial solve phase
evaluate fm(1:ni,1:nx)
case(12)
Adaptive solve phase
evaluate fm(needi(1:ni)=1,1:nx)
case(0)
investigate fail
end select
end while
Diagnostic phase
d01zlc('option',ivalue,rvalue,cvalue,optype,iopts,opts,fail)
...
During the initial solve phase, the first estimation of the definite integral and error estimate is constructed over the interval $\left[a,b\right]$. This will have been divided into ${s}_{\mathit{pri}}$ level $1$ segments, where ${s}_{\mathit{pri}}$ is the number of ${\mathbf{Primary\; Divisions}}$, and will use any provided breakpoints if ${\mathbf{Primary\; Division\; Mode}}=\mathrm{MANUAL}$.
Once a complete integral estimate over $\left[a,b\right]$ is available, i.e., after all the estimates for the level $1$ segments have been evaluated, the function enters the adaptive phase. The estimated errors are tested against the requested tolerances ${\epsilon}_{a}$ and ${\epsilon}_{r}$, corresponding to the ${\mathbf{Absolute\; Tolerance}}$ and ${\mathbf{Relative\; Tolerance}}$ respectively. Should this test fail, and additional subdivision be allowed, a segment is selected for subdivision, and is subsequently replaced by two new segments at the next level of refinement. How this segment is chosen may be altered by setting ${\mathbf{Prioritize\; Error}}$ to either favour the segment with the maximum error, or the segment with the lowest level supporting an unacceptable (although potentially nonmaximal) error. Up to ${\mathrm{max}}_{\mathit{sdiv}}$ subdivisions are allowed if sufficient memory is provided, where ${\mathrm{max}}_{\mathit{sdiv}}$ is the value of ${\mathbf{Maximum\; Subdivisions}}$.
Once a sufficient number of error estimates have been constructed for a particular integral, the function may optionally use ${\mathbf{Extrapolation}}$ to attempt to accelerate convergence. This may significantly lower the amount of work required for a given integration. To minimize the risk of premature convergence from extrapolation, a safeguard ${\epsilon}_{\mathit{safe}}$ can be set using ${\mathbf{Extrapolation\; Safeguard}}$, and the extrapolated solution will only be considered if ${\epsilon}_{\mathit{safe}}{\epsilon}_{q}\le {\epsilon}_{\mathit{ex}}$, where ${\epsilon}_{q}$ and ${\epsilon}_{\mathit{ex}}$ are the estimated error directly from the quadrature and from the extrapolation respectively. If extrapolation is successful for the computation of integral $j$, the extrapolated solution will be returned in ${\mathbf{dinest}}\left[j1\right]$ on completion of d01rac. Otherwise the direct solution will be returned in ${\mathbf{dinest}}\left[j1\right]$. This is indicated by the value of ${\mathbf{needi}}\left[j1\right]$ on completion.
4
References
Malcolm M A and Simpson R B (1976) Local versus global strategies for adaptive quadrature ACM Trans. Math. Software 1 129–146
Piessens R (1973) An algorithm for automatic integration Angew. Inf. 15 399–401
Piessens R, de Doncker–Kapenga E, Überhuber C and Kahaner D (1983) QUADPACK, A Subroutine Package for Automatic Integration Springer–Verlag
Wynn P (1956) On a device for computing the ${e}_{m}\left({S}_{n}\right)$ transformation Math. Tables Aids Comput. 10 91–96
5
Arguments
Note: this function uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the argument
irevcm. Between intermediate exits and reentries,
all arguments other than irevcm, needi and fm must remain unchanged.
where ${\mathbf{FM}}\left(j,i\right)$ appears in this document it refers to the array element ${\mathbf{fm}}\left[\left(i1\right)\times {\mathbf{ldfm}}+j1\right]$.

1:
$\mathbf{irevcm}$ – Integer *
Input/Output

On initial entry:
${\mathbf{irevcm}}=1$.
 ${\mathbf{irevcm}}=1$
 Sets up data structures in icom and com and starts a new integration.
Constraint:
${\mathbf{irevcm}}=1$ on initial entry.
On intermediate exit:
${\mathbf{irevcm}}=11$ or
$12$.
irevcm requests the integrands
${f}_{j}\left({x}_{i}\right)$ be evaluated for all required
$j\in 1,\dots ,{n}_{i}$ as indicated by
needi, and at all the points
${x}_{\mathit{i}}$, for
$\mathit{i}=1,2,\dots ,{n}_{x}$. Abscissae
${x}_{i}$ are provided in
${\mathbf{x}}\left[i1\right]$ and
${f}_{j}\left({x}_{i}\right)$ must be returned in
${\mathbf{FM}}\left(j,i\right)$.
During the initial solve phase:
 ${\mathbf{irevcm}}=11$
 Function values are required to construct the initial estimates of the definite integrals.
If ${\mathbf{needi}}\left[j1\right]=1$, ${f}_{j}\left({x}_{i}\right)$ must be supplied in ${\mathbf{FM}}\left(j,i\right)$. This will be the case unless you have abandoned the evaluation of specific integrals on a previous call.
If ${\mathbf{needi}}\left[j1\right]=0$, you have previously abandoned the evaluation of integral $j$, and hence should not supply the value of ${f}_{j}\left({x}_{i}\right)$.
dinest and
errest contain incomplete information during this phase. As such you should not abandon the evaluation of any integrals during this phase unless you do not require their estimate.
If
irevcm is set to a negative value during this phase,
${\mathbf{needi}}\left[\mathit{j}1\right]$, for
$\mathit{j}=1,2,\dots ,{n}_{i}$, will be set to this negative value and
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_USER_STOP will be returned.
During the adaptive solve phase:
 ${\mathbf{irevcm}}=12$
 Function values are required to improve the estimates of the definite integrals.
If ${\mathbf{needi}}\left[j1\right]=0$, any evaluation of ${f}_{j}\left({x}_{i}\right)$ will be discarded, so there is no need to provide them.
If ${\mathbf{needi}}\left[j1\right]=1$, ${f}_{j}\left({x}_{i}\right)$ must be provided in ${\mathbf{FM}}\left(j,i\right)$.
If ${\mathbf{needi}}\left[j1\right]=2$, $3$ or $4$, the current error estimate of integral $j$ does not require integrand $j$ to be evaluated and provided in ${\mathbf{FM}}\left(j,i\right)$. Should you choose to, integrand $j$ can be evaluated in which case ${\mathbf{needi}}\left[j1\right]$ must be set to $1$.
dinest and
errest contain complete information during this phase.
If
irevcm is set to a negative value during this phase
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_ACCURACY or
NE_QUAD_BAD_SUBDIV_INT will be returned and the elements of
needi will reflect the current state of the adaptive process.
On intermediate reentry:
irevcm should normally be left unchanged. However, if
irevcm is set to a negative value,
d01rac will terminate, (see
${\mathbf{irevcm}}=11$ and
${\mathbf{irevcm}}=12$ above).
On final exit:
${\mathbf{irevcm}}=0$.
 ${\mathbf{irevcm}}=0$
 Indicates the algorithm has completed.
Note: any values you return to d01rac as part of the reverse communication procedure should not include floatingpoint NaN (Not a Number) or infinity values, since these are not handled by d01rac. If your code inadvertently does return any NaNs or infinities, d01rac is likely to produce unexpected results.

2:
$\mathbf{ni}$ – Integer
Input

On entry: ${n}_{i}$, the number of integrands.

3:
$\mathbf{a}$ – double
Input

On entry: $a$, the lower bound of the domain.

4:
$\mathbf{b}$ – double
Input

On entry:
$b$, the upper bound of the domain.
If
$\leftba\right<10\epsilon $, where
$\epsilon $ is the
machine precision (see
X02AJC),
d01rac will return
${\mathbf{dinest}}\left[\mathit{j}1\right]={\mathbf{errest}}\left[\mathit{j}1\right]=0.0$, for
$\mathit{j}=1,2,\dots ,{n}_{i}$.

5:
$\mathbf{sid}$ – Integer *
Output

For advanced users.
On intermediate exit:
sid identifies a specific set of abscissae,
$\mathbf{x}$, returned during the integration process. When a new set of abscissae are generated the value of
sid is incremented by
$1$. Advanced users may store calculations required for an identified set
$\mathbf{x}$, and reuse them should
d01rac return the same value of
sid, i.e., the same set of abscissae was used.

6:
$\mathbf{needi}\left[{\mathbf{ni}}\right]$ – Integer
Input/Output

On initial entry: need not be set.
On intermediate exit:
${\mathbf{needi}}\left[j1\right]$ indicates what action must be taken for integral
$j=1,2,\dots {n}_{i}$ (see
irevcm).
 ${\mathbf{needi}}\left[j1\right]=0$
 Do not provide ${f}_{j}\left({x}_{i}\right)$. Any provided values will be ignored.
 ${\mathbf{needi}}\left[j1\right]=1$
 The values
${f}_{j}\left({x}_{\mathit{i}}\right)$ must be provided in ${\mathbf{FM}}\left(j,\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{x}$.
 ${\mathbf{needi}}\left[j1\right]=2$
 The values ${f}_{j}\left({x}_{i}\right)$ are not required, however the error estimate for integral $j$ is still above the requested tolerance. If you wish to provide values for the evaluation of integral $j$, set ${\mathbf{needi}}\left[j1\right]=1$, and supply
${f}_{j}\left({x}_{\mathit{i}}\right)$ in ${\mathbf{FM}}\left(j,\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{x}$.
 ${\mathbf{needi}}\left[j1\right]=3$
 The error estimate for integral $j$ cannot be improved to below the requested tolerance directly, either because no more new splits may be performed due to exhaustion, or due to the detection of extremely bad integrand behaviour. However, providing the values ${f}_{j}\left({x}_{i}\right)$ may still lead to some improvement, and may lead to an acceptable error estimate indirectly using Wynn's epsilon algorithm. If you wish to provide values for the evaluation of integral $j$, set ${\mathbf{needi}}\left[j1\right]=1$, and supply
${f}_{j}\left({x}_{\mathit{i}}\right)$ in ${\mathbf{FM}}\left(j,\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{x}$.
 ${\mathbf{needi}}\left[j1\right]=4$
 The error estimate of integral $j$ is below the requested tolerance. If you believe this to be false, if for example the result in ${\mathbf{dinest}}\left[j1\right]$ is greatly different to what you may expect, you may force the algorithm to reevaluate this conclusion by including the evaluations of integrand $j$ at
${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{n}_{x}$, and setting ${\mathbf{needi}}\left[j1\right]=1$. Integral and error estimation will be performed again during the next iteration.
On intermediate reentry:
${\mathbf{needi}}\left[j1\right]$ may be used to indicate what action you have taken for integral
$j$.
 ${\mathbf{needi}}\left[j1\right]=1$
 You have provided the values
${f}_{j}\left({x}_{\mathit{i}}\right)$ in ${\mathbf{FM}}\left(j,\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{x}$.
 ${\mathbf{needi}}\left[j1\right]<0$
 You are abandoning the evaluation of integral $j$. The current values of ${\mathbf{dinest}}\left[j1\right]$ and ${\mathbf{errest}}\left[j1\right]$ will be returned on final completion.
Otherwise you have not provided the value ${f}_{j}\left({x}_{i}\right)$.
On final exit:
${\mathbf{needi}}\left[j1\right]$ indicates the final state of integral
$j$.
 ${\mathbf{needi}}\left[j1\right]=0$
 The error estimate for ${F}_{j}$ is below the requested tolerance.
 ${\mathbf{needi}}\left[j1\right]=1$
 The error estimate for ${F}_{j}$ is below the requested tolerance after extrapolation.
 ${\mathbf{needi}}\left[j1\right]=2$
 The error estimate for ${F}_{j}$ is above the requested tolerance.
 ${\mathbf{needi}}\left[j1\right]=3$
 The error estimate for ${F}_{j}$ is above the requested tolerance, and extremely bad behaviour of integral $j$ has been detected.
 ${\mathbf{needi}}\left[j1\right]<0$
 You prohibited further evaluation of integral $j$.

7:
$\mathbf{x}\left[{\mathbf{lenx}}\right]$ – double
Input/Output

On initial entry: if
${\mathbf{Primary\; Division\; Mode}}=\mathrm{AUTOMATIC}$,
x need not be set. This is the default behaviour.
If
${\mathbf{Primary\; Division\; Mode}}=\mathrm{MANUAL}$,
x is used to supply a set of initial ‘breakpoints’ inside the domain of integration. Specifically,
${\mathbf{x}}\left[i1\right]$ must contain a breakpoint
${x}_{\mathit{i}}^{0}$, for
$\mathit{i}=1,2,\dots ,\left({s}_{\mathit{pri}}1\right)$, where
${s}_{\mathit{pri}}$ is the number of
${\mathbf{Primary\; Divisions}}$.
Constraint:
if breakpoints are supplied, ${x}_{\mathit{i}}^{0}\in \left(a,b\right)$, $\left{x}_{\mathit{i}}^{0}a\right>10.0\epsilon $, $\left{x}_{\mathit{i}}^{0}b\right>10.0\epsilon $, for $\mathit{i}=1,2,\dots ,\left({s}_{\mathit{pri}}1\right)$.
On intermediate exit:
${\mathbf{x}}\left[\mathit{i}1\right]$ is the abscissa ${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{n}_{x}$, at which the appropriate integrals must be evaluated.

8:
$\mathbf{lenx}$ – Integer
Input

On entry: the dimension of the array
x. Currently
${\mathbf{lenx}}=\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(122,{s}_{\mathit{pri}}1\right)$ will be sufficient for all cases.
Constraint:
${\mathbf{lenx}}\ge \mathit{lenxrq}$, where
$\mathit{lenxrq}$ is dependent upon the options currently set (see
Section 11).
$\mathit{lenxrq}$ is returned as
lenxrq from
d01rcc.

9:
$\mathbf{nx}$ – Integer *
Input/Output

On initial entry: need not be set.
On intermediate exit:
${n}_{x}$, the number of abscissae at which integrands are required.
On intermediate reentry: must not be changed.

10:
$\mathbf{fm}\left[\mathit{dim}\right]$ – const double
Input

Note: the dimension,
dim, of the array
fm
must be at least
${\mathbf{ldfm}}\times \mathit{sdfmrq}$, where
$\mathit{sdfmrq}$ is dependent upon
${n}_{i}$ and the options currently set.
$\mathit{sdfmrq}$ is returned as
sdfmrq from
d01rcc. If default options are chosen,
$\mathit{sdfmrq}=\mathit{lenxrq}$.
On initial entry: need not be set.
On intermediate reentry: if indicated by ${\mathbf{needi}}\left[j1\right]$ you must supply the values ${f}_{j}\left({x}_{i}\right)$ in ${\mathbf{FM}}\left(\mathit{j},\mathit{i}\right)$, for $\mathit{i}=1,2,\dots ,{n}_{x}$ and $\mathit{j}=1,2,\dots ,{n}_{i}$.

11:
$\mathbf{ldfm}$ – Integer
Input

On entry: the stride separating matrix row elements in the array
fm.
Constraint:
${\mathbf{ldfm}}\ge \mathit{ldfmrq}$, where
$\mathit{ldfmrq}$ is dependent upon
${n}_{i}$ and the options currently set.
$\mathit{ldfmrq}$ is returned as
ldfmrq from
d01rcc. If default options are chosen,
$\mathit{ldfmrq}={n}_{i}$, implying
${\mathbf{ldfm}}\ge {\mathbf{ni}}$.

12:
$\mathbf{dinest}\left[{\mathbf{ni}}\right]$ – double
Input/Output

${\mathbf{dinest}}\left[j1\right]$ contains the current estimate of the definite integral ${F}_{j}$.
On initial entry: need not be set.
On intermediate reentry: must not be altered.
On exit: contains the current estimates of the
ni integrals. If
${\mathbf{irevcm}}=0$, this will be the final solution.

13:
$\mathbf{errest}\left[{\mathbf{ni}}\right]$ – double
Input/Output

${\mathbf{errest}}\left[j1\right]$ contains the current error estimate of the definite integral ${F}_{j}$.
On initial entry: need not be set.
On intermediate reentry: must not be altered.
On exit: contains the current error estimates for the
ni integrals. If
${\mathbf{irevcm}}=0$,
errest contains the final error estimates of the
ni integrals.

14:
$\mathbf{iopts}\left[100\right]$ – const Integer
Communication Array

15:
$\mathbf{opts}\left[100\right]$ – const double
Communication Array

The arrays
iopts and
opts MUST NOT be altered between calls to any of the functions
d01rac,
d01rcc,
d01zkc and
d01zlc.

16:
$\mathbf{icom}\left[{\mathbf{licom}}\right]$ – Integer
Communication Array

icom contains details of the integration procedure, including information on the integration of the
${n}_{i}$ integrals over individual segments. This data is stored sequentially in the order that segments are created. For further information see
Section 9.1.

17:
$\mathbf{licom}$ – Integer
Input

On entry: the dimension of the array
icom.
Constraint:
${\mathbf{licom}}\ge \mathit{licmin}$, where
$\mathit{licmin}$ is dependent upon
ni and the current options set.
$\mathit{licmin}$ is returned as
licmin from
d01rcc. If the default options are set,
$\mathit{licmin}=55+6\times {\mathbf{ni}}$. Larger values than
$\mathit{licmin}$ are recommended if you anticipate that any integrals will require the domain to be further subdivided.
The maximum value that may be required,
$\mathit{licmax}$, is returned as
licmax from
d01rcc. If default options are chosen, except for possibly increasing the value of
${s}_{\mathit{pri}}$,
$\mathit{licmax}=50+5\times {\mathbf{ni}}+\left({s}_{\mathit{pri}}+100\right)\times \left(5+{\mathbf{ni}}\right)$.

18:
$\mathbf{com}\left[{\mathbf{lcom}}\right]$ – double
Communication Array

com contains details of the integration procedure, including information on the integration of the
${n}_{i}$ integrals over individual segments. This data is stored sequentially in the order that segments are created. For further information see
Section 9.1.

19:
$\mathbf{lcom}$ – Integer
Input

On entry: the dimension of the array
com.
Constraint:
${\mathbf{lcom}}>\mathit{lcmin}$, where
$\mathit{lcmin}$ is dependent upon
ni,
${s}_{\mathit{pri}}$ and the current options set.
$\mathit{lcmin}$ is returned as
lcmin from
d01rcc. If default options are set,
$\mathit{lcmin}=96+12\times {\mathbf{ni}}$. Larger values are recommended if you anticipate that any integrals will require the domain to be further subdivided.
Given the current options and arguments, the maximum value,
$\mathit{lcmax}$, of
lcom that may be required, is returned as
lcmax from
d01rcc. If default options are chosen,
$\mathit{lcmax}=94+9\times {\mathbf{ni}}+\u2308{\mathbf{ni}}/2\u2309+\left({s}_{\mathit{pri}}+100\right)\times \left(2+\u2308{\mathbf{ni}}/2\u2309+2\times {\mathbf{ni}}\right)$.

20:
$\mathbf{fail}$ – NagError *
Input/Output

The NAG error argument (see
Section 7 in the Introduction to the NAG Library CL Interface).
6
Error Indicators and Warnings
 NE_ACCURACY

At least one error estimate exceeded the requested tolerances.
 NE_ALLOC_FAIL

Dynamic memory allocation failed.
See
Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
 NE_BAD_PARAM

On entry, argument $\u2329\mathit{\text{value}}\u232a$ had an illegal value.
 NE_INT

irevcm had an illegal value.
On entry,
${\mathbf{irevcm}}=\u2329\mathit{\text{value}}\u232a$.
On entry, ${\mathbf{ni}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ni}}\ge 1$.
 NE_INT_2

${\mathbf{ldfm}}<\mathit{ldfmrq}$. If default options are chosen, this implies ${\mathbf{ldfm}}<{\mathbf{ni}}$.
On entry, ${\mathbf{ldfm}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{ldfm}}\ge \u2329\mathit{\text{value}}\u232a$.
 NE_INTERNAL_ERROR

An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact
NAG for assistance.
See
Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
 NE_INVALID_ARRAY

On entry, one of
icom and
com has become corrupted.
 NE_INVALID_OPTION

Either the option arrays
iopts and
opts have not been initialized for
d01rac, or they have become corrupted.
 NE_NO_LICENCE

Your licence key may have expired or may not have been installed correctly.
See
Section 8 in the Introduction to the NAG Library CL Interface for further information.
 NE_QUAD_BAD_SUBDIV_INT

Extremely bad behaviour was detected for at least one integral.
Extremely bad behaviour was detected for at least one integral. At least one other integral error estimate was above the requested tolerance.
 NE_QUAD_BRKPTS_INVAL

On entry,
${\mathbf{Primary\; Division\; Mode}}=\mathrm{MANUAL}$ and at least one supplied breakpoint in
x is outside of the domain of integration.
 NE_TOO_SMALL

lcom is insufficient for additional subdivision.
On entry,
${\mathbf{lcom}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{lcom}}\ge \u2329\mathit{\text{value}}\u232a$.
lenx is insufficient for the chosen options.
On entry,
${\mathbf{lenx}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{lenx}}\ge \u2329\mathit{\text{value}}\u232a$.
licom is insufficient for additional subdivision.
On entry,
${\mathbf{licom}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{licom}}\ge \u2329\mathit{\text{value}}\u232a$.
 NE_USER_STOP

Evaluation of all integrals has been stopped during the initial phase.
7
Accuracy
d01rac cannot guarantee, but in practice usually achieves, the following accuracy for each integral
${F}_{j}$:
where
${\epsilon}_{a}$ and
${\epsilon}_{r}$ are the error tolerances
${\mathbf{Absolute\; Tolerance}}$ and
${\mathbf{Relative\; Tolerance}}$ respectively. Moreover, it returns
errest, the entries of which in normal circumstances satisfy,
8
Parallelism and Performance
d01rac is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
d01rac 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 function. Please also consult the
Users' Note for your implementation for any additional implementationspecific information.
The time required by d01rac is usually dominated by the time required to evaluate the values of the integrands ${f}_{j}$.
d01rac will be most efficient if any badly behaved integrands provided have irregularities over similar subsections of the domain. For example, evaluation of the integrals,
will be quite efficient, as the irregular behaviour of the first two integrands is at
$x=0$. On the contrary, the evaluation of the integrals,
will be less efficient, as the two integrands have singularities at opposite ends of the domain, which will result in subdivisions which are only of use to one integrand. In such cases, it will be more efficient to use two sets of calls to
d01rac.
d01rac will flag extremely bad behaviour if a subinterval $\overline{\mathit{k}}$ with bounds $\left[{a}_{\overline{\mathit{k}}},{b}_{\overline{\mathit{k}}}\right]$ satisfying $\left{b}_{\overline{\mathit{k}}}{a}_{\overline{\mathit{k}}}\right<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\delta}_{a},{\delta}_{r}\times \leftba\right\right)$ has a local error estimate greater than the requested tolerance for at least one integral. The values ${\delta}_{a}$ and ${\delta}_{r}$ can be set through the optional parameters ${\mathbf{Absolute\; Interval\; Minimum}}$ and ${\mathbf{Relative\; Interval\; Minimum}}$ respectively.
9.1
Details of the Computation
This section is recommended for expert users only. It describes the contents of the arrays
com and
icom upon exit from
d01rac with
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_NOERROR,
NE_ACCURACY or
NE_QUAD_BAD_SUBDIV_INT, and provided at least one iteration completed, failure due to insufficient
licom or
lcom.
The arrays
icom and
com contain details of the integration, including various scalars, onedimensional arrays, and (effectively) twodimensional arrays. The dimensions of these arrays vary depending on the arguments and options used and the progress of the algorithm. Here we describe some of these details, including how and where they are stored in
icom and
com.
Scalar quantities:
The indices in
icom including the following scalars are available via query only options, see
Section 11.2. For example,
${I}_{\mathit{ldi}}$ is the integer value returned by the option
${\mathbf{Index\; LDI}}$.
Note the indices returned start at
$1$ and so the corresponding zero based indices require
$1$ to be subtracted as shown below. This is also true for the location of arrays within
icom and
com and consequently the indices of the elements of the one and twodimensional arrays must be modified as detailed below.
$\mathit{ldi}$ 
The leading dimension of the twodimensional integer arrays stored in icom detailed below.
$\mathit{ldi}={\mathbf{icom}}\left[{I}_{\mathit{ldi}}1\right]$. 
$\mathit{ldr}$ 
The leading dimension of the twodimensional real arrays stored in com detailed below.
$\mathit{ldr}={\mathbf{icom}}\left[{I}_{\mathit{ldr}}1\right]$. 
$\mathit{nsdiv}$ 
The number of segments that have been subdivided during the adaptive process.
$\mathit{nsdiv}={\mathbf{icom}}\left[{I}_{\mathit{nsdiv}}1\right]$. 
$\mathit{nseg}$ 
The total number of segments formed.
$\mathit{nseg}=2\mathit{nsdiv}+{s}_{\mathit{pri}}$.
$\mathit{nseg}={\mathbf{icom}}\left[{I}_{\mathit{nseg}}1\right]$. 
$\mathit{dsp}$ 
The reference of the first element of the array $\mathit{ds}$ stored in com.
$\mathit{dsp}={\mathbf{icom}}\left[{I}_{\mathit{dsp}}1\right]$. 
$\mathit{esp}$ 
The reference of the first element of the array $\mathit{es}$ stored in com.
$\mathit{esp}={\mathbf{icom}}\left[{I}_{\mathit{esp}}1\right]$. 
$\mathit{evalsp}$ 
The reference of the first element of the array $\mathit{evals}$ stored in icom.
$\mathit{evalsp}={\mathbf{icom}}\left[{I}_{\mathit{evalsp}}1\right]$. 
$\mathit{fcp}$ 
The reference of the first element of the array $\mathit{fcount}$ stored in icom.
$\mathit{fcp}={\mathbf{icom}}\left[{I}_{\mathit{fcp}}1\right]$. 
$\mathit{sinforp}$ 
The reference of the first element of the array $\mathit{sinfor}$ stored in com.
$\mathit{sinforp}={\mathbf{icom}}\left[{I}_{\mathit{sinforp}}1\right]$. 
$\mathit{sinfoip}$ 
The reference of the first element of the array $\mathit{sinfoi}$ stored in icom.
$\mathit{sinfoip}={\mathbf{icom}}\left[{I}_{\mathit{sinfoip}}1\right]$. 
Onedimensional arrays:

$\mathit{fcount}\left[{n}_{i}\right]$

$\mathit{fcount}\left[0\right]={\mathbf{icom}}\left[\mathit{fcp}1\right]$.
$\mathit{fcount}\left[j1\right]$ contains the number of different approximations of integral $\mathit{j}$ calculated, for $\mathit{j}=1,2,\dots ,{n}_{i}$.
Twodimensional arrays:

$\mathit{sinfoi}\left[5\times \mathit{nseg}\right]$

$\mathit{sinfoi}\left[0\right]={\mathbf{icom}}\left[\mathit{sinfoip}1\right]$.
$\mathit{sinfoi}$ contains information about the hierarchy of splitting.
$\mathit{sinfoi}\left[\left(\mathit{k}1\right)\times \mathit{ldi}\right]$ contains the split identifier for segment $\mathit{k}$, for $\mathit{k}=1,2,\dots ,\mathit{nseg}$.
$\mathit{sinfoi}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+1\right]$ contains the parent segment number of segment $\mathit{k}$ (i.e., the segment was split to create segment $\mathit{k}$), for $\mathit{k}=1,2,\dots ,\mathit{nseg}$.
$\mathit{sinfoi}\left[\left(k1\right)\times \mathit{ldi}+2\right]$ and
$\mathit{sinfoi}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+3\right]$
contain the segment numbers of the two child segments formed from segment $\mathit{k}$, if segment $\mathit{k}$ has been split. If segment $\mathit{k}$ has not been split, these will be negative.
$\mathit{sinfoi}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+4\right]$ contains the level at which the segment exists, corresponding to ${n}_{a}+1$, where ${n}_{a}$ is the number of ancestor segments of segment $\mathit{k}$, for $\mathit{k}=1,2,\dots ,\mathit{nseg}$. A negative level indicates that segment $\mathit{k}$ will not be split further, the level is then given by the absolute value of $\mathit{sinfoi}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+4\right]$.

$\mathit{sinfor}\left[2\times \mathit{nseg}\right]$

$\mathit{sinfor}\left[0\right]={\mathbf{com}}\left[\mathit{sinforp}1\right]$.
$\mathit{sinfor}$ contains the bounds of each segment.
$\mathit{sinfor}\left[\left(\mathit{k}1\right)\times \mathit{ldr}\right]$ contains the lower bound of segment $\mathit{k}$, for $\mathit{k}=1,2,\dots ,\mathit{nseg}$.
$\mathit{sinfor}\left[\left(\mathit{k}1\right)\times \mathit{ldr}+1\right]$ contains the upper bound of segment $\mathit{k}$, for $\mathit{k}=1,2,\dots ,\mathit{nseg}$.

$\mathit{evals}\left[{n}_{i}\times \mathit{nseg}\right]$

$\mathit{evals}\left[0\right]={\mathbf{icom}}\left[\mathit{evalsp}1\right]$.
$\mathit{evals}$ contains information to indicate whether an estimate of the integral $\mathit{j}$ has been obtained over segment $\mathit{k}$, and if so whether this evaluation still contributes to the direct estimate of ${F}_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,{n}_{i}$ and $\mathit{k}=1,2,\dots ,\mathit{nseg}$.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=0$ indicates that integral $j$ has not been evaluated over segment $\mathit{k}$.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=1$ indicates that integral $j$ has been evaluated over segment $\mathit{k}$, and that this evaluation contributes to the direct estimate of ${F}_{j}$.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=2$ indicates that integral $j$ has been evaluated over segment $\mathit{k}$, that this evaluation contributes to the direct estimate of ${F}_{j}$, and that you have requested no further evaluation of this integral at this segment by setting ${\mathbf{needi}}\left[j1\right]<0$.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=3$ indicates that integral $j$ has been evaluated over segment $\mathit{k}$, and this evaluation no longer contributes to the direct estimate of ${F}_{j}$.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=4$ indicates that integral
$j$ has been evaluated over segment
$\mathit{k}$, that this evaluation contributes to the direct estimate of
${F}_{j}$, and that this segment is too small for any further splitting to be performed. Integral
$j$ also has a local error estimate over this segment above the requested tolerance. Such segments cause
d01rac to return
${\mathbf{fail}}\mathbf{.}\mathbf{code}=$ NE_QUAD_BAD_SUBDIV_INT, indicating extremely bad behaviour.
$\mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+\mathit{j}1\right]=5$ indicates that integral $j$ has been evaluated over segment $\mathit{k}$, that this evaluation contributes to the direct estimate of ${F}_{j}$, and that this segment is too small for any further splitting to be performed. The local error estimate is however below the requested tolerance.

$\mathit{ds}\left[{n}_{i}\times \mathit{nseg}\right]$

$\mathit{ds}\left[0\right]={\mathbf{com}}\left[\mathit{dsp}1\right]$.
$\mathit{ds}\left[\left(\mathit{k}1\right)\times \mathit{ldr}+j1\right]$ contains the definite integral estimate of the $j$th integral over the $\mathit{k}$th segment, ${\mathit{ds}}_{j,\mathit{k}}$, provided it has been evaluated, for $\mathit{j}=1,2,\dots ,{n}_{i}$ and $\mathit{k}=1,2,\dots ,\mathit{nseg}$.

$\mathit{es}\left[{n}_{i}\times \mathit{nseg}\right]$

$\mathit{es}\left[0\right]={\mathbf{com}}\left[\mathit{esp}1\right]$.
$\mathit{es}\left[\left(\mathit{k}1\right)\times \mathit{ldr}+j1\right]$ contains the definite integral error estimate of the $j$th integral over the $\mathit{k}$th segment, ${\mathit{es}}_{j,\mathit{k}}$, provided it has been evaluated, for $\mathit{j}=1,2,\dots ,{n}_{i}$ and $\mathit{k}=1,2,\dots ,\mathit{nseg}$.
For each integral
$j$, the direct approximation
${D}_{j}$ of
${F}_{j}$, and its error estimate
${E}_{j}$, may be constructed as,
where
${K}_{j}$ is the set of all contributing segments,
${K}_{j}=\left\{\mathit{k}\mid \mathit{evals}\left[\left(\mathit{k}1\right)\times \mathit{ldi}+j1\right]=1\text{,}2\text{,}4\text{ or}5\text{,}1\le \mathit{k}\le \mathit{nseg}\right\}$.
${D}_{j}$ will have been returned in
${\mathbf{dinest}}\left[j1\right]$, unless extrapolation was successful, as indicated by
${\mathbf{needi}}\left[j1\right]$.
Similarly,
${E}_{j}$ will have been returned in
${\mathbf{errest}}\left[j1\right]$ unless extrapolation was successful, in which case the error estimate from the extrapolation will have been returned. If for a given integral
$j$ one or more contributing segments have unacceptable error estimates, it may be possible to improve the direct approximation by replacing the contributions from these segments with more accurate estimates should these be calculable by some means. Indeed for any segment
$\overline{\mathit{k}}\in \mathit{k}$, with lower bound
${a}_{\overline{\mathit{k}}}=\mathit{sinfor}\left(1,\overline{\mathit{k}}\right)$ and upper bound
${b}_{\overline{\mathit{k}}}=\mathit{sinfor}\left(2,\overline{\mathit{k}}\right)$, one may alter the direct approximation
${D}_{j}$ by the following,
The error estimate ${E}_{j}$ may be altered similarly.
10
Example
10.1
Program Text
10.2
Program Data
None.
10.3
Program Results
11
Optional Parameters
This section can be skipped if you wish to use the default values for all optional parameters, otherwise, the following is a list of the optional parameters available. A full description of each optional parameter is provided in
Section 11.1.
The following optional parameters, see
Section 11.2, may be utilized by expert users in conjunction with the information provided in
Section 9.1.
 ${\mathbf{Index\; LDI}}$
 ${\mathbf{Index\; LDR}}$
 ${\mathbf{Index\; NSDIV}}$
 ${\mathbf{Index\; NSEG}}$
 ${\mathbf{Index\; FCP}}$
 ${\mathbf{Index\; EVALSP}}$
 ${\mathbf{Index\; DSP}}$
 ${\mathbf{Index\; ESP}}$
 ${\mathbf{Index\; SINFOIP}}$
 ${\mathbf{Index\; SINFORP}}$
11.1
Description of the Optional Parameters
For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
 the keywords, where the minimum abbreviation of each keyword is underlined;
 a parameter value,
where the letters $a$, $i$ and $r$ denote options that take character, integer and real values respectively;
 the default value.
The following
symbols represent
various machine constants:
 $\epsilon $ represents the machine precision (see X02AJC);
 ${e}_{\mathrm{max}}$ is the maximum exponent argument of the model of floatingpoint arithmetic, (see X02BLC);
 ${R}_{\mathrm{max}}$ represents the largest representable real value (see X02ALC).
All options accept the value ‘DEFAULT’ in order to return single options to their default states.
Keywords and character values are case insensitive, however they must be separated by at least one space.
Unsetable options will return the appropriate value when calling
d01zlc. They will have no effect if passed to
d01zkc.
For
d01rac the maximum length of the argument
cvalue used by
d01zlc is
$15$.
Absolute Interval Minimum  $r$  Default $\text{}=128.0\epsilon $ 
$r={\delta}_{a}$, the absolute lower limit for a segment to be considered for subdivision. See also
${\mathbf{Relative\; Interval\; Minimum}}$ and
Section 9.
Constraint: $r\ge 128\epsilon $.
Absolute Tolerance  $r$  Default $\text{}=1024\epsilon $ 
$r={\epsilon}_{a}$, the absolute tolerance required. See also
${\mathbf{Relative\; Tolerance}}$ and
Section 3.
Constraint: $r\ge 0.0$.
$a$  Default $\text{}=\mathrm{ON}$ 
Activate or deactivate the use of the
$\epsilon $ algorithm (
Wynn (1956)).
${\mathbf{Extrapolation}}$ often reduces the number of iterations required to achieve the desired solution, but it can occasionally lead to premature convergence towards an incorrect answer.
 $\mathrm{ON}$
 Use extrapolation.
 $\mathrm{OFF}$
 Disable extrapolation.
$r$  Default $\text{}=\text{1.0e\u221212}$ 
$r={\epsilon}_{\mathit{safe}}$. If ${\epsilon}_{q}$ is the estimated error from the quadrature evaluation alone, and ${\epsilon}_{\mathit{ex}}$ is the error estimate determined using extrapolation, then the extrapolated solution will only be accepted if ${\epsilon}_{\mathit{safe}}{\epsilon}_{q}\le {\epsilon}_{\mathit{ex}}$.
Maximum Subdivisions  $i$  Default $\text{}=50$ 
$i={\mathrm{max}}_{\mathit{sdiv}}$, the maximum number of subdivisions the algorithm may use in the adaptive phase, forming at most an additional $\left(2\times {\mathrm{max}}_{\mathit{sdiv}}\right)$ segments.
Primary Divisions  $i$  Default $\text{}=1$ 
$i={s}_{\mathit{pri}}$, the number of initial segments of the domain $\left[a,b\right]$. By default the initial segment is the entire domain.
Constraint: $0<i<1000000$.
Primary Division Mode  $a$  Default $\text{}=\mathrm{AUTOMATIC}$ 
Determines how the initial set of ${s}_{\mathit{pri}}$ segments will be generated.
 AUTOMATIC
 d01rac will automatically generate ${s}_{\mathit{pri}}$ segments of equal size covering the interval $\left[a,b\right]$.
 MANUAL
 d01rac will use the breakpoints ${x}_{\mathit{i}}^{0}$, for $\mathit{i}=1,2,\dots ,{s}_{\mathit{pri}}1$, supplied in x on initial entry to generate the initial segments covering $\left[a,b\right]$. These may be supplied in any order, however it will be more efficient to supply them in ascending (or descending if $a>b$) order. Repeated breakpoints are allowed, although this will generate fewer initial segments.
Note: an absolute bound on the size of an initial segment of $10.0\epsilon $ is automatically applied in all cases, and will result in fewer initial subdivisions being generated if automatically generated or supplied breakpoints result in segments smaller than this.
Prioritize Error  $a$  Default $\text{}=\mathrm{LEVEL}$ 
Indicates how new subdivisions of segments sustaining unacceptable local errors for integrals should be prioritized.
 $\text{LEVEL}$
 Segments with lower level with unsatisfactory error estimates will be chosen over segments with greater error on higher levels. This will probably lead to more integrals being improved in earlier iterations of the algorithm, and hence will probably lead to fewer repeated returns (see argument sid), and to more integrals being satisfactorily estimated if computational exhaustion occurs.
 $\text{MAXERR}$
 The segment with the worst overall error will be split, regardless of level. This will more rapidly improve the worst integral estimates, although it will probably result in the fewest integrals being improved in earlier iterations, and may hence lead to more repeated returns (see argument sid), and potentially fewer integrals satisfying the requested tolerances if computational exhaustion occurs.
Quadrature Rule  $a$  Default $\text{}=\text{GK15}$ 
The basic quadrature rule to be used during the integration. Currently
$6$ Gauss–Kronrod rules are available, all identifiable by the letters GK followed by the number of points required by the Kronrod rule. Higher order rules generally provide higher accuracy with fewer subdivisons. However, for integrands with sharp singularities, lower order rules may be more efficient, particularly if the integrand away from the singularity is well behaved. With higher order rules, you may need to increase the
${\mathbf{Absolute\; Interval\; Minimum}}$ and the
${\mathbf{Relative\; Interval\; Minimum}}$ to maintain numerical difference between the abscissae and the segment bounds.
 GK15
 The Gauss–Kronrod rule based on $7$ Gauss points and $15$ Kronrod points.
 GK21
 The Gauss–Kronrod rule based on $10$ Gauss points and $21$ Kronrod points. This is the rule used by
d01sjc.
 GK31
 The Gauss–Kronrod rule based on $15$ Gauss points and $31$ Kronrod points.
 GK41
 The Gauss–Kronrod rule based on $20$ Gauss points and $41$ Kronrod points.
 GK51
 The Gauss–Kronrod rule based on $25$ Gauss points and $51$ Kronrod points.
 GK61
 The Gauss–Kronrod rule based on $30$ Gauss points and $61$ Kronrod points. This is the highest order rule, most suitable for highly oscilliatory integrals.
Relative Interval Minimum  $r$  Default $\text{}=\text{1.0e\u22126}$ 
$r={\delta}_{r}$, the relative factor in the lower limit,
${\delta}_{r}\leftba\right$, for a segment to be considered for subdivision. See also
${\mathbf{Absolute\; Interval\; Minimum}}$ and
Section 9.
Constraint: $r\ge 0.0$.
Relative Tolerance  $r$  Default $\text{}=\sqrt{\epsilon}$

$r={\epsilon}_{r}$, the required relative tolerance. See also
${\mathbf{Absolute\; Tolerance}}$ and
Section 3.
Constraint: $r\ge 0.0$.
Note: setting both ${\epsilon}_{r}={\epsilon}_{a}=0.0$ is possible, although it will most likely result in an excessive amount of computational effort.
11.2
Diagnostic Options
These options are provided for expert users who wish to examine and modify the precise details of the computation. They should only be used
after d01rac returns, as opposed to the options listed in
Section 11.1 which must be used
before the first call to
d01rac.
${I}_{\mathit{ldi}}$, the index of
icom required for obtaining
$\mathit{ldi}$. See
Section 9.1.
${I}_{\mathit{ldr}}$, the index of
icom required for obtaining
$\mathit{ldr}$. See
Section 9.1.
${I}_{\mathit{nsdiv}}$, the index of
icom required for obtaining
$\mathit{nsdiv}$. See
Section 9.1.
${I}_{\mathit{nseg}}$, the index of
icom required for obtaining
$\mathit{nseg}$. See
Section 9.1.
${I}_{\mathit{fcp}}$, the index of
icom required for obtaining
$\mathit{fcp}$. See
Section 9.1.
Index EVALSP  $i$  query only 
${I}_{\mathit{evalsp}}$, the index of
icom required for obtaining
$\mathit{evalsp}$. See
Section 9.1.
${I}_{\mathit{dsp}}$, the index of
icom required for obtaining
$\mathit{dsp}$. See
Section 9.1.
${I}_{\mathit{esp}}$, the index of
icom required for obtaining
$\mathit{esp}$. See
Section 9.1.
Index SINFOIP  $i$  query only 
${I}_{\mathit{sinfoip}}$, the index of
icom required for obtaining
$\mathit{sinfoip}$. See
Section 9.1.
Index SINFORP  $i$  query only 
${I}_{\mathit{sinforp}}$, the index of
icom required for obtaining
$\mathit{sinforp}$. See
Section 9.1.