# NAG CL Interfaced01gdc (md_​numth_​vec)

Settings help

CL Name Style:

## 1Purpose

d01gdc calculates an approximation to a definite integral in up to $20$ dimensions, using the Korobov–Conroy number theoretic method.

## 2Specification

 #include
void  d01gdc (Integer ndim,
 void (*vecfun)(Integer ndim, const double x[], double fv[], Integer m, Nag_Comm *comm),
 void (*vecreg)(Integer ndim, const double x[], Integer j, double c[], double d[], Integer m, Nag_Comm *comm),
Integer npts, double vk[], Integer nrand, Nag_Boolean transform, double *res, double *err, Nag_Comm *comm, NagError *fail)
The function may be called by the names: d01gdc or nag_quad_md_numth_vec.

## 3Description

d01gdc calculates an approximation to the integral
 $I= ∫ c1 d1 ⋯ ∫ cn dn f (x1,…,xn) dxn … dx1$ (1)
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}_{i-1}$, 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
 $xi = ci + (di-ci) yi , i= 1 , 2 ,…, n .$
The method then uses as its basis the number theoretic formula for the $n$-cube, ${\left[0,1\right]}^{n}$:
 $∫01 ⋯ ∫01 g (x1,…,xn) dxn ⋯ dx1 = 1p ∑k=1p g ({ka1p},…,{kanp}) - E$ (2)
where $\left\{x\right\}$ denotes the fractional part of $x$, ${a}_{1},\dots ,{a}_{n}$ are the so-called 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
 $E≤CK p-α ln αβ ⁡p$ (3)
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
 (4)
and gives a procedure for calculating the argument, $a$, to satisfy the optimal conditions.
In this function the periodisation is achieved by the simple transformation
 $xi = yi2 (3-2⁢yi) , i= 1 , 2 ,…, n .$
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 function provides built-in sets of optimal coefficients, corresponding to six different values of $p$. Alternatively, the optimal coefficients may be supplied by you. Functions d01gyc and d01gzc compute the optimal coefficients for the cases where $p$ is a prime number or $p$ is a product of two primes, respectively.
Conroy H (1967) Molecular Shroedinger equation VIII. A new method for evaluting multi-dimensional 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

## 5Arguments

1: $\mathbf{ndim}$Integer Input
On entry: $n$, the number of dimensions of the integral.
Constraint: $1\le {\mathbf{ndim}}\le 20$.
2: $\mathbf{vecfun}$function, supplied by the user External Function
vecfun must evaluate the integrand at a specified set of points.
The specification of vecfun is:
 void vecfun (Integer ndim, const double x[], double fv[], Integer m, Nag_Comm *comm)
1: $\mathbf{ndim}$Integer Input
On entry: $n$, the number of dimensions of the integral.
2: $\mathbf{x}\left[{\mathbf{m}}×{\mathbf{ndim}}\right]$const double Input
Note: where ${\mathbf{X}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{x}}\left[\left(j-1\right)×{\mathbf{m}}+i-1\right]$.
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[\mathit{dim}\right]$double Output
On exit: ${\mathbf{fv}}\left[\mathit{i}-1\right]$ must contain the value of the integrand of the $\mathit{i}$th point, i.e., ${\mathbf{fv}}\left[\mathit{i}-1\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}$Integer Input
On entry: the number of points $m$ at which the integrand is to be evaluated.
5: $\mathbf{comm}$Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to vecfun.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling d01gdc you may allocate memory and initialize these pointers with various quantities for use by vecfun when called from d01gdc (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: vecfun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01gdc. If your code inadvertently does return any NaNs or infinities, d01gdc is likely to produce unexpected results.
3: $\mathbf{vecreg}$function, supplied by the user External Function
vecreg must evaluate the limits of integration in any dimension for a set of points.
The specification of vecreg is:
 void vecreg (Integer ndim, const double x[], Integer j, double c[], double d[], Integer m, Nag_Comm *comm)
1: $\mathbf{ndim}$Integer Input
On entry: $n$, the number of dimensions of the integral.
2: $\mathbf{x}\left[{\mathbf{m}}×{\mathbf{ndim}}\right]$const double Input
Note: where ${\mathbf{X}}\left(i,j\right)$ appears in this document, it refers to the array element ${\mathbf{x}}\left[\left(j-1\right)×{\mathbf{m}}+i-1\right]$.
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,j-1\right)$ contain the current values of the first $\left(j-1\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}$Integer Input
On entry: the index $j$ for which the limits of the range of integration are required.
4: $\mathbf{c}\left[\mathit{dim}\right]$double Output
On exit: ${\mathbf{c}}\left[\mathit{i}-1\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[\mathit{dim}\right]$double Output
On exit: ${\mathbf{d}}\left[\mathit{i}-1\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}$Integer Input
On entry: the number of points $m$ at which the limits of integration must be specified.
7: $\mathbf{comm}$Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to vecreg.
userdouble *
iuserInteger *
pPointer
The type Pointer will be void *. Before calling d01gdc you may allocate memory and initialize these pointers with various quantities for use by vecreg when called from d01gdc (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
Note: vecreg should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d01gdc. If your code inadvertently does return any NaNs or infinities, d01gdc is likely to produce unexpected results.
4: $\mathbf{npts}$Integer Input
On entry: the Korobov rule to be used. There are two alternatives depending on the value of npts.
1. (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$.
2. (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]$double Input/Output
On entry: if ${\mathbf{npts}}>6$, vk must contain the $n$ optimal coefficients (which may be calculated using d01gyc or d01gzc).
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}$Integer Input
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 function 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}}×{\mathbf{npts}}$. If $1\le {\mathbf{npts}}\le 6$, the number of integrand evaluations will be ${\mathbf{nrand}}×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{transform}$Nag_Boolean Input
On entry: indicates whether the periodising transformation is to be used.
${\mathbf{transform}}=\mathrm{Nag_TRUE}$
The transformation is to be used.
${\mathbf{transform}}=\mathrm{Nag_FALSE}$
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{transform}}=\mathrm{Nag_TRUE}$.
8: $\mathbf{res}$double * Output
On exit: the approximation to the integral $I$.
9: $\mathbf{err}$double * Output
On exit: the standard error as computed from nrand sample values. If ${\mathbf{nrand}}=1$, err contains zero.
10: $\mathbf{comm}$Nag_Comm *
The NAG communication argument (see Section 3.1.1 in the Introduction to the NAG Library CL Interface).
11: $\mathbf{fail}$NagError * Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

## 6Error Indicators and Warnings

NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
On entry, argument $⟨\mathit{\text{value}}⟩$ had an illegal value.
NE_INT
On entry, ${\mathbf{ndim}}=⟨\mathit{\text{value}}⟩$.
Constraint: $1\le {\mathbf{ndim}}\le 20$.
On entry, ${\mathbf{npts}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{npts}}\ge 1$.
On entry, ${\mathbf{nrand}}=⟨\mathit{\text{value}}⟩$.
Constraint: ${\mathbf{nrand}}\ge 1$.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.

## 7Accuracy

If ${\mathbf{nrand}}>1$, an estimate of the absolute standard error is given by the value, on exit, of err.

## 8Parallelism and Performance

d01gdc 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 function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

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 functions, which must be supplied by you, may account for a significant part of the total computation time.
The time taken will be approximately proportional to ${\mathbf{nrand}}×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 d01gdc by calls to g05sac. Separate runs will produce identical answers.

## 10Example

This example calculates the integral
 $∫01 ∫01 ∫01 ∫01 cos(0.5+2⁢(x1+x2+x3+x4)-4) dx1 dx2 dx3 dx4 .$

### 10.1Program Text

Program Text (d01gdce.c)

### 10.2Program Data

Program Data (d01gdce.d)

### 10.3Program Results

Program Results (d01gdce.r)