NAG Toolbox: nag_quad_md_sgq_multi_vec (d01es)

Purpose

Syntax

Description

Comparing Quadrature Over Full and Sparse Grids

Sparse Grid Quadrature

Using nag_quad_md_sgq_multi_vec (d01es)

References

Parameters

Compulsory Input Parameters

1:     $\mathrm{ni}$ – int64int32nag_int scalar

$n_i$, the number of integrands.

Constraint: $\mathbf{ni}\ge 1$.

2:     $\mathrm{f}$ – function handle or string containing name of m-file

f must return the value of the integrands $f_j$ at a set of $n_x$ $d$-dimensional points ${\mathbf{x}}_i$, implicitly supplied as columns of a matrix $X\left(d,n_x\right)$. If $X$ was supplied explicitly you would find that most of the elements attain the same value, $x_{\mathrm{tr}}$; the larger the number of dimensions, the greater the proportion of elements of $X$ would be equal to $x_{\mathrm{tr}}$. So, $X$ is effectively a sparse matrix, except that the ‘zero’ elements are replaced by elements that are all equal to the value $x_{\mathrm{tr}}$. 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 $x_{\mathit{k},i}$ of $X$, for $\mathit{k}=1,2,\dots ,d$, are either trivially valued, in which case $x_{k,i}=x_{\mathrm{tr}}$, or are non-trivially valued. For point $i$, the non-trivial row indices and corresponding abscissae values are supplied in elements $c\left(\mathit{i}\right)=\mathbf{icolzp}\left(\mathit{i}\right),\dots ,\mathbf{icolzp}\left(\mathit{i}+1\right)-1$, for $\mathit{i}=1,2,\dots ,n_x$, of the arrays irowix and xs, respectively. Hence the $i$th column of the matrix $X$ is retrievable as

$$\begin{array}{l}X\left(\mathbf{irowix}\left(c\left(i\right)\right),i\right)=\mathbf{xs}\left(c\left(i\right)\right)\text{,}\\ \\ X\left(k\notin \mathbf{irowix}\left(c\left(i\right)\right),i\right)=x_{\mathrm{tr}}\text{.}\end{array}$$

An equivalent integer valued matrix $Q$ is also implicitly provided. This contains the unique indices $q_{k,i}$ of the underlying one-dimensional quadrature rule corresponding to the individual abscissae $x_{k,i}$. For trivial abscissae, the implicit index $q_{k,i}=1$. $Q$ is supplied in the same CCS format as $X$, with the non-trivial values supplied in qs.

[fm, iflag, user] = f(ni, ndim, nx, xtr, nntr, icolzp, irowix, xs, qs, iflag, user)

Input Parameters

1:     $\mathrm{ni}$ – int64int32nag_int scalar

$n_i$, the number of integrands.

2:     $\mathrm{ndim}$ – int64int32nag_int scalar

$d$, the number of dimensions.

3:     $\mathrm{nx}$ – int64int32nag_int scalar

$n_x$, the number of points $x_i$, corresponding to the number of columns of $X$, at which the set of integrands must be evaluated.

4:     $\mathrm{xtr}$ – double scalar

$x_{\mathrm{tr}}$, the value of the trivial elements of $X$.

5:     $\mathrm{nntr}$ – int64int32nag_int scalar

If $\mathbf{iflag}>0$, the number of non-trivial elements of $X$.

If $\mathbf{iflag}=0$, the total number of abscissae from the underlying one-dimensional quadrature.

6:     $\mathrm{icolzp}\left(\mathbf{nx}+1\right)$ – int64int32nag_int array

The set $\left\{\mathbf{icolzp}\left(\mathit{i}\right),\dots ,\mathbf{icolzp}\left(\mathit{i}+1\right)-1\right\}$ contains the indices of irowix and xs corresponding to the non-trivial elements of column $\mathit{i}$ of $X$ and hence of the point ${\mathbf{x}}_{\mathit{i}}$, for $\mathit{i}=1,2,\dots ,n_x$.

7:     $\mathrm{irowix}\left(\mathbf{nntr}\right)$ – int64int32nag_int array

The row indices corresponding to the non-trivial elements of $X$.

8:     $\mathrm{xs}\left(\mathbf{nntr}\right)$ – double array

$x_{k,i}\ne x_{\mathrm{tr}}$, the non-trivial entries of $X$.

9:     $\mathrm{qs}\left(\mathbf{nntr}\right)$ – int64int32nag_int array

$q_{k,i}\ne 1$, the indices of the underlying quadrature rules corresponding to $x_{k,i}\ne x_{\mathrm{tr}}$.

10:   $\mathrm{iflag}$ – int64int32nag_int scalar

If $\mathbf{iflag}=0$, this is the first call to f. $n_x=1$, and the entire point ${\mathbf{x}}_1$ will satisfy $x_{\mathit{k},1}=x_{\mathrm{tr}}$, for $\mathit{k}=1,2,\dots ,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 $\mathbf{xs}\left(1\right)=x_{\mathrm{tr}}$ and $\mathbf{qs}\left(1\right)=1$. This will always be called in serial.

In subsequent calls to f, $\mathbf{iflag}=1$. Subsequent calls may be made from within an OpenMP parallel region. See Parallelism and Performance for details.

11:   $\mathrm{user}$ – Any MATLAB object

f is called from nag_quad_md_sgq_multi_vec (d01es) with the object supplied to nag_quad_md_sgq_multi_vec (d01es).

Output Parameters

1:     $\mathrm{fm}\left(\mathbf{ni},\mathbf{nx}\right)$ – double array

$\mathbf{fm}\left(\mathit{p},\mathit{i}\right)=f_{\mathit{p}}\left({\mathbf{x}}_{\mathit{i}}\right)$, for $\mathit{i}=1,2,\dots ,n_x$ and $\mathit{p}=1,2,\dots ,n_{\mathit{i}}$.

2:     $\mathrm{iflag}$ – int64int32nag_int scalar

Set $\mathbf{iflag}<0$ if you wish to force an immediate exit from nag_quad_md_sgq_multi_vec (d01es) with $\mathbf{ifail}=-\mathbf{1}$.

3:     $\mathrm{user}$ – Any MATLAB object

3:     $\mathrm{maxdlv}\left(\mathbf{ndim}\right)$ – int64int32nag_int array

Suggested value: $\mathbf{maxdlv}\left(j\right)=0$ for all $j=1,2,\dots ,d$.

$\mathbf{m}$, the array of maximum levels for each dimension. $\mathbf{maxdlv}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,d$, contains $m_{\mathit{j}}$, the maximum level of quadrature rule dimension j will support.

The default value, $\min \left(m_q,L\right)$ will be used if either $\mathbf{maxdlv}\left(j\right)\le 0$ or $\mathbf{maxdlv}\left(j\right)\ge \min \left(m_q,L\right)$ (for details on the default values for $m_q$ and $L$ and on how to change these values see the options Maximum Level, Maximum Quadrature Level and Quadrature Rule).

If $\mathbf{maxdlv}\left(j\right)=1$ for all $j$, only one evaluation will be performed, and as such no error estimation will be possible.

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.

4:     $\mathrm{iopts}\left(100\right)$ – int64int32nag_int array

5:     $\mathrm{opts}\left(100\right)$ – double array

The arrays iopts and opts must not be altered between calls to any of the functions nag_quad_md_sgq_multi_vec (d01es), nag_quad_opt_set (d01zk) and nag_quad_opt_get (d01zl).

Optional Input Parameters

1:     $\mathrm{ndim}$ – int64int32nag_int scalar

Default: the dimension of the array maxdlv.

$d$, the number of dimensions.

Constraint: $\mathbf{ndim}\ge 1$.

2:     $\mathrm{user}$ – Any MATLAB object

user is not used by nag_quad_md_sgq_multi_vec (d01es), but is passed to f. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Output Parameters

1:     $\mathrm{dinest}\left(\mathbf{ni}\right)$ – double array

$\mathbf{dinest}\left(\mathit{p}\right)$ contains the final estimate ${\hat{F}}_{\mathit{p}}$ of the definite integral $F_p$, for $\mathit{p}=1,2,\dots ,n_i$.

2:     $\mathrm{errest}\left(\mathbf{ni}\right)$ – double array

$\mathbf{errest}\left(\mathit{p}\right)$ contains the final error estimate $E_p$ of the definite integral $F_p$, for $\mathit{p}=1,2,\dots ,n_i$.

3:     $\mathrm{ivalid}\left(\mathbf{ni}\right)$ – int64int32nag_int array

$\mathbf{ivalid}\left(\mathit{p}\right)$ indicates the final state of integral $\mathit{p}$, for $\mathit{p}=1,2,\dots ,n_i$.

$\mathbf{ivalid}\left(p\right)=0$

The error estimate for integral $p$ was below the requested tolerance.

$\mathbf{ivalid}\left(p\right)=1$

The error estimate for integral $p$ was below the requested tolerance. The final level used was non-isotropic.

$\mathbf{ivalid}\left(p\right)=2$

The error estimate for integral $p$ was above the requested tolerance.

$\mathbf{ivalid}\left(p\right)=3$

The error estimate for integral $p$ was above $\max \left(0.1\left|{\hat{F}}_p\right|,0.01\right)$.

$\mathbf{ivalid}\left(p\right)<0$

You aborted the evaluation before an error estimate could be made.

4:     $\mathrm{user}$ – Any MATLAB object

5:     $\mathrm{ifail}$ – int64int32nag_int scalar

$\mathbf{ifail}=\mathbf{0}$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

   $\mathbf{ifail}=1$

The requested accuracy was not achieved for at least one integral.

   $\mathbf{ifail}=2$

No accuracy was achieved for at least one integral.

   $\mathbf{ifail}=11$

Constraint: $\mathbf{ni}\ge 1$.

   $\mathbf{ifail}=21$

Constraint: $\mathbf{ndim}\ge 1$.

   $\mathbf{ifail}=1001$

Either the option arrays iopts and opts have not been initialized for nag_quad_md_sgq_multi_vec (d01es), or they have become corrupted.

   $\mathbf{ifail}=-1$

Exit requested from f with $\mathbf{iflag}=\_$.

   $\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

   $\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

   $\mathbf{ifail}=-999$

Dynamic memory allocation failed.

Accuracy

Further Comments

Example

Open in the MATLAB editor: d01es_example

function d01es_example


fprintf('d01es example results\n\n');

% Evaluate an estimate to the set of integrals of sin(n+|x|)log(|x|),
% with |x| = sum(j*x(j)), for n=1:10 on a hypercube of dimension 4.

ni = int64(10);
d  = int64(4);

% Initialize d01es and set options
iopts = zeros(100,1,'int64');
opts  = zeros(100,1);
[iopts, opts, ifail] = d01zk('Initialize = d01es', iopts, opts);
[iopts, opts, ifail] = d01zk('Absolute Tolerance = 0.0', iopts, opts);
[iopts, opts, ifail] = d01zk('Relative Tolerance = 1.0e-3', iopts, opts);
[iopts, opts, ifail] = d01zk('Maximum Level = 6', iopts, opts);
[iopts, opts, ifail] = d01zk('Index Level = 5', iopts, opts);

maxdlv = zeros(d,1,'int64');
[dinest, errest, ivalid, user, ifail] = d01es(ni, @f, maxdlv, iopts, opts);

fprintf('%10s%13s%13s%14s\n','integrand', 'integral','error','state');
for i = 1:ni
  fprintf('%7d%16.6f%15.2e%10d\n',i,dinest(i),errest(i),ivalid(i));
end



function [fm, iflag, user] = ...
         f(ni, ndim, nx, xtr, nntr, icolzp, irowix, xs, qs, iflag, user)

  fm = zeros(ni,nx);

  % Decompose |x| as |x_{tr}| + |(x_{non-tr} - x_{tr})|:
  % Here is |x_{tr}|:
  s_tr = xtr*double(ndim*(ndim+1))/2;

  for i = 1:nx
    i1 = icolzp(i);
    i2 = icolzp(i+1)-1;
    
    % Here is |(x_{non-tr} - x_{tr})|:
    s_ntr = sum(double(irowix(i1:i2)).*(xs(i1:i2)-xtr));
    
    for j = 1:ni
      fm(j,i) = sin(double(j)+s_ntr+s_tr)*log(s_ntr+s_tr);
    end
  end

d01es example results

 integrand     integral        error         state
      1        0.038352       2.40e-05         0
      2        0.401177       1.70e-05         0
      3        0.395161       5.66e-06         0
      4        0.025836       2.31e-05         0
      5       -0.367242       1.93e-05         0
      6       -0.422680       2.25e-06         0
      7       -0.089508       2.17e-05         0
      8        0.325958       2.12e-05         0
      9        0.441739       1.21e-06         0
     10        0.151388       1.99e-05         0

Optional Parameters

• Absolute Tolerance
• Index Level
• Maximum Level
• Maximum Nx
• Maximum Quadrature Level
• Minimum Level
• Quadrature Rule
• Relative Tolerance
• Serial Levels
• Summation Precision

Description of the Optional Parameters

• the keywords, where the minimum abbreviation of each keyword is underlined;
• a parameter value, where the letters $a$, $i$ and $r$ denote options that take character, integer and real values respectively.
• the default value.