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

NAG Library Routine Document

D01ESF

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.

▸▿ Contents

    1  Purpose
    2  Specification
▸▿ 3  Description
    4  References
    5  Parameters
    6  Error Indicators and Warnings
    7  Accuracy
▸▿ 8  Parallelism and Performance
    9  Further Comments
▸▿ 10  Example
▸▿ 11  Optional Parameters

1  Purpose

D01ESF approximates a vector of definite integrals F over the unit hypercube Ω=0,1d, given the vector of integrands fx.
F = ∫ Ω fx dx = ∫ 0 1 ∫ 0 1 … ∫ 0 1 f x1,x2,…,xd dx1 dx2 … dxd .  
The routine uses a sparse grid discretisation, allowing for computationally feasible estimations of integrals of high dimension (d∼O100).

2  Specification

SUBROUTINE D01ESF ( NI, NDIM, F, MAXDLV, DINEST, ERREST, IVALID, IOPTS, OPTS, IUSER, RUSER, IFAIL)
INTEGER  NI, NDIM, MAXDLV(NDIM), IVALID(NI), IOPTS(100), IUSER(*), IFAIL
REAL (KIND=nag_wp)  DINEST(NI), ERREST(NI), OPTS(100), RUSER(*)
EXTERNAL  F

3  Description

D01ESF uses a sparse grid to generate a vector of approximations F^ to a vector of integrals F over the unit hypercube Ω=0,1d, that is,
F^ ≈ F = ∫ 0,1 d fx dx .  

3.1  Comparing Quadrature Over Full and Sparse Grids

Before illustrating the sparse grid construction, it is worth comparing integration over a sparse grid to integration over a full grid.
Given a one-dimensional quadrature rule with N abscissae, which accurately evaluates a polynomial of order PN, a full tensor product over d dimensions, a full grid, may be constructed with Nd multidimensional abscissae. Such a product will accurately integrate a polynomial where the maximum power of any dimension is PN. For example if d=2 and PN=3, such a rule will accurately integrate any polynomial whose highest order term is x13x23. Such a polynomial may be said to have a maximum combined order of PNd, provided no individual dimension contributes a power greater than PN. However, the number of multidimensional abscissae, or points, required increases exponentially with the dimension, rapidly making such a construction unusable.
The sparse grid technique was developed by Smolyak (Smolyak (1963)). In this, multiple one-dimensional quadrature rules of increasing accuracy are combined in such a way as to provide a multidimensional quadrature rule which will accurately evaluate the integral of a polynomial whose maximum order appears as a monomial. Hence a sparse grid construction whose highest level quadrature rule has polynomial order PN will accurately integrate a polynomial whose maximum combined order is also PN. Again taking PN=3, one may, theoretically, accurately integrate a polynomial such as x3+x2y+y3, but not a polynomial such as x3y3+xy. Whilst this has a lower maximum combined order than the full tensor product, the number of abscissae required increases significantly slower than the equivalent full grid, making some classes of integrals of dimension d∼O100 tractable. Specifically, if a one-dimensional quadrature rule of level ℓ has N∼O2ℓ abscissae, the corresponding full grid will have O 2ℓ d  multidimensional abscissae, whereas the sparse grid will have O 2ℓ d ℓ-1 . Figure 1 demonstrates this using a Gauss–Patterson rule with 15 points in 3 dimensions. The full grid requires 3375 points, whereas the sparse grid only requires 111.
Three-dimensional full (left) and sparse (right) grids, constructed from the 15 point Gauss–Patterson rule
Figure 1: Three-dimensional full (left) and sparse (right) grids, constructed from the 15 point Gauss–Patterson rule

3.2  Sparse Grid Quadrature

We now include a brief description of the sparse grid construction, sufficient for the understanding of the use of this routine. For a more detailed analysis, see Gerstner and Griebel (1998).
Consider a one-dimensional nℓ-point quadrature rule of level ℓ, Qℓ. The action of this rule on a integrand f is to approximate its definite one-dimensional integral F 1  as,
F 1 = ∫ 0 1 fx dx ≈ Qℓf = ∑ i=1 nℓ wℓ,i × fxℓ,i ,  
using weights wℓ,i and abscissae xℓ,i, for i=1,2,…,nℓ.
Now construct a set of one-dimensional quadrature rules, Qℓ ∣ ℓ=1,…,L , such that the accuracy of the quadrature rule increases with the level number. In this routine we exclusively use quadrature rules which are completely nested, implying that if an abscissae xℓ,k is in level ℓ, it is also in level ℓ+1. The quantity L denotes some maximum level appropriate to the rules that have been selected.
Now define the action of the tensor product of d rules as,
Qℓ1 ⊗ ⋯ ⊗ Qℓd f = ∑ i1=1 n ℓ1 ⋯ ∑ i d =1 n ℓ d w ℓ1 , i1 ⋯ w ℓd , id f x ℓ1 , i1 ,…, x ℓd , id ,  
where the individual level indices ℓj are not necessarily ordered or unique. Each tensor product of d rules defines an action of the quadrature rules Qℓ, ℓ=ℓ1,ℓ2,…,ℓd over a subspace, which is given a level ℓ = ∑ j=1 d ℓj. If all rule levels are equal, this is the full tensor product of that level.
The sparse grid construction of level ℓ  can then be declared as the sum of all actions of the quadrature differences Δk = Qk - Qk-1 , over all subspaces having a level at most ℓ-d+1,
F d ≈ Q ℓ d f = ∑ level at most ​ℓ-d+1 Δk1 ⊗ ⋯ ⊗ Δkd f . (1)
By definition, all subspaces used for level ℓ-1 must also be used for level ℓ, and as such the difference between the result of all actions over subsequent sparse grid constructions may be used as an error estimate.
Let L be the maximum level allowable in a sparse grid construction. The classical sparse grid construction of ℓ=L allows each dimension to support a one-dimensional quadrature rule of level at most L, with such a quadrature rule being used in every dimension at least once. Such a construction lends equal weight to each dimension of the integration, and is termed here ‘isotropic’.
Define the set m=mj,j=1,2,…,d, where mj is the maximum quadrature rule allowed in the jth dimension, and mq  to be the maximum quadrature rule used by any dimension. Let a subspace be identified by its quadrature difference levels, k = k1,k2,…,kd.
The classical construction may be extended by allowing different dimensions to have different values mj, and by allowing mq≤L. This creates non-isotropic constructions. These are especially useful in higher dimensions, where some dimensions contribute more than others to the result, as they can drastically reduce the number of function evaluations required.
For example, consider the two-dimensional construction with L=4. The classical isotropic construction would have the following subspaces.
Subspaces generated by a classical sparse grid with L=4.
Level Subspaces
1 1,1 
2 2,1, 1,2 
3 3,1, 2,2, 1,3 
4 4,1, 3,2, 2,3, 1,4 
If the variation in the second dimension is sufficiently accurately described by a quadrature rule of level 2, the contributions of the subspaces 1,3 and 1,4 are probably negligible. Similarly, if the variation in the first dimension is sufficiently accurately described by a quadrature rule of level 3, the subspace 4,1 is probably negligible. Furthermore the subspace 2,3 would also probably have negligible impact, whereas the subspaces 2,2 and 3,2 would not. Hence restricting the first dimension to a maximum level of 3, and the second dimension to a maximum level of 2 would probably give a sufficiently acceptable estimate, and would generate the following subspaces.
Subspaces generated by a non-isotropic sparse grid with L=4, mq=3 and m=3,2.
Level Subspaces
1 1,1 
2 2,1, 1,2 
3 3,1, 2,2 
4 4,1, 3,2 
Taking this to the extreme, if the variation in the first and second dimensions are sufficiently accurately described by a level 2 quadrature rule, restricting the maximum level of both dimensions to 2 would generate the following subspaces.
Subspaces generated by a sparse grid construction with L=4, mq=2 and m=2,2.
Level Subspaces
1 1,1 
2 2,1, 1,2 
3 2,2 
4 None
Hence one subspace is generated at level 3, and no subspaces are generated at level 4. The level 3 subspace 2,2 actually indicates that this is the full grid of level 2.

3.3  Using D01ESF

D01ESF uses optional parameters, supplied in the option arrays IOPTS and OPTS. Before calling D01ESF, these option arrays must be initialized using D01ZKF. Once initialized, the required options may be set and queried using D01ZKF and D01ZLF respectively. A complete list of the options available may be found in Section 11.
You may control the maximum level required, L, using the optional parameter Maximum Level. Furthermore, you may control the first level at which the error comparison will take place using the optional parameter Minimum Level, allowing for the forced evaluation of a predetermined number of levels before the routine attempts to complete. Completion is flagged when the error estimate is sufficiently small:
F^ d k - F^ d k-1 ≤ maxεa, εr × F^ d k ,  
where εa and εr are the absolute and relative error tolerances, respectively, and k≤L is the highest level at which computation was performed. The tolerances εa and εr can be controlled by setting the optional parameters Absolute Tolerance and Relative Tolerance.
Owing to the interlacing nature of the quadrature rules used herein, abscissae x required in lower level subspaces will also appear in higher-level subspaces. This allows for calculations which will be repeated later to be stored and re-used. However, this is naturally at the expense of memory. It may also be at the expense of computational time, depending on the complexity of the integrands, as the lookup time for a given value is (at worst) Od. Furthermore, as the sparse grid level increases, fewer subsequent levels will require values from the current level. You may control the number of levels for which values are stored by setting the optional parameter Index Level.
Two different sets of interlacing quadrature rules are selectable using the optional parameter Quadrature Rule: Gauss–Patterson and Clenshaw–Curtis. Gauss–Patterson rules offer greater polynomial accuracy, whereas Clenshaw–Curtis rules are often effective for oscillatory integrands. Clenshaw–Curtis rules require function values to be evaluated on the boundary of the hypercube, whereas Gauss–Patterson rules do not. Both of these rules use precomputed weights, and as such there is an effective limit on mq; see the description of the optional parameter Quadrature Rule. The value of mq is returned by the queriable optional parameter Maximum Quadrature Level.
D01ESF also allows for non-isotropic sparse grids to be constructed. This is done by appropriately setting the array MAXDLV. It should be emphasised that a non-isometric construction should only be used if the integrands behave in a suitable way. For example, they may decay toward zero as the lesser dimensions approach their bounds of Ω. It should also be noted that setting MAXDLVk=1 will not reduce the dimension of the integrals, it will simply indicate that only one point in dimension k should be used. It is also advisable to approximate the integrals several times, once with an isometric construction of some level, and then with a non-isometric construction with higher levels in various dimensions. If the difference between the solutions is significantly more than the returned error estimates, the assumptions of dimensional importance are probably incorrect.
The abscissae in each subspace are generally expressible in a sparse manner, because many of the elements of each abscissa will in fact be the centre point of the dimension, which is termed here the ‘trivial’ element. In this routine the trivial element is always returned as 0.5 owing to the restriction to the 0,1  hypercube. As such, the subroutine F returns the abscissae in Compressed Column Storage (CCS) format (see the F11 Chapter Introduction). This has particular advantages when using accelerator hardware to evaluate the required functions, as much less data must be forwarded. It also, potentially, allows for calculations to be computed faster, as any sub-calculations dependent upon the trivial value may be potentially re-used. See the example in Section 10.

4  References

Caflisch R E, Morokoff W and Owen A B (1997) Valuation of mortgage backed securities using Brownian bridges to reduce effective dimension Journal of Computational Finance 1 27–46
Gerstner T and Griebel M (1998) Numerical integration using sparse grids Numerical Algorithms 18 209–232
Smolyak S A (1963) Quadrature and interpolation formulas for tensor products of certain classes of functions Dokl. Akad. Nauk SSSR 4 240–243

5  Parameters

1:     NI – INTEGERInput
On entry: ni, the number of integrands.
Constraint: NI≥1.
2:     NDIM – INTEGERInput
On entry: d, the number of dimensions.
Constraint: NDIM≥1.
3:     F – SUBROUTINE, supplied by the user.External Procedure
F must return the value of the integrands fj at a set of nx d-dimensional points xi, implicitly supplied as columns of a matrix Xd,nx. If X was supplied explicitly you would find that most of the elements attain the same value, xtr; the larger the number of dimensions, the greater the proportion of elements of X would be equal to xtr. So, X is effectively a sparse matrix, except that the ‘zero’ elements are replaced by elements that are all equal to the value xtr. For this reason X is supplied, as though it were a sparse matrix, in compressed column storage (CCS) format (see the F11 Chapter Introduction).
Individual entries xk,i of X, for k=1,2,…,d, are either trivially valued, in which case xk,i=xtr, or are non-trivially valued. For point i, the non-trivial row indices and corresponding abscissae values are supplied in elements ci=ICOLZPi,…,ICOLZPi+1-1, for i=1,2,…,nx, of the arrays IROWIX and XS, respectively. Hence the ith column of the matrix X is retrievable as
X IROWIX ci ,i = XS ci , X k ∉ IROWIX ci ,i = xtr .  
An equivalent integer valued matrix Q is also implicitly provided. This contains the unique indices qk,i of the underlying one-dimensional quadrature rule corresponding to the individual abscissae xk,i. For trivial abscissae, the implicit index qk,i=1. Q is supplied in the same CCS format as X, with the non-trivial values supplied in QS.
The specification of F is:
SUBROUTINE F ( NI, NDIM, NX, XTR, NNTR, ICOLZP, IROWIX, XS, QS, FM, IFLAG, IUSER, RUSER)
INTEGER  NI, NDIM, NX, NNTR, ICOLZP(NX+1), IROWIX(NNTR), QS(NNTR), IFLAG, IUSER(*)
REAL (KIND=nag_wp)  XTR, XS(NNTR), FM(NI,NX), RUSER(*)
1:     NI – INTEGERInput
On entry: ni, the number of integrands.
2:     NDIM – INTEGERInput
On entry: d, the number of dimensions.
3:     NX – INTEGERInput
On entry: nx, the number of points xi, corresponding to the number of columns of X, at which the set of integrands must be evaluated.
4:     XTR – REAL (KIND=nag_wp)Input
On entry: xtr, the value of the trivial elements of X.
5:     NNTR – INTEGERInput
On entry: if IFLAG>0, the number of non-trivial elements of X.
If IFLAG=0, the total number of abscissae from the underlying one-dimensional quadrature.
6:     ICOLZPNX+1 – INTEGER arrayInput
On entry: the set ICOLZPi,…,ICOLZPi+1 - 1 contains the indices of IROWIX and XS corresponding to the non-trivial elements of column i of X and hence of the point xi, for i=1,2,…,nx.
7:     IROWIXNNTR – INTEGER arrayInput
On entry: the row indices corresponding to the non-trivial elements of X.
8:     XSNNTR – REAL (KIND=nag_wp) arrayInput
On entry: xk,i≠xtr, the non-trivial entries of X.
9:     QSNNTR – INTEGER arrayInput
On entry: qk,i≠1, the indices of the underlying quadrature rules corresponding to xk,i≠xtr.
10:   FMNINX – REAL (KIND=nag_wp) arrayOutput
On exit: FMpi = fp xi , for i=1,2,…,nx and p=1,2,…,ni.
11:   IFLAG – INTEGERInput/Output
On entry: if IFLAG=0, this is the first call to F. nx=1, and the entire point x1 will satisfy xk,1=xtr, for k=1,2,…,d. In addition, NNTR contains the total number of abscissae from the underlying one-dimensional quadrature; XS contains the complete set of abscissae and QS contains the corresponding quadrature indices, with XS1 = xtr  and QS1 = 1 . This will always be called in serial.
In subsequent calls to F, IFLAG=1. Subsequent calls may be made from within an OpenMP parallel region. See Section 8 for details.
On exit: set IFLAG<0 if you wish to force an immediate exit from D01ESF with IFAIL=-1.
12:   IUSER* – INTEGER arrayUser Workspace
13:   RUSER* – REAL (KIND=nag_wp) arrayUser Workspace
F is called with the parameters IUSER and RUSER as supplied to D01ESF. You are free to use the arrays IUSER and RUSER to supply information to F as an alternative to using COMMON global variables.
F must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which D01ESF is called. Parameters denoted as Input must not be changed by this procedure.
4:     MAXDLVNDIM – INTEGER arrayInput
On entry: m, the array of maximum levels for each dimension. MAXDLVj, for j=1,2,…,d, contains mj, the maximum level of quadrature rule dimension j will support.
The default value, minmq,L  will be used if either MAXDLVj≤0 or MAXDLVj≥minmq,L (for details on the default values for mq and L and on how to change these values see the options Maximum Level, Maximum Quadrature Level and Quadrature Rule).
If MAXDLVj=1 for all j, only one evaluation will be performed, and as such no error estimation will be possible.
Suggested value: MAXDLVj=0 for all j=1,2,…,d.
Note: setting non-default values for some dimensions makes the assumption that the contribution from the omitted subspaces is 0. The integral and error estimates will only be based on included subspaces, which if the 0 contribution assumption is not valid will be erroneous.
5:     DINESTNI – REAL (KIND=nag_wp) arrayOutput
On exit: DINESTp  contains the final estimate F^p of the definite integral Fp, for p=1,2,…,ni.
6:     ERRESTNI – REAL (KIND=nag_wp) arrayOutput
On exit: ERRESTp  contains the final error estimate Ep of the definite integral Fp, for p=1,2,…,ni.
7:     IVALIDNI – INTEGER arrayOutput
On exit: IVALIDp indicates the final state of integral p, for p=1,2,…,ni.
IVALIDp=0
The error estimate for integral p was below the requested tolerance.
IVALIDp=1
The error estimate for integral p was below the requested tolerance. The final level used was non-isotropic.
IVALIDp=2
The error estimate for integral p was above the requested tolerance.
IVALIDp=3
The error estimate for integral p was above max0.1F^p,0.01.
IVALIDp<0
You aborted the evaluation before an error estimate could be made.
8:     IOPTS100 – INTEGER arrayCommunication Array
9:     OPTS100 – REAL (KIND=nag_wp) arrayCommunication Array
The arrays IOPTS and OPTS must not be altered between calls to any of the routines D01ESF, D01ZKF and D01ZLF.
10:   IUSER* – INTEGER arrayUser Workspace
11:   RUSER* – REAL (KIND=nag_wp) arrayUser Workspace
IUSER and RUSER are not used by D01ESF, but are passed directly to F and may be used to pass information to this routine as an alternative to using COMMON global variables.
12:   IFAIL – INTEGERInput/Output
On 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 IFAIL≠0 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 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: D01ESF may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the routine:
IFAIL=1
The requested accuracy was not achieved for at least one integral.
IFAIL=2
No accuracy was achieved for at least one integral.
IFAIL=11
On entry, NI=value.
Constraint: NI≥1.
IFAIL=21
On entry, NDIM=value.
Constraint: NDIM≥1.
IFAIL=1001
Either the option arrays IOPTS and OPTS have not been initialized for D01ESF, or they have become corrupted.
IFAIL=-1
Exit requested from F with IFLAG=value.
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

For each integral p, an error estimate Ep is returned, where,
Ep = F^ p k - F^ p k-1 ≈ F^p - Fp ,  
where k ≤ L is the highest level at which computation was performed.

8  Parallelism and Performance

8.1  Direct Threading

D01ESF is directly threaded for parallel execution. For each level, at most nt threads will evaluate the integrands over independent subspaces of the construction, and will construct a partial sum of the level's contribution. Once all subspaces from a given level have been processed, the partial sums are combined to give the total contribution of the level, which is in turn added to the total solution. For a given number of threads, the division of subspaces between the threads, and the order in which a single thread operates over its assigned subspaces, is fixed. However, the order in which all subspaces are combined will necessarily be different to the single threaded case, which may result in some discrepency in the results between parallel and serial execution.
To mitigate this discrepency, it is recommended that D01ESF be instructed to use higher-than-working precision to accumulate the actions over the subspaces. This is done by setting the option Summation Precision=HIGHER, which is the default behaviour. This has some computational cost, although this is typically negligible in comparison to the overall runtime, particularly for non-trivial integrands.
If Summation Precision=WORKING, then the accumulation will be performed using working precision, which may provide some increase in performance. Again, this is probably negligible in comparison to the overall runtime.
For some problems, typically of lower dimension, there may be insufficient work to warrant direct threading at lower levels. For example, a three-dimensional problem will require at most 3 subspaces to be evaluated at level 2, and at most 6 subspaces at level 3. Furthermore, level 2 subspaces typically contain only 2 new multidimensional abscissae, while level 3 subspaces typically contain 2 or 4 new multidimensional abscissae depending on the Quadrature Rule. If there are more threads than the number of available subspaces at a given level, or the amount of work in each subspace is outweighed by the amount of work required to generate the parallel region, parallel efficiency will be decreased. This may be mitigated to some extent by evaluating the first sl levels in serial. The value of sl may be altered using the optional parameter Serial Levels. If sl≥L, then all levels will be evaluated in serial and no direct threading will be utilized.
If you use direct threading in the manner just described, you must ensure any access to the user arrays IUSER and RUSER is done in a thread-safe manner. These are classed as OpenMP SHARED, and are passed directly to the subroutine F for every thread.

8.2  Parallelization of F

The vectorized interface also allows for parallelization inside the subroutine F by evaluating the required integrands in parallel. Provided the values returned by F match those that would be returned without parallelizing F, the final result should match the serial result, barring any discrepencies in accumulation. If you wish to parallelize F, it is advisable to set a large value for Maximum Nx, although be aware that increasing Maximum Nx will increase the memory requirement of the routine. In general, parallelization of F should not be necessary, as the higher-level parallelism over different subspaces scales well for many problems.

9  Further Comments

Not applicable.

10  Example

The example program evaluates an estimate to the set of integrals
F = ∫ Ω sin1+x ⋮ sinni+x logx d x  
where x = ∑ j=1 d jxj. It also demonstrates a simple method to safely use IUSER and RUSER as workspace for sub-calculations when running in parallel.

10.1  Program Text

Program Text (d01esfe.f90)

10.2  Program Data

None.

10.3  Program Results

Program Results (d01esfe.r)

11  Optional Parameters

Several optional parameters in D01ESF control aspects of the algorithm, methodology used, logic or output. Their values are contained in the arrays IOPTS and OPTS; these must be initialized before calling D01ESF by first calling D01ZKF with OPTSTR set to "Initialize = D01ESF".
Each optional parameter has an associated default value; to set any of them to a non-default value, or to reset any of them to the default value, use D01ZKF. The current value of an optional parameter can be queried using D01ZLF.
The remainder of this section can be skipped if you wish to use the default values for all optional parameters.
The following is a list of the optional parameters available. A full description of each optional parameter is provided in Section 11.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 represent 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.
Queriable options will return the appropriate value when queried by calling D01ZLF. They will have no effect if passed to D01ZKF.
For D01ESF the maximum length of the parameter CVALUE used by D01ZLF is 15.
Absolute TolerancerDefault =ε
r=εa, the absolute tolerance required.
Index LeveliDefault =4
The maximum level at which function values are stored for later use. Larger values use increasingly more memory, and require more time to access specific elements. Lower levels result in more repeated computation.
Constraint: i≥1.
Maximum LeveliDefault =5
i=L, the maximum number of levels to evaluate.
Constraint: 1<i≤20.
Note: the maximum allowable level in any single dimension, mq, is governed by the Quadrature Rule selected. If a value greater than mq is set, only a subset of subspaces from higher levels will be used. Should this subset be empty for a given level, computation will consider the preceding level to be the maximum level and will terminate.
Maximum NxiDefault =128
i=maxnx, the maximum number of points to evaluate in a single call to F.
Constraint: 1≤i≤16384.
Maximum Quadrature LeveliQueriable only
i=mq, the maximum level of the underlying one-dimensional quadrature rule (see Quadrature Rule).
Minimum LeveliDefault =2
The minimum number of levels which must be evaluated before an error estimate is used to determine convergence.
Constraint: i>1.
Note: if the minimum level is greater than the maximum computable level, the maximum level will be used.
Quadrature RuleaDefault =Gauss–Patterson
The underlying one-dimensional quadrature rule to be used in the construction. Open rules do not require evaluations at boundaries.
Quadrature Rule=Gauss–Patterson or GP
The interlacing Gauss–Patterson rules. Level ℓ uses 2ℓ-1 abscissae. All levels are open. These rules provide high order accuracy. mq=9.
Quadrature Rule=Clenshaw–Curtis or CC
The interlacing Clenshaw–Curtis rules. Level ℓ uses 2ℓ-1+1 abscissae. All levels above level 1 are closed. mq=12.
Relative TolerancerDefault =ε
r=εa, the relative tolerance required.
Summation PrecisionaDefault =HIGHER
Determines whether D01ESF uses working precision or higher-than-working precision to accumulate the actions over subspaces.
Summation Precision=HIGHER or H
Higher-than-working precision is used to accumulate the action over a subspace, and for the accumulation of all such actions. This is more expensive computationally, although this is probably negligible in comparison to the cost of evaluating the integrands and the overall runtime. This significantly reduces variation in the result when changing the number of threads.
Summation Precision=WORKING or W
Working precision is used to accumulate the actions over subspaces. This may provide some speedup, particularly if ni or nt is large. The results of parallel simulations will however be more prone to variation.
Note: the following option is relevant only to multithreaded implementations of the NAG Library..
Serial LevelsiDefault =1
i=sl, the number of levels to be evaluated in serial before initializing parallelization. For relatively trivial integrands, this may need to be set greater than the default to reduce parallel overhead.
D01ESF (PDF version)
D01 Chapter Contents
D01 Chapter Introduction
NAG Library Manual
© The Numerical Algorithms Group Ltd, Oxford, UK. 2015