# NAG FL Interfacee01zaf (dimn_​grid)

## 1Purpose

e01zaf interpolates data at a point in $n$-dimensional space, that is defined by a set of gridded data points. It offers three methods to interpolate the data: Linear Interpolation, Cubic Interpolation and Weighted Average.

## 2Specification

Fortran Interface
 Subroutine e01zaf ( d, narr, axis, lx, v, k, wf, ans,
 Integer, Intent (In) :: d, narr(d), lx, method, k Integer, Intent (Inout) :: ifail Real (Kind=nag_wp), Intent (In) :: axis(lx), v(*), point(d), wf Real (Kind=nag_wp), Intent (Out) :: ans Logical, Intent (In) :: uniform
#include <nag.h>
 void e01zaf_ (const Integer *d, const Integer narr[], const logical *uniform, const double axis[], const Integer *lx, const double v[], const double point[], const Integer *method, const Integer *k, const double *wf, double *ans, Integer *ifail)
The routine may be called by the names e01zaf or nagf_interp_dimn_grid.

## 3Description

e01zaf interpolates an $n$-dimensional point within a set of gridded data points, $\mathbf{Z}=\left\{{z}_{1{j}_{1}},{z}_{2{j}_{2}},\dots ,{z}_{d{j}_{d}}\right\}$, with corresponding data values $\mathbf{F}=\left\{{f}_{1{j}_{1}},{f}_{2{j}_{2}},\dots ,{f}_{d{j}_{d}}\right\}$ where, for the $i$th dimension, ${j}_{i}=1,\dots ,{n}_{i}$ and ${n}_{i}$ is the number of ordinates in the $i$th dimension.
A hypercube of ${\left(2k\right)}^{d}$ data points $\left[{\mathbf{h}}_{1},{\mathbf{h}}_{2},\dots ,{\mathbf{h}}_{{\left(2k\right)}^{d}}\right]\subset \mathbf{Z}$, where ${\mathbf{h}}_{i}=\left({h}_{i1},{h}_{i2},\dots ,{h}_{id}\right)$ and the corresponding data values are $f\left({\mathbf{h}}_{i}\right)\in \mathbf{F}$, around the given point, $\mathbf{x}=\left({x}_{1},{x}_{2},\dots ,{x}_{d}\right)$, is found and then used to interpolate using one of the following three methods.
1. (i)Weighted Average, that is a modification of Shepard's method (Shepard (1968)) as used for scattered data in e01zmf. This method interpolates the data with the weighted mean
 $Q x = ∑ r=1 2kd wr x fr ∑ r=1 2kd wr x ,$
where ${f}_{r}=f\left({\mathbf{h}}_{r}\right)$, ${w}_{r}\left(\mathbf{x}\right)=\frac{1}{D\left(\left|\mathbf{x}-{\mathbf{h}}_{r}\right|\right)}$ and $D\left(\mathbf{y}\right)={y}_{1}^{\rho }+{y}_{2}^{\rho }+\cdots +{y}_{d}^{\rho }$, for a given value of $\rho$.
2. (ii)Linear Interpolation, which takes ${2}^{d}$ surrounding data points ($k=1$) and performs two one-dimensional linear interpolations in each dimension on data points ${\mathbf{h}}_{a}$ and ${\mathbf{h}}_{b}$, reducing the dimension every iteration until it has reached an answer. The formula for linear interpolation in dimension $i$ is simply
 $f = fa + xi-hai fb-fa hbi-hai ,$
where ${f}_{r}=f\left({\mathbf{h}}_{r}\right)$ and ${h}_{ai}<{x}_{i}<{h}_{bi}$.
3. (iii)Cubic Interpolation, based on cubic convolution (Keys (1981)). In a similar way to the Linear Interpolation method, it performs the interpolations in one dimension reducing it each time, however it requires four surrounding data points in each dimension ($k=2$), two in each direction $\left({\mathbf{h}}_{-1},{\mathbf{h}}_{0},{\mathbf{h}}_{1},{\mathbf{h}}_{2}\right)$. The following is used to calculate the one-dimensional interpolant in dimension $i$
 $f = 12 1 t t2 t3 0 2 0 0 -1 0 1 0 2 -5 4 -1 -1 3 -3 1 f-1 f0 f1 f2$
where $t={x}_{i}-{h}_{0i}$ and ${f}_{r}=f\left({\mathbf{h}}_{r}\right)$.

## 4References

Keys R (1981) Cubic Convolution Interpolation for Digital Image Processing IEEE Transactions on Acoutstics, Speech, and Signal Processing Vol ASSP-29 No. 6 1153–1160 http://hmi.stanford.edu/doc/Tech_Notes/HMI-TN-2004-004-Interpolation/Keys_cubic_interp.pdf
Shepard D (1968) A two-dimensional interpolation function for irregularly spaced data Proc. 23rd Nat. Conf. ACM 517–523 Brandon/Systems Press Inc., Princeton

## 5Arguments

1: $\mathbf{d}$Integer Input
On entry: $d$, the number of dimensions.
Constraint: ${\mathbf{d}}\ge 2$.
2: $\mathbf{narr}\left({\mathbf{d}}\right)$Integer array Input
On entry: the number of data ordinates in each dimension, with ${\mathbf{narr}}\left(\mathit{i}\right)={n}_{i}$, for $\mathit{i}=1,2,\dots ,{\mathbf{d}}$.
Constraint: ${\mathbf{narr}}\left(i\right)\ge 2$.
3: $\mathbf{uniform}$Logical Input
On entry: states whether the data points are uniformly spaced.
${\mathbf{uniform}}=\mathrm{.TRUE.}$
The data points are uniformly spaced.
${\mathbf{uniform}}=\mathrm{.FALSE.}$
The data points are not uniformly spaced.
Constraint: if ${\mathbf{method}}=3$, uniform must be .TRUE..
4: $\mathbf{axis}\left({\mathbf{lx}}\right)$Real (Kind=nag_wp) array Input
On entry: defines the axis. If the data points are uniformly spaced (see argument uniform) axis should contain the start and end of each dimension $\left({x}_{1 1},{x}_{1 {n}_{1}},\dots ,{x}_{d 1},{x}_{d {n}_{d}}\right)$.
If the data points are not uniformly spaced, axis should contain all the data ordinates for each dimension $\left({x}_{1 1},{x}_{1 2},\dots ,{x}_{1 {n}_{1}},\dots ,{x}_{d 1},{x}_{d 2},\dots ,{x}_{d {n}_{d}}\right)$.
Constraint: axis must be strictly increasing in each dimension.
5: $\mathbf{lx}$Integer Input
On entry: the dimension of the array axis as declared in the (sub)program from which e01zaf is called.
Constraints:
• if ${\mathbf{uniform}}=\mathrm{.TRUE.}$, ${\mathbf{lx}}=2{\mathbf{d}}$;
• if ${\mathbf{uniform}}=\mathrm{.FALSE.}$, ${\mathbf{lx}}={\sum }_{i=1}^{{\mathbf{d}}}{\mathbf{narr}}\left(i\right)$.
6: $\mathbf{v}\left(*\right)$Real (Kind=nag_wp) array Input
Note: the dimension of the array v must be at least $\prod _{\mathit{i}=1}^{{\mathbf{d}}}{\mathbf{narr}}\left(\mathit{i}\right)$.
On entry: holds the values of the data points in such an order that the index of a data value with coordinates $\left({z}_{1},{z}_{2},\dots ,{z}_{d}\right)$ is
 $∑ i=1 d zi∏n∈Sin,$
where ${\mathbf{S}}_{i}=\left\{{\mathbf{narr}}\left(l\right):l=1,\dots ,i-1\right\}$ e.g., $\left(\left({x}_{11},{x}_{21},\dots ,{x}_{d1}\right),\left({x}_{12},{x}_{21},\dots ,{x}_{d1}\right),\dots ,\left({x}_{1{n}_{d}},{x}_{21},\dots ,{x}_{d1}\right),\left({x}_{11},{x}_{22},\dots ,{x}_{d1}\right),\left({x}_{12},{x}_{22},\dots ,{x}_{d1}\right),\dots ,\left({x}_{1{n}_{d}},{x}_{2{n}_{d}},\dots ,{x}_{d{n}_{d}}\right)\right)$.
7: $\mathbf{point}\left({\mathbf{d}}\right)$Real (Kind=nag_wp) array Input
On entry: $\mathbf{x}$, the point at which the data value is to be interpolated.
Constraint: the point must lie inside the limits of the data values in each dimension supplied in axis.
8: $\mathbf{method}$Integer Input
On entry: the method to be used.
${\mathbf{method}}=1$
Weighted Average.
${\mathbf{method}}=2$
Linear Interpolation.
${\mathbf{method}}=3$
Cubic Interpolation.
Constraint: ${\mathbf{method}}=1$, $2$ or $3$.
9: $\mathbf{k}$Integer Input
On entry: if ${\mathbf{method}}=1$, k controls the number of data points used in the Weighted Average method, with k points used in each dimension, either side of the interpolation point. The total number of data points used for the interpolation will therefore be ${\left(2{\mathbf{k}}\right)}^{{\mathbf{d}}}$.
If ${\mathbf{method}}\ne 1$, then k is not referenced and need not be set.
Constraint: if ${\mathbf{method}}=1$, ${\mathbf{k}}\ge 1$.
10: $\mathbf{wf}$Real (Kind=nag_wp) Input
On entry: the power used for the weighted average such that a high power will cause closer points to be more heavily weighted.
Constraint: if ${\mathbf{method}}=1$, $1.0\le {\mathbf{wf}}\le 15.0$.
11: $\mathbf{ans}$Real (Kind=nag_wp) Output
On exit: holds the result of the interpolation.
12: $\mathbf{ifail}$Integer Input/Output
On entry: ifail must be set to $0$, $-1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $-1$ means that an error message is printed while a value of $1$ means that it is not.
If halting is not appropriate, the value $-1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. When the value $-\mathbf{1}$ 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).

## 6Error 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{d}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{d}}\ge 2$.
${\mathbf{ifail}}=2$
On entry, ${\mathbf{narr}}\left(〈\mathit{\text{value}}〉\right)=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{narr}}\left(\mathit{i}\right)\ge 2$.
${\mathbf{ifail}}=4$
On entry, axis decreases in dimension $〈\mathit{\text{value}}〉$.
Constraint: axis definition must be strictly increasing.
${\mathbf{ifail}}=5$
On entry, ${\mathbf{lx}}=〈\mathit{\text{value}}〉$, ${\mathbf{d}}=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{uniform}}=\mathrm{.TRUE.}$, ${\mathbf{lx}}=2{\mathbf{d}}$.
On entry, ${\mathbf{lx}}=〈\mathit{\text{value}}〉$, $\sum$narr$=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{uniform}}=\mathrm{.FALSE.}$, ${\mathbf{lx}}=\sum {\mathbf{narr}}$.
${\mathbf{ifail}}=7$
On entry, ${\mathbf{k}}=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{method}}=1$, ${\mathbf{k}}\ge 1$.
${\mathbf{ifail}}=8$
On entry, ${\mathbf{point}}\left(〈\mathit{\text{value}}〉\right)=〈\mathit{\text{value}}〉$ and data range $=\left[〈\mathit{\text{value}}〉,〈\mathit{\text{value}}〉\right]$.
Constraint: point must be within the data range.
${\mathbf{ifail}}=9$
On entry, ${\mathbf{method}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{method}}=1$, $2$ or $3$.
On entry, ${\mathbf{method}}=3$ and ${\mathbf{uniform}}=\mathrm{.FALSE.}$.
Constraint: if ${\mathbf{method}}=3$, uniform must be .TRUE..
${\mathbf{ifail}}=10$
On entry, ${\mathbf{wf}}=〈\mathit{\text{value}}〉$.
Constraint: if ${\mathbf{method}}=1$, $1.0\le {\mathbf{wf}}\le 15.0$.
${\mathbf{ifail}}=101$
Cubic Interpolation method does not have enough data surrounding point; interpolation not possible.
${\mathbf{ifail}}=201$
Warning: the size of k has been reduced, due to too few data points around point.
${\mathbf{ifail}}=-99$
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

## 7Accuracy

For most data the Cubic Interpolation method will provide the best interpolation but it is data dependent. If the data is linear, the Linear Interpolation method will be best. For noisy data the Weighted Average method is advised with ${\mathbf{wf}}<2.0$ and ${\mathbf{k}}>1$. This will include more data points and give them a greater influence to the answer.

## 8Parallelism and Performance

e01zaf is not threaded in any implementation.

None.

## 10Example

This program takes a set of uniform three-dimensional grid data points which come from the function
 $fx= x13- x22+ x3 .$
e01zaf then interpolates the data at the point $\left(1.10,0.25,0.75\right)$ using all three methods. The answers and the absolute errors are then printed.

### 10.1Program Text

Program Text (e01zafe.f90)

### 10.2Program Data

Program Data (e01zafe.d)

### 10.3Program Results

Program Results (e01zafe.r)