NAG Library Routine Document
D01RAF
Note: this routine 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 9 of this document. If, however, you wish to reset some or all of the settings please refer to Section 10 for a detailed description of the specification of the optional parameters.
1 Purpose
D01RAF 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’. D01RAF is particularly effective for the integration of a vector of similar functions.
2 Specification
SUBROUTINE D01RAF ( 
IREVCM, NI, A, B, SID, NEEDI, X, LENX, NX, FM, LDFM, DINEST, ERREST, IOPTS, OPTS, ICOMM, LICOMM, COMM, LCOMM, IFAIL) 
INTEGER 
IREVCM, NI, SID, NEEDI(NI), LENX, NX, LDFM, IOPTS(100), ICOMM(LICOMM), LICOMM, LCOMM, IFAIL 
REAL (KIND=nag_wp) 
A, B, X(LENX), FM(LDFM,*), DINEST(NI), ERREST(NI), OPTS(100), COMM(LCOMM) 

3 Description
D01RAF is an extension to various QUADPACK routines, including QAG, QAGS and QAWSE. The extensions made allow multiple integrands to be evaluated simultaneously, using a vectorized interface and reverse communication.
The quadrature scheme employed by D01RAF 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).
D01RAF is the integration routine in the suite of routines
D01RAF,
D01RBF,
D01RCF and
D01ZKF.
First, the option arrays
IOPTS and
OPTS must be initialized using
D01ZKF. Thereafter any required options must be set before calling D01RAF, or the routines
D01RBF and
D01RCF.
The options available are described in
Section 10.
A typical usage of this suite of routines is (in psuedocode for clarity),
Setup phase
liopts = 100
lopts = 100
allocate(iopts(liopts),opts(lopts))
call D01ZKF('initialize = d01raf',iopts,liopts,opts,lopts,
ifail)
call D01ZKF('option = value',iopts,liopts,opts,lopts,
ifail)
...
call D01RCF(ni,lenxrq,ldfmrq,sdfmrq,licmin,licmax,lcmin,lcmax,iopts,opts,
ifail)
lenx = lenxrq
ldfm = ldfmrq
sdfm = sdfmrq
licomm = licmax
lcomm = lcmax
allocate(icomm(licomm),comm(lcomm),x(lenx),fm(ldfm,sdfm),needi(ni),
dinest(ni),errest(ni))
Solve phase
irevcm = 1
while irevcm≠0
call D01RAF(irevcm,ni,a,b,sid,needi,x,lenx,nx,fm,ldfm, &
dinest,errest,iopts,opts,icomm,licomm,comm, &
lcomm,ifail)
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 ifail
end select
end while
Diagnostic phase
call D01ZLF('option',ivalue,rvalue,cvalue,optype,iopts,opts,ifail)
...
call D01RBF(monit,ni,dinest,errest,iopts,opts,icomm,licomm,comm,lcomm, &
iuser,ruser,ifail)
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
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 routine enters the adaptive phase. The estimated errors are tested against the requested tolerances
${\epsilon}_{a}$ and
${\epsilon}_{r}$, corresponding to the
Absolute Tolerance and
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
Prioritize Error to favour segments with a lower level over segments with a high level and potentially worse error estimate or not. Up to
${\mathrm{max}}_{\mathit{sdiv}}$ subdivisions are allowed provided sufficient memory, where
${\mathrm{max}}_{\mathit{sdiv}}$ is the value of
Maximum Subdivisions.
Once a sufficient number of error estimates have been constructed for a particular integral, the routine may optionally use
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
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.
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 Parameters
Note: this routine uses
reverse communication. Its use involves an initial entry, intermediate exits and reentries, and a final exit, as indicated by the parameter
IREVCM. Between intermediate exits and reentries,
all parameters other than IREVCM, NEEDI and FM must remain unchanged.
 1: IREVCM – INTEGERInput/Output

On initial entry:
${\mathbf{IREVCM}}=1$.
 ${\mathbf{IREVCM}}=1$
 Set up data structures in ICOMM and COMM and start a new integration.
Constraint:
${\mathbf{IREVCM}}=1$ on initial entry.
On intermediate exit:
${\mathbf{IREVCM}}=11$ or
$12$.
IREVCM requests the values of
${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(i\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 inital estimates of the definite integrals.
If ${\mathbf{NEEDI}}\left(j\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(j\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}\right)$, for
$\mathit{j}=1,2,\dots ,{n}_{i}$, will be set to this negative value and
${\mathbf{IFAIL}}={{\mathbf{1}}}$ 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(j\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(j\right)=1$, ${f}_{j}\left({x}_{i}\right)$ must be returned.
If ${\mathbf{NEEDI}}\left(j\right)=2$, $3$ or $4$, the current error estimate of integral $j$ does not require integrand $j$ to be evaluated. Should you choose to, integrand $j$ can be evaluated in which case ${\mathbf{NEEDI}}\left(j\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 either
${\mathbf{IFAIL}}={\mathbf{1}}$,
${\mathbf{2}}$ or
${\mathbf{3}}$ will be returned and the elements of
NEEDI will reflect the current state of of the adaptive process.
On intermediate reentry:
IREVCM should normally be left unchanged. However, if
IREVCM is set to a negative value, D01RAF 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.
 2: NI – INTEGERInput
On entry: ${n}_{i}$, number of integrands.
 3: A – REAL (KIND=nag_wp)Input
On entry: $a$, the lower bound of the domain.
 4: B – REAL (KIND=nag_wp)Input
On entry:
$b$, the upper bound of the domain.
If
$\leftba\right<10\epsilon $, where
$\epsilon $ is the
machine precision (see
X02AJF), then D01RAF will return
${\mathbf{DINEST}}\left(\mathit{j}\right)={\mathbf{ERREST}}\left(\mathit{j}\right)=0.0$, for
$\mathit{j}=1,2,\dots ,{n}_{i}$.
 5: SID – INTEGEROutput
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 D01RAF return the same value of
SID, i.e., the same set of abscissae was used.
 6: NEEDI(NI) – INTEGER arrayInput/Output
On initial entry: need not be set.
On intermediate exit:
${\mathbf{NEEDI}}\left(j\right)$ indicates what action must be taken for integral
$j=1,2,\dots {n}_{i}$ (see
IREVCM).
 ${\mathbf{NEEDI}}\left(j\right)=0$
 Do not provide ${f}_{j}\left({x}_{i}\right)$. Any provided values will be ignored.
 ${\mathbf{NEEDI}}\left(j\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(j\right)=2$
 The values ${f}_{j}\left({x}_{i}\right)$ are not definitely 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(j\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(j\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 integral 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(j\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(j\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(j\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(j\right)=1$. Integral and error estimation will be performed again during the next iteration.
On intermediate reentry:
${\mathbf{NEEDI}}\left(j\right)$ may be used to indicate what action you have taken for integral
$j$.
 ${\mathbf{NEEDI}}\left(j\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(j\right)<0$
 You are abandoning the evaluation of integral $j$. The current values of ${\mathbf{DINEST}}\left(j\right)$ and ${\mathbf{ERREST}}\left(j\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(j\right)$ indicates the final state of integral
$j$.
 ${\mathbf{NEEDI}}\left(j\right)=0$
 The error estimate for ${F}_{j}$ is below the requested tolerance.
 ${\mathbf{NEEDI}}\left(j\right)=1$
 The error estimate for ${F}_{j}$ is below the requested tolerance after extrapolation.
 ${\mathbf{NEEDI}}\left(j\right)=2$
 The error estimate for ${F}_{j}$ is above the requested tolerance.
 ${\mathbf{NEEDI}}\left(j\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(j\right)<0$
 You prohibited further evaluation of integral $j$.
 7: X(LENX) – REAL (KIND=nag_wp) arrayInput/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 may be used to supply a set of initial ‘break points’ inside the domain of integration. Specifically,
${\mathbf{X}}\left(i\right)$ must contain a break point
${x}_{\mathit{i}}^{0}$, for
$\mathit{i}=1,2,\dots ,\left({s}_{\mathit{pri}}1\right)$, where
${s}_{\mathit{pri}}$ is the number of
Primary Divisions.
Constraint:
if break points 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}\right)$ is the abscissa ${x}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,{n}_{x}$, at which the appropriate integrals must be evaluated.
 8: LENX – INTEGERInput
On entry: the dimension of the array
X as declared in the (sub)program from which D01RAF is called. 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 10).
$\mathit{lenxrq}$ is returned as
LENXRQ from
D01RCF.
 9: NX – INTEGERInput/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: FM(LDFM,$*$) – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array
FM
must be at least
$\mathit{sdfmrq}$, where
$\mathit{sdfmrq}$ is dependent upon
${n}_{i}$ and the options currently set.
$\mathit{sdfmrq}$ is returned as
SDFMRQ from
D01RCF. 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(j\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: LDFM – INTEGERInput
On entry: the first dimension of the array
FM as declared in the (sub)program from which D01RAF is called.
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
D01RCF. If default options are chosen,
$\mathit{ldfmrq}={n}_{i}$, implying
${\mathbf{LDFM}}\ge {\mathbf{NI}}$.
 12: DINEST(NI) – REAL (KIND=nag_wp) arrayInput/Output
${\mathbf{DINEST}}\left(j\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: ERREST(NI) – REAL (KIND=nag_wp) arrayInput/Output
${\mathbf{ERREST}}\left(j\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: IOPTS($100$) – INTEGER arrayInput
On entry: integer option array generated by
D01ZKF.
Constraint:
IOPTS must not be altered between calls to D01RAF,
D01RCF,
D01ZKF and
D01ZLF.
 15: OPTS($100$) – REAL (KIND=nag_wp) arrayInput
On entry: real option array generated by
D01ZKF.
Constraint:
OPTS must not be altered between calls to D01RAF,
D01RCF,
D01ZKF and
D01ZLF.
 16: ICOMM(LICOMM) – INTEGER arrayCommunication Array
The contents of this array
must not be changed directly once the integration process has started.
ICOMM 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. You are free to examine the contents of
ICOMM using the diagnostic routine
D01RBF at any time during, or after, the integration process.
 17: LICOMM – INTEGERInput
On entry: the dimension of the array
ICOMM.
Constraint:
${\mathbf{LICOMM}}\ge \mathit{licmin}$, where
$\mathit{licmin}$ is dependent upon
NI and the current options set.
$\mathit{licmin}$ is returned as
LICMIN from
D01RCF. If the default options are set, then
$\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
D01RCF. If default options are chosen, except for possibly increasing the value of
${s}_{\mathit{pri}}$, then
$\mathit{licmax}=50+5\times {\mathbf{NI}}+\left({s}_{\mathit{pri}}+100\right)\times \left(5+{\mathbf{NI}}\right)$.
 18: COMM(LCOMM) – REAL (KIND=nag_wp) arrayCommunication Array
The contents of this array
must not be changed directly once the integration process has started.
COMM 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. You are free to examine the contents of
COMM using the diagnostic routine
D01RBF at any time during, or after, the solution phase.
 19: LCOMM – INTEGERInput
On entry: the dimension of the array
COMM.
Constraint:
${\mathbf{LCOMM}}>\mathit{lcmin}$, where
$\mathit{lcmin}$ is dependent upon
NI,
${s}_{\mathit{pri}}$ and the current options set.
$\mathit{lcmin}$ is returned as
LCMIN from
D01RCF. If default options are set, then
$\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 parameters, the maximum value,
$\mathit{lcmax}$, of
LCOMM that may be required, is returned as
LCMAX from
D01RCF. 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: IFAIL – INTEGERInput/Output
On initial entry:
IFAIL must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction 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 parameters 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 final 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}}={\mathbf{0}}$ or
${{\mathbf{1}}}$, explanatory error messages are output on the current error message unit (as defined by
X04AAF).
Note: D01RAF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
 ${\mathbf{IFAIL}}=1$

At least one error estimate exceeded the requested tolerances.
 ${\mathbf{IFAIL}}=2$

Extremely bad behaviour was detected for at least one integral.
 ${\mathbf{IFAIL}}=3$

Extremely bad behaviour was detected for at least one integral. At least one other integral error estimate was above the requested tolerance.
 ${\mathbf{IFAIL}}=11$

IREVCM had an illegal value.
On entry,
${\mathbf{IREVCM}}=\u27e8\mathit{\text{value}}\u27e9$.
 ${\mathbf{IFAIL}}=21$

On entry, ${\mathbf{NI}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{NI}}\ge 1$.
 ${\mathbf{IFAIL}}=71$

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

LENX is insufficient for the chosen options.
On entry,
${\mathbf{LENX}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint:
${\mathbf{LENX}}\ge \u27e8\mathit{\text{value}}\u27e9$.
 ${\mathbf{IFAIL}}=111$

${\mathbf{LDFM}}<\mathit{ldfmrq}$. If default options are chosen, this implies ${\mathbf{LDFM}}<{\mathbf{NI}}$.
On entry, ${\mathbf{LDFM}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint: ${\mathbf{LDFM}}\ge \u27e8\mathit{\text{value}}\u27e9$.
 ${\mathbf{IFAIL}}=171$

The length of
LICOMM is insufficient for additional subdivision.
On entry,
${\mathbf{LICOMM}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint:
${\mathbf{LICOMM}}\ge \u27e8\mathit{\text{value}}\u27e9$.
 ${\mathbf{IFAIL}}=191$

The length of
LCOMM is insufficient for additional subdivision.
On entry,
${\mathbf{LCOMM}}=\u27e8\mathit{\text{value}}\u27e9$.
Constraint:
${\mathbf{LCOMM}}\ge \u27e8\mathit{\text{value}}\u27e9$.
 ${\mathbf{IFAIL}}=1001$

Either the option arrays
IOPTS and
OPTS have not been initialized for D01RAF, or they have become corrupted.
 ${\mathbf{IFAIL}}=1101$

On entry, one of
ICOMM and
COMM has become corrupted.
 ${\mathbf{IFAIL}}=1$

Evaluation of all integrals has been stopped during the initial phase.
7 Accuracy
D01RAF 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
Absolute Tolerance and
Relative Tolerance respectively. Moreover, it returns
ERREST, the entries of which in normal circumstances satisfy,
The time required by D01RAF is usually dominated by the time required to evaluate the values of the integrands ${f}_{j}$.
D01RAF 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 D01RAF.
The diagnostic routine
D01RBF may be used to examine the subdivision strategy used by D01RAF both during intermediate stages and after completion. In particular it may be used to examine how individual integrals behaved over individual subintervals, and so may be used to determine where integrals behaved extremely badly should D01RAF return with
${\mathbf{IFAIL}}={\mathbf{2}}$ or
${\mathbf{3}}$.
D01RAF will flag extremely bad behaviour if a subinterval
$\stackrel{}{k}$ with bounds
$\left[{a}_{\stackrel{}{k}},{b}_{\stackrel{}{k}}\right]$ satisfying
$\left{b}_{\stackrel{}{k}}{a}_{\stackrel{}{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
Absolute Interval Minimum and
Relative Interval Minimum respectively. You may try using D01RAF again over such subintervals with only the integrals that exhibited bad behaviour and a new set of communication arrays, or you may use another quadrature routine designed to deal with any known specific singular behaviour over this subinterval. If these provide sufficient accuracy over the specific subintervals, you may correct the definite integral evaluations over
$\left[a,b\right]$ by subtracting the inadequate subinterval estimations from
DINEST and
ERREST and adding the acceptable subinterval estimations generated by the direct integration of the subinterval where bad behaviour was detected. See
Section 8 in D01RBF.
9 Example
9.1 Program Text
Program Text (d01rafe.f90)
9.2 Program Data
None.
9.3 Program Results
Program Results (d01rafe.r)
10 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 10.1.
10.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\text{ and}r$ denote options that take character, integer and real values respectively;
 the default value.
The following
symbol represents
various machine constants:
 $\epsilon $ represents the machine precision, typically ${2}^{53}$ for $64$bit real numbers (see X02AJF).
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.
For D01RAF the maximum length of the parameter
CVALUE used by
D01ZLF 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
Relative Interval Minimum and
Section 8.
Constraint: $r\ge 128\epsilon $.
Absolute Tolerance  $r$  Default $\text{}=1024\epsilon $ 
$r={\epsilon}_{a}$, the absolute tolerance required. See also
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)).
Extrapolation can occasionally lead to premature convergence of the integrations.
 $\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{}=\text{AUTOMATIC}$ 
Determines how the initial set of ${s}_{\mathit{pri}}$ segments will be generated.
 AUTOMATIC
 D01RAF will automatically generate ${s}_{\mathit{pri}}$ segments of equal size covering the interval $\left[a,b\right]$.
 MANUAL
 D01RAF will use the break points ${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 break points 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{}=\text{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 parameter 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 parameter 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
Absolute Interval Minimum and the
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
D01ATF
 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
Absolute Interval Minimum and
Section 8.
Constraint: $r\ge 0.0$.
Relative Tolerance  $r$  Default $\text{}=\sqrt{\epsilon}$ 
$r={\epsilon}_{r}$, the required relative tolerance. See also
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.