D01RAF (PDF version)
D01 Chapter Contents
D01 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

D01RAF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.
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 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.

 Contents

    1  Purpose
    7  Accuracy

1  Purpose

D01RAF is a general purpose adaptive integrator which calculates an approximation to a vector of definite integrals F over a finite range a,b , given the vector of integrands fx.
F = a b fx d x  
If the same subdivisions of the range are equally good for functions f1x and f2x, because f1x and f2x have common areas of the range where they vary slowly and where they vary quickly, then we say that f1x and f2x 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, ICOM, LICOM, COM, LCOM, IFAIL)
INTEGER  IREVCM, NI, SID, NEEDI(NI), LENX, NX, LDFM, IOPTS(100), ICOM(LICOM), LICOM, LCOM, IFAIL
REAL (KIND=nag_wp)  A, B, X(LENX), FM(LDFM,*), DINEST(NI), ERREST(NI), OPTS(100), COM(LCOM)

3  Description

D01RAF 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 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 and D01RCF. It also uses optional parameters, which can be set and queried using the routines D01ZKF and D01ZLF respectively. The options available are described in Section 11.
First, the option arrays IOPTS and OPTS must be initialized using D01ZKF. Thereafter any required options must be set before calling D01RAF, or the routine D01RCF.
A typical usage of this suite of routines is (in pseudo-code 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
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
    call D01RAF(irevcm,ni,a,b,sid,needi,x,lenx,nx,fm,ldfm,   &
                dinest,errest,iopts,opts,icom,licom,com, &
                lcom,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)
...
During the initial solve phase, the first estimation of the definite integral and error estimate is constructed over the interval a,b . This will have been divided into spri  level 1 segments, where spri  is the number of Primary Divisions, and will use any provided break-points if Primary Division Mode=MANUAL.
Once a complete integral estimate over a,b  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 εa  and ε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 either favour the segment with the maximum error, or the segment with the lowest level supporting an unacceptable (although potentially non-maximal) error. Up to maxsdiv  subdivisions are allowed if sufficient memory is provided, where maxsdiv  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 εsafe  can be set using Extrapolation Safeguard, and the extrapolated solution will only be considered if εsafe εq εex , where εq  and ε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 DINESTj on completion of D01RAF. Otherwise the direct solution will be returned in DINESTj. This is indicated by the value of NEEDIj 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 emSn 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 re-entries, and a final exit, as indicated by the parameter IREVCM. Between intermediate exits and re-entries, all parameters other than IREVCM, NEEDI and FM must remain unchanged.
1:     IREVCM – INTEGERInput/Output
On initial entry: IREVCM=1.
IREVCM=1
Sets up data structures in ICOM and COM and starts a new integration.
Constraint: IREVCM=1 on initial entry.
On intermediate exit: IREVCM=11 or 12.
IREVCM requests the integrands fjxi be evaluated for all required j1,,ni as indicated by NEEDI, and at all the points xi, for i=1,2,,nx. Abscissae xi are provided in Xi and fjxi must be returned in FMji.
During the initial solve phase:
IREVCM=11
Function values are required to construct the initial estimates of the definite integrals.
If NEEDIj=1, fjxi must be supplied in FMji. This will be the case unless you have abandoned the evaluation of specific integrals on a previous call.
If NEEDIj=0, you have previously abandoned the evaluation of integral j, and hence should not supply the value of fjxi.
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, NEEDIj, for j=1,2,,ni, will be set to this negative value and IFAIL=-1 will be returned.
During the adaptive solve phase:
IREVCM=12
Function values are required to improve the estimates of the definite integrals.
If NEEDIj=0, any evaluation of fjxi will be discarded, so there is no need to provide them.
If NEEDIj=1, fjxi must be provided in FMji.
If NEEDIj=2, 3 or 4, the current error estimate of integral j does not require integrand j to be evaluated and provided in FMji. Should you choose to, integrand j can be evaluated in which case NEEDIj 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 IFAIL=1, 2 or 3 will be returned and the elements of NEEDI will reflect the current state of the adaptive process.
On intermediate re-entry: IREVCM should normally be left unchanged. However, if IREVCM is set to a negative value, D01RAF will terminate, (see IREVCM=11 and IREVCM=12 above).
On final exit: IREVCM=0.
IREVCM=0
Indicates the algorithm has completed.
2:     NI – INTEGERInput
On entry: ni, the 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 b-a<10ε, where ε is the machine precision (see X02AJF), then D01RAF will return DINESTj=ERRESTj=0.0, for j=1,2,,ni.
5:     SID – INTEGEROutput
For advanced users.
On intermediate exit: SID identifies a specific set of abscissae, 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 x, and reuse them should D01RAF return the same value of SID, i.e., the same set of abscissae was used.
6:     NEEDINI – INTEGER arrayInput/Output
On initial entry: need not be set.
On intermediate exit: NEEDIj indicates what action must be taken for integral j=1,2,ni (see IREVCM).
NEEDIj=0
Do not provide fjxi. Any provided values will be ignored.
NEEDIj=1
The values fjxi must be provided in FMji, for i=1,2,,nx.
NEEDIj=2
The values fjxi 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 NEEDIj=1, and supply fjxi in FMji, for i=1,2,,nx.
NEEDIj=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 fjxi 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 NEEDIj=1, and supply fjxi in FMji, for i=1,2,,nx.
NEEDIj=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 DINESTj is greatly different to what you may expect, you may force the algorithm to re-evaluate this conclusion by including the evaluations of integrand j at xi, for i=1,2,,nx, and setting NEEDIj=1. Integral and error estimation will be performed again during the next iteration.
On intermediate re-entry: NEEDIj may be used to indicate what action you have taken for integral j.
NEEDIj=1
You have provided the values fjxi in FMji, for i=1,2,,nx.
NEEDIj<0
You are abandoning the evaluation of integral j. The current values of DINESTj and ERRESTj will be returned on final completion.
Otherwise you have not provided the value fjxi.
On final exit: NEEDIj indicates the final state of integral j.
NEEDIj=0
The error estimate for Fj is below the requested tolerance.
NEEDIj=1
The error estimate for Fj is below the requested tolerance after extrapolation.
NEEDIj=2
The error estimate for Fj is above the requested tolerance.
NEEDIj=3
The error estimate for Fj is above the requested tolerance, and extremely bad behaviour of integral j has been detected.
NEEDIj<0
You prohibited further evaluation of integral j.
7:     XLENX – REAL (KIND=nag_wp) arrayInput/Output
On initial entry: if Primary Division Mode=AUTOMATIC, X need not be set. This is the default behaviour.
If Primary Division Mode=MANUAL, X is used to supply a set of initial ‘break-points’ inside the domain of integration. Specifically, Xi must contain a break-point xi0, for i=1,2,,spri-1, where spri is the number of Primary Divisions.
Constraint: if break-points are supplied, xi0a,b, xi0-a>10.0ε, xi0-b>10.0ε, for i=1,2,,spri-1.
On intermediate exit: Xi is the abscissa xi, for i=1,2,,nx, 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 LENX=max122,spri-1 will be sufficient for all cases.
Constraint: LENXlenxrq, where lenxrq is dependent upon the options currently set (see Section 11). lenxrq is returned as LENXRQ from D01RCF.
9:     NX – INTEGERInput/Output
On initial entry: need not be set.
On intermediate exit: nx, the number of abscissae at which integrands are required.
On intermediate re-entry: must not be changed.
10:   FMLDFM* – REAL (KIND=nag_wp) arrayInput
Note: the second dimension of the array FM must be at least sdfmrq, where sdfmrq is dependent upon ni and the options currently set. sdfmrq is returned as SDFMRQ from D01RCF. If default options are chosen, sdfmrq=lenxrq.
On initial entry: need not be set.
On intermediate re-entry: if indicated by NEEDIj you must supply the values fjxi in FMji, for i=1,2,,nx and j=1,2,,ni.
11:   LDFM – INTEGERInput
On entry: the first dimension of the array FM as declared in the (sub)program from which D01RAF is called.
Constraint: LDFMldfmrq, where ldfmrq is dependent upon ni and the options currently set. ldfmrq is returned as LDFMRQ from D01RCF. If default options are chosen, ldfmrq=ni, implying LDFMNI.
12:   DINESTNI – REAL (KIND=nag_wp) arrayInput/Output
DINESTj contains the current estimate of the definite integral Fj.
On initial entry: need not be set.
On intermediate re-entry: must not be altered.
On exit: contains the current estimates of the NI integrals. If IREVCM=0, this will be the final solution.
13:   ERRESTNI – REAL (KIND=nag_wp) arrayInput/Output
ERRESTj contains the current error estimate of the definite integral Fj.
On initial entry: need not be set.
On intermediate re-entry: must not be altered.
On exit: contains the current error estimates for the NI integrals. If IREVCM=0, ERREST contains the final error estimates of the NI integrals.
14:   IOPTS100 – INTEGER arrayCommunication Array
15:   OPTS100 – REAL (KIND=nag_wp) arrayCommunication Array
The arrays IOPTS and OPTS must not be altered between calls to any of the routines D01RAF, D01RCF, D01ZKF and D01ZLF.
16:   ICOMLICOM – INTEGER arrayCommunication Array
ICOM contains details of the integration procedure, including information on the integration of the ni integrals over individual segments. This data is stored sequentially in the order that segments are created. For further information see Section 9.1.
17:   LICOM – INTEGERInput
On entry: the dimension of the array ICOM.
Constraint: LICOMlicmin, where licmin is dependent upon NI and the current options set. licmin is returned as LICMIN from D01RCF. If the default options are set, then licmin=55+6×NI. Larger values than licmin are recommended if you anticipate that any integrals will require the domain to be further subdivided.
The maximum value that may be required, licmax, is returned as LICMAX from D01RCF. If default options are chosen, except for possibly increasing the value of spri, then licmax=50+5×NI+spri+100×5+NI.
18:   COMLCOM – REAL (KIND=nag_wp) arrayCommunication Array
COM contains details of the integration procedure, including information on the integration of the ni integrals over individual segments. This data is stored sequentially in the order that segments are created. For further information see Section 9.1.
19:   LCOM – INTEGERInput
On entry: the dimension of the array COM.
Constraint: LCOM>lcmin, where lcmin is dependent upon NI, spri and the current options set. lcmin is returned as LCMIN from D01RCF. If default options are set, then lcmin =96+12×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, lcmax, of LCOM that may be required, is returned as LCMAX from D01RCF. If default options are chosen, lcmax=94+9×NI+NI/2 +spri+100×2+NI/2+2×NI.
20:   IFAIL – INTEGERInput/Output
On initial entry: IFAIL must be set to 0, -1​ 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​ 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 IFAIL0 on exit, the recommended value is -1. When the value -1​ or ​1 is used it is essential to test the value of IFAIL on exit.
On final exit: IFAIL=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6  Error Indicators and Warnings

If on entry IFAIL=0 or -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:
IFAIL=1
At least one error estimate exceeded the requested tolerances.
IFAIL=2
Extremely bad behaviour was detected for at least one integral.
IFAIL=3
Extremely bad behaviour was detected for at least one integral. At least one other integral error estimate was above the requested tolerance.
IFAIL=11
IREVCM had an illegal value.
On entry, IREVCM=value.
IFAIL=21
On entry, NI=value.
Constraint: NI1.
IFAIL=71
On entry, Primary Division Mode=MANUAL and at least one supplied break-point in X is outside of the domain of integration.
IFAIL=81
LENX is insufficient for the chosen options.
On entry, LENX=value.
Constraint: LENXvalue.
IFAIL=111
LDFM<ldfmrq. If default options are chosen, this implies LDFM<NI.
On entry, LDFM=value.
Constraint: LDFMvalue.
IFAIL=171
LICOM is insufficient for additional subdivision.
On entry, LICOM=value.
Constraint: LICOMvalue.
IFAIL=191
LCOM is insufficient for additional subdivision.
On entry, LCOM=value.
Constraint: LCOMvalue.
IFAIL=1001
Either the option arrays IOPTS and OPTS have not been initialized for D01RAF, or they have become corrupted.
IFAIL=1101
On entry, one of ICOM and COM has become corrupted.
IFAIL=-1
Evaluation of all integrals has been stopped during the initial phase.
IFAIL=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.8 in the Essential Introduction for further information.
IFAIL=-399
Your licence key may have expired or may not have been installed correctly.
See Section 3.7 in the Essential Introduction for further information.
IFAIL=-999
Dynamic memory allocation failed.
See Section 3.6 in the Essential Introduction for further information.

7  Accuracy

D01RAF cannot guarantee, but in practice usually achieves, the following accuracy for each integral Fj:
Fj - DINESTj tol  
where
tol = maxεa, εr × Fj  
εa and εr are the error tolerances Absolute Tolerance and Relative Tolerance respectively. Moreover, it returns ERREST, the entries of which in normal circumstances satisfy,
Fj - DINESTj ERRESTj tol .  

8  Parallelism and Performance

D01RAF is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
D01RAF 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.

9  Further Comments

The time required by D01RAF is usually dominated by the time required to evaluate the values of the integrands fj.
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,
0 1 logx x-12 x2 dx  
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,
0 1 logx log1-x dx  
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.
D01RAF will flag extremely bad behaviour if a sub-interval k-  with bounds ak-,bk-  satisfying bk- - ak- < maxδa, δr × b-a  has a local error estimate greater than the requested tolerance for at least one integral. The values δa and δr can be set through the optional parameters Absolute Interval Minimum and 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 D01RAF with IFAIL=0, 1, 2 or 3, 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, one-dimensional arrays, and (effectively) two-dimensional arrays. The dimensions of these arrays vary depending on the parameters 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, Ildi is the integer value returned by the option Index LDI.
ldi The leading dimension of the two-dimensional integer arrays stored in ICOM detailed below.
ldi=ICOMIldi.
ldr The leading dimension of the two-dimensional real arrays stored in COM detailed below.
ldr=ICOMIldr.
nsdiv The number of segments that have been subdivided during the adaptive process.
nsdiv=ICOMInsdiv.
nseg The total number of segments formed.
nseg=2nsdiv+spri.
nseg=ICOMInseg.
dsp The reference of the first element of the array ds stored in COM.
dsp=ICOMIdsp.
esp The reference of the first element of the array es stored in COM.
esp=ICOMIesp.
evalsp The reference of the first element of the array evals stored in ICOM.
evalsp=ICOMIevalsp.
fcp The reference of the first element of the array fcount stored in ICOM.
fcp=ICOMIfcp.
sinforp The reference of the first element of the array sinfor stored in COM.
sinforp=ICOMIsinforp.
sinfoip The reference of the first element of the array sinfoi stored in ICOM.
sinfoip=ICOMIsinfoip.
One-dimensional arrays:
fcountni 
fcount1=ICOMfcp.
fcountj contains the number of different approximations of integral j calculated, for j=1,2,,ni.
Two-dimensional arrays:
sinfoi5nseg 
sinfoi11=ICOMsinfoip.
sinfoi contains information about the hierarchy of splitting.
sinfoi1k contains the split identifier for segment k, for k=1,2,,nseg.
sinfoi2k contains the parent segment number of segment k (i.e., the segment was split to create segment k), for k=1,2,,nseg.
sinfoi3k and sinfoi4k contain the segment numbers of the two child segments formed from segment k, if segment k has been split. If segment k has not been split, these will be negative.
sinfoi5k contains the level at which the segment exists, corresponding to na+1, where na is the number of ancestor segments of segment k, for k=1,2,,nseg. A negative level indicates that segment k will not be split further, the level is then given by the absolute value of sinfoi5k.
sinfor2nseg 
sinfor11=COMsinforp.
sinfor contains the bounds of each segment.
sinfor1k contains the lower bound of segment k, for k=1,2,,nseg.
sinfor2k contains the upper bound of segment k, for k=1,2,,nseg.
evalsninseg 
evals11=ICOMevalsp.
evals contains information to indicate whether an estimate of the integral j has been obtained over segment k, and if so whether this evaluation still contributes to the direct estimate of Fj, for j=1,2,,ni and k=1,2,,nseg.
evalsjk=0 indicates that integral j has not been evaluated over segment k.
evalsjk=1 indicates that integral j has been evaluated over segment k, and that this evaluation contributes to the direct estimate of Fj.
evalsjk=2 indicates that integral j has been evaluated over segment k, that this evaluation contributes to the direct estimate of Fj, and that you have requested no further evaluation of this integral at this segment by setting NEEDIj<0.
evalsjk=3 indicates that integral j has been evaluated over segment k, and this evaluation no longer contributes to the direct estimate of Fj.
evalsjk=4 indicates that integral j has been evaluated over segment k, that this evaluation contributes to the direct estimate of Fj, 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 D01RAF to return IFAIL=2 or 3, indicating extremely bad behaviour.
evalsjk=5 indicates that integral j has been evaluated over segment k, that this evaluation contributes to the direct estimate of Fj, and that this segment is too small for any further splitting to be performed. The local error estimate is however below the requested tolerance.
dsninseg 
ds11=COMdsp.
dsjk contains the definite integral estimate of the jth integral over the kth segment, dsj,k, provided it has been evaluated, for j=1,2,,ni and k=1,2,,nseg.
esninseg 
es11=COMesp.
esjk contains the definite integral error estimate of the jth integral over the kth segment, esj,k, provided it has been evaluated, for j=1,2,,ni and k=1,2,,nseg.
For each integral j, the direct approximation Dj of Fj, and its error estimate Ej, may be constructed as,
Fj Dj = Kj dsj,k , Fj-Dj Ej = Kj esj,k ,  
where Kj is the set of all contributing segments, Kj=kevalsjk=1, ​2, ​4​ or ​5, ​1knseg. Dj will have been returned in DINESTj, unless extrapolation was successful, as indicated by NEEDIj.
Similarly, Ej will have been returned in ERRESTj 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 k-k, with lower bound ak-=sinfor1,k- and upper bound bk-=sinfor2,k-, one may alter the direct approximation Dj by the following,
ds j,k- new ak- bk- fjx ​ ​ dx Dj = Kj ds j,k - ds j , k- + ds j,k- new .  
The error estimate Ej may be altered similarly.

10  Example

This example integrates
F = 0 π x sin2x cos15x x2 sin2x cos50x dx .  

10.1  Program Text

Program Text (d01rafe.f90)

10.2  Program Data

None.

10.3  Program Results

Program Results (d01rafe.r)

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.

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 following symbol represents various machine constants:
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 D01ZLF. They will have no effect if passed to D01ZKF.
For D01RAF the maximum length of the parameter CVALUE used by D01ZLF is 15.
Absolute Interval MinimumrDefault =128.0ε
r=δa, the absolute lower limit for a segment to be considered for subdivision. See also Relative Interval Minimum and Section 9.
Constraint: r128ε.
Absolute TolerancerDefault =1024ε
r=εa, the absolute tolerance required. See also Relative Tolerance and Section 3.
Constraint: r0.0.
ExtrapolationaDefault =ON
Activate or deactivate the use of the ε algorithm (Wynn (1956)). 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.
ON
Use extrapolation.
OFF
Disable extrapolation.
Extrapolation SafeguardrDefault =1.0E−12
r=εsafe. If εq is the estimated error from the quadrature evaluation alone, and εex is the error estimate determined using extrapolation, then the extrapolated solution will only be accepted if εsafe εq εex.
Maximum SubdivisionsiDefault =50
i=maxsdiv, the maximum number of subdivisions the algorithm may use in the adaptive phase, forming at most an additional 2×maxsdiv  segments.
Primary DivisionsiDefault =1
i=spri, the number of initial segments of the domain a,b. By default the initial segment is the entire domain.
Constraint: 0<i<1000000.
Primary Division ModeaDefault =AUTOMATIC
Determines how the initial set of spri segments will be generated.
AUTOMATIC
D01RAF will automatically generate spri segments of equal size covering the interval a,b.
MANUAL
D01RAF will use the break-points xi0, for i=1,2,,spri-1, supplied in X on initial entry to generate the initial segments covering a,b. 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ε is automatically applied in all cases, and will result in fewer initial subdivisions being generated if automatically generated or supplied break-points result in segments smaller than this..
Prioritize ErroraDefault=LEVEL
Indicates how new subdivisions of segments sustaining unacceptable local errors for integrals should be prioritized.
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.
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 RuleaDefault =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 MinimumrDefault =1.0E−6
r=δr, the relative factor in the lower limit, δrb-a, for a segment to be considered for subdivision. See also Absolute Interval Minimum and Section 9.
Constraint: r0.0.
Relative TolerancerDefault =ε 
r=εr, the required relative tolerance. See also Absolute Tolerance and Section 3.
Constraint: r0.0.
Note: setting both εr=ε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 D01RAF returns, as opposed to the options listed in Section 11.1 which must be used before the first call to D01RAF.
Index LDIiquery only
Ildi, the index of ICOM required for obtaining ldi. See Section 9.1.
Index LDRiquery only
Ildr, the index of ICOM required for obtaining ldr. See Section 9.1.
Index NSDIViquery only
Insdiv, the index of ICOM required for obtaining nsdiv. See Section 9.1.
Index NSEGiquery only
Inseg, the index of ICOM required for obtaining nseg. See Section 9.1.
Index FCPiquery only
Ifcp, the index of ICOM required for obtaining fcp. See Section 9.1.
Index EVALSPiquery only
Ievalsp, the index of ICOM required for obtaining evalsp. See Section 9.1.
Index DSPiquery only
Idsp, the index of ICOM required for obtaining dsp. See Section 9.1.
Index ESPiquery only
Iesp, the index of ICOM required for obtaining esp. See Section 9.1.
Index SINFOIPiquery only
Isinfoip, the index of ICOM required for obtaining sinfoip. See Section 9.1.
Index SINFORPiquery only
Isinforp, the index of ICOM required for obtaining sinforp. See Section 9.1.

D01RAF (PDF version)
D01 Chapter Contents
D01 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2015