NAG Library Routine Document
d01gdf (md_numth_vec)
1
Purpose
d01gdf calculates an approximation to a definite integral in up to $20$ dimensions, using the Korobov–Conroy number theoretic method. This routine is designed to be particularly efficient on vector processors.
2
Specification
Fortran Interface
Integer, Intent (In)  ::  ndim, npts, nrand, itrans  Integer, Intent (Inout)  ::  ifail  Real (Kind=nag_wp), Intent (Inout)  ::  vk(ndim)  Real (Kind=nag_wp), Intent (Out)  ::  res, err  External  ::  vecfun, vecreg 

C Header Interface
#include <nagmk26.h>
void 
d01gdf_ (const Integer *ndim, void (NAG_CALL *vecfun)(const Integer *ndim, const double x[], double fv[], const Integer *m), void (NAG_CALL *vecreg)(const Integer *ndim, const double x[], const Integer *j, double c[], double d[], const Integer *m), const Integer *npts, double vk[], const Integer *nrand, const Integer *itrans, double *res, double *err, Integer *ifail) 

3
Description
d01gdf calculates an approximation to the integral
using the Korobov–Conroy number theoretic method (see
Korobov (1957),
Korobov (1963) and
Conroy (1967)). The region of integration defined in
(1) is such that generally
${c}_{i}$ and
${d}_{i}$ may be functions of
${x}_{1},{x}_{2},\dots ,{x}_{i1}$, for
$i=2,3,\dots ,n$, with
${c}_{1}$ and
${d}_{1}$ constants. The integral is first of all transformed to an integral over the
$n$cube
${\left[0,1\right]}^{n}$ by the change of variables
The method then uses as its basis the number theoretic formula for the
$n$cube,
${\left[0,1\right]}^{n}$:
where
$\left\{x\right\}$ denotes the fractional part of
$x$,
${a}_{1},\dots ,{a}_{n}$ are the socalled optimal coefficients,
$E$ is the error, and
$p$ is a prime integer. (It is strictly only necessary that
$p$ be relatively prime to all
${a}_{1},\dots ,{a}_{n}$ and is in fact chosen to be even for some cases in
Conroy (1967).) The method makes use of properties of the Fourier expansion of
$g\left({x}_{1},\dots ,{x}_{n}\right)$ which is assumed to have some degree of periodicity. Depending on the choice of
${a}_{1},\dots ,{a}_{n}$ the contributions from certain groups of Fourier coefficients are eliminated from the error,
$E$. Korobov shows that
${a}_{1},\dots ,{a}_{n}$ can be chosen so that the error satisfies
where
$\alpha $ and
$C$ are real numbers depending on the convergence rate of the Fourier series,
$\beta $ is a constant depending on
$n$, and
$K$ is a constant depending on
$\alpha $ and
$n$. There are a number of procedures for calculating these optimal coefficients. Korobov imposes the constraint that
and gives a procedure for calculating the argument,
$a$, to satisfy the optimal conditions.
In this routine the periodisation is achieved by the simple transformation
More sophisticated periodisation procedures are available but in practice the degree of periodisation does not appear to be a critical requirement of the method.
An easily calculable error estimate is not available apart from repetition with an increasing sequence of values of
$p$ which can yield erratic results. The difficulties have been studied by
Cranley and Patterson (1976) who have proposed a Monte–Carlo error estimate arising from converting
(2) into a stochastic integration rule by the inclusion of a random origin shift which leaves the form of the error
(3) unchanged; i.e., in the formula
(2),
$\left\{k\frac{{a}_{i}}{p}\right\}$ is replaced by
$\left\{{\alpha}_{i}+k\frac{{a}_{i}}{p}\right\}$, for
$i=1,2,\dots ,n$, where each
${\alpha}_{i}$, is uniformly distributed over
$\left[0,1\right]$. Computing the integral for each of a sequence of random vectors
$\alpha $ allows a ‘standard error’ to be estimated.
This routine provides builtin sets of optimal coefficients, corresponding to six different values of
$p$. Alternatively, the optimal coefficients may be supplied by you. Routines
d01gyf and
d01gzf compute the optimal coefficients for the cases where
$p$ is a prime number or
$p$ is a product of two primes, respectively.
This routine is designed to be particularly efficient on vector processors, although it is very important that you also code
vecfun and
vecreg efficiently.
4
References
Conroy H (1967) Molecular Shroedinger equation VIII. A new method for evaluting multidimensional integrals J. Chem. Phys. 47 5307–5318
Cranley R and Patterson T N L (1976) Randomisation of number theoretic methods for mulitple integration SIAM J. Numer. Anal. 13 904–914
Korobov N M (1957) The approximate calculation of multiple integrals using number theoretic methods Dokl. Acad. Nauk SSSR 115 1062–1065
Korobov N M (1963) Number Theoretic Methods in Approximate Analysis Fizmatgiz, Moscow
5
Arguments
 1: $\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
Constraint:
$1\le {\mathbf{ndim}}\le 20$.
 2: $\mathbf{vecfun}$ – Subroutine, supplied by the user.External Procedure

vecfun must evaluate the integrand at a specified set of points.
The specification of
vecfun is:
Fortran Interface
Integer, Intent (In)  ::  ndim, m  Real (Kind=nag_wp), Intent (In)  ::  x(m,ndim)  Real (Kind=nag_wp), Intent (Out)  ::  fv(m) 

C Header Interface
#include <nagmk26.h>
void 
vecfun (const Integer *ndim, const double x[], double fv[], const Integer *m) 

 1: $\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
 2: $\mathbf{x}\left({\mathbf{m}},{\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the coordinates of the $m$ points at which the integrand must be evaluated. ${\mathbf{x}}\left(i,j\right)$ contains the $j$th coordinate of the $i$th point.
 3: $\mathbf{fv}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: ${\mathbf{fv}}\left(\mathit{i}\right)$ must contain the value of the integrand of the $\mathit{i}$th point, i.e., ${\mathbf{fv}}\left(\mathit{i}\right)=f\left({\mathbf{x}}\left(\mathit{i},1\right),{\mathbf{x}}\left(\mathit{i},2\right),\dots ,{\mathbf{x}}\left(\mathit{i},{\mathbf{ndim}}\right)\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{m}}$.
 4: $\mathbf{m}$ – IntegerInput

On entry: the number of points $m$ at which the integrand is to be evaluated.
vecfun must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d01gdf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: vecfun should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d01gdf. If your code inadvertently
does return any NaNs or infinities,
d01gdf is likely to produce unexpected results.
 3: $\mathbf{vecreg}$ – Subroutine, supplied by the user.External Procedure

vecreg must evaluate the limits of integration in any dimension for a set of points.
The specification of
vecreg is:
Fortran Interface
Integer, Intent (In)  ::  ndim, j, m  Real (Kind=nag_wp), Intent (In)  ::  x(m,ndim)  Real (Kind=nag_wp), Intent (Out)  ::  c(m), d(m) 

C Header Interface
#include <nagmk26.h>
void 
vecreg (const Integer *ndim, const double x[], const Integer *j, double c[], double d[], const Integer *m) 

 1: $\mathbf{ndim}$ – IntegerInput

On entry: $n$, the number of dimensions of the integral.
 2: $\mathbf{x}\left({\mathbf{m}},{\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: for $i=1,2,\dots ,m$, ${\mathbf{x}}\left(i,1\right)$, ${\mathbf{x}}\left(i,2\right),\dots ,{\mathbf{x}}\left(i,j1\right)$ contain the current values of the first $\left(j1\right)$ coordinates of the $i$th point, which may be used if necessary in calculating the $m$ values of ${c}_{j}$ and ${d}_{j}$.
 3: $\mathbf{j}$ – IntegerInput

On entry: the index $j$ for which the limits of the range of integration are required.
 4: $\mathbf{c}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: ${\mathbf{c}}\left(\mathit{i}\right)$ must be set to the lower limit of the range for ${\mathbf{x}}\left(\mathit{i},j\right)$, for $\mathit{i}=1,2,\dots ,m$.
 5: $\mathbf{d}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: ${\mathbf{d}}\left(\mathit{i}\right)$ must be set to the upper limit of the range for ${\mathbf{x}}\left(\mathit{i},j\right)$, for $\mathit{i}=1,2,\dots ,m$.
 6: $\mathbf{m}$ – IntegerInput

On entry: the number of points $m$ at which the limits of integration must be specified.
vecreg must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which
d01gdf is called. Arguments denoted as
Input must
not be changed by this procedure.
Note: vecreg should not return floatingpoint NaN (Not a Number) or infinity values, since these are not handled by
d01gdf. If your code inadvertently
does return any NaNs or infinities,
d01gdf is likely to produce unexpected results.
 4: $\mathbf{npts}$ – IntegerInput

On entry: the Korobov rule to be used. There are two alternatives depending on the value of
npts.
(i) 
$1\le {\mathbf{npts}}\le 6$.
In this case one of six preset rules is chosen using $2129$, $5003$, $10007$, $20011$, $40009$ or $80021$ points depending on the respective value of npts being $1$, $2$, $3$, $4$, $5$ or $6$. 
(ii) 
${\mathbf{npts}}>6$.
npts is the number of actual points to be used with corresponding optimal coefficients supplied in the array vk. 
Constraint:
${\mathbf{npts}}\ge 1$.
 5: $\mathbf{vk}\left({\mathbf{ndim}}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: if
${\mathbf{npts}}>6$,
vk must contain the
$n$ optimal coefficients (which may be calculated using
d01gyf or
d01gzf).
If
${\mathbf{npts}}\le 6$,
vk need not be set.
On exit: if
${\mathbf{npts}}>6$,
vk is unchanged.
If
${\mathbf{npts}}\le 6$,
vk contains the
$n$ optimal coefficients used by the preset rule.
 6: $\mathbf{nrand}$ – IntegerInput

On entry: the number of random samples to be generated (generally a small value, say
$3$ to
$5$, is sufficient). The estimate,
res, of the value of the integral returned by the routine is then the average of
nrand calculations with different random origin shifts. If
${\mathbf{npts}}>6$, the total number of integrand evaluations will be
${\mathbf{nrand}}\times {\mathbf{npts}}$. If
$1\le {\mathbf{npts}}\le 6$, the number of integrand evaluations will be
${\mathbf{nrand}}\times p$, where
$p$ is the number of points corresponding to the six preset rules. For reasons of efficiency, these values are calculated a number at a time in
vecfun.
Constraint:
${\mathbf{nrand}}\ge 1$.
 7: $\mathbf{itrans}$ – IntegerInput

On entry: indicates whether the periodising transformation is to be used.
 ${\mathbf{itrans}}=0$
 The transformation is to be used.
 ${\mathbf{itrans}}\ne 0$
 The transformation is to be suppressed (to cover cases where the integrand may already be periodic or where you want to specify a particular transformation in the definition of vecfun).
Suggested value:
${\mathbf{itrans}}=0$.
 8: $\mathbf{res}$ – Real (Kind=nag_wp)Output

On exit: the approximation to the integral $I$.
 9: $\mathbf{err}$ – Real (Kind=nag_wp)Output

On exit: the standard error as computed from
nrand sample values. If
${\mathbf{nrand}}=1$,
err contains zero.
 10: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation 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, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry, ${\mathbf{ndim}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: $1\le {\mathbf{ndim}}\le 20$.
 ${\mathbf{ifail}}=2$

On entry, ${\mathbf{npts}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{npts}}\ge 1$.
 ${\mathbf{ifail}}=3$

On entry, ${\mathbf{nrand}}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{nrand}}\ge 1$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
If
${\mathbf{nrand}}>1$, an estimate of the absolute standard error is given by the value, on exit, of
err.
8
Parallelism and Performance
d01gdf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
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 implementationspecific information.
d01gdf performs the same computation as
d01gcf. However, the interface has been modified so that it can perform more efficiently on machines with vector processing capabilities. In particular,
vecfun and
vecreg must calculate the integrand and limits of integration at a
set of points. For some problems the amount of time spent in these two subroutines, which must be supplied by you, may account for a significant part of the total computation time.
For this reason it is vital that you consider the possibilities for vectorization in the code supplied for these two subroutines.
The time taken will be approximately proportional to
${\mathbf{nrand}}\times p$, where
$p$ is the number of points used, but may depend significantly on the efficiency of the code provided by you in
vecfun and
vecreg.
The exact values of
res and
err on return will depend (within statistical limits) on the sequence of random numbers generated within
d01gdf by calls to
g05saf. Separate runs will produce identical answers.
10
Example
This example calculates the integral
10.1
Program Text
Program Text (d01gdfe.f90)
10.2
Program Data
None.
10.3
Program Results
Program Results (d01gdfe.r)