e02bfc evaluates a cubic spline and up to its first three derivatives from its B-spline representation at a vector of points. e02bfc can be used to compute the values and derivatives of cubic spline fits and interpolants produced by reference to e01bac, e02bac and e02bec.

2 Specification

#include <nag.h>

void	e02bfc (Nag_SplineVectorSort start, Nag_Spline spline, Nag_DerivType deriv, Nag_Boolean xord, const double x[], Integer ixloc[], Integer nx, double s[], Integer pds, Integer iwrk[], Integer liwrk, NagError fail)

The function may be called by the names: e02bfc, nag_fit_dim1_spline_deriv_vector or nag_fit_1dspline_deriv_vector.

3 Description

e02bfc evaluates the cubic spline

s (x)

and optionally derivatives up to order

3

for a vector of points

x_{j}

, for

j = 1, 2, \dots, n_{x}

. It is assumed that

s (x)

is represented in terms of its B-spline coefficients

c_{i}

, for

i = 1, 2, \dots, \bar{n} + 3

, and (augmented) ordered knot set

λ_{i}

, for

i = 1, 2, \dots, \bar{n} + 7

, (see e02bac and e02bec), i.e.,

s (x) = \sum_{i = 1}^{q} c_{i} N_{i} (x) .

Here

q = \bar{n} + 3

\bar{n}

is the number of intervals of the spline and

N_{i} (x)

denotes the normalized B-spline of degree

3

(order

4

) defined upon the knots

λ_{i}, λ_{i + 1}, \dots, λ_{i + 4}

. The knots

λ_{5}, λ_{6}, \dots, λ_{\bar{n} + 3}

are the interior knots. The remaining knots,

λ_{1}

λ_{2}

λ_{3}

λ_{4}

and

λ_{\bar{n} + 4}

λ_{\bar{n} + 5}

λ_{\bar{n} + 6}

λ_{\bar{n + 7}}

are the exterior knots. The knots

λ_{4}

and

λ_{\bar{n} + 4}

are the boundaries of the spline.

Only abscissae satisfying,

λ_{4} \leq x_{j} \leq λ_{\bar{n} + 4},

will be evaluated. At a simple knot

λ_{i}

(i.e., one satisfying

λ_{i - 1} < λ_{i} < λ_{i + 1}

), the third derivative of the spline is, in general, discontinuous. At a multiple knot (i.e., two or more knots with the same value), lower derivatives, and even the spline itself, may be discontinuous. Specifically, at a point

x = u

where (exactly)

r

knots coincide (such a point is termed a knot of multiplicity

r

), the values of the derivatives of order

4 - j

, for

j = 1, 2, \dots, r

, are, in general, discontinuous. (Here

1 \leq r \leq 4

;

r > 4

is not meaningful.) The maximum order of the derivatives to be evaluated

D_{ord}

, and the left- or right-handedness of the computation when an abscissa corresponds exactly to an interior knot, are determined by the value of deriv.

Each abscissa (point at which the spline is to be evaluated)

x_{j}

contained in x has an associated enclosing interval number,

{ixloc}_{j}

either supplied or returned in ixloc (see argument start). A simple call to e02bfc would set

start = Nag_SplineVectorSort_Sorted

and the contents of ixloc need never be set nor referenced, and the following description on modes of operation can be ignored. However, where efficiency is an important consideration, the following description will help to choose the appropriate mode of operation.

The interval numbers are used to determine which B-splines must be evaluated for a given abscissa, and are defined as

{ixloc}_{j} = (\begin{array}{l} \leq 0 & x_{j} < λ_{1} \\ 4 & λ_{4} = x_{j} \\ k & λ_{k} < x_{j} < λ_{k + 1} \\ k & λ_{4} < λ_{k} = x_{j} & left derivatives \\ k & x_{j} = λ_{k + 1} < λ_{\bar{n} + 4} & right derivatives or no derivatives \\ \bar{n} + 4 & λ_{\bar{n} + 4} = x_{j} \\ > \bar{n} + 7 & x_{j} > λ_{\bar{n} + 7} \end{array})

(1)

The algorithm has two modes of vectorization, termed here sorted and unsorted, which are selectable by the argument start.

Furthermore, if the supplied abscissae are sufficiently ordered, as indicated by the argument xord, the algorithm will take advantage of significantly faster methods for the determination of both the interval numbers and the subsequent spline evaluations.

The sorted mode has two phases, a sorting phase and an evaluation phase. This mode is recommended if there are many abscissae to evaluate relative to the number of intervals of the spline, or the abscissae are distributed relatively densely over a subsection of the spline. In the first phase,

{ixloc}_{j}

is determined for each

x_{j}

and a permutation is calculated to sort the

x_{j}

by interval number. The first phase may be either partially or completely by-passed using the argument start if the enclosing segments and/or the subsequent ordering are already known a priori, for example if multiple spline coefficients

spline \to c

are to be evaluated over the same set of knots

spline \to lamda

In the second phase of the sorted mode, spline approximations are evaluated by segment, so that non-abscissa dependent calculations over a segment may be reused in the evaluation for all abscissae belonging to a specific segment. For example, all third derivatives of all abscissae in the same segment will be identical.

In the unsorted mode of vectorization, no a priori segment sorting is performed, and if the abscissae are not sufficiently ordered, the evaluation at an abscissa will be independent of evaluations at other abscissae; also non-abscissa dependent calculations over a segment will be repeated for each abscissa in a segment. This may be quicker if the number of abscissa is small in comparison to the number of knots in the spline, and they are distributed sparsely throughout the domain of the spline. This is effectively a direct vectorization of e02bbc and e02bcc, although if the enclosing interval numbers

{ixloc}_{j}

are known, these may again be provided.

If the abscissae are sufficiently ordered, then once the first abscissa in a segment is known, an efficient algorithm will be used to determine the location of the final abscissa in this segment. The spline will subsequently be evaluated in a vectorized manner for all the abscissae indexed between the first and last of the current segment.

If no derivatives are required, the spline evaluation is calculated by taking convex combinations due to de Boor (1972). Otherwise, the calculation of

s (x)

and its derivatives is based upon,

(i)evaluating the nonzero B-splines of orders $1$ , $2$ , $3$ and $4$ by recurrence (see Cox (1972) and Cox (1978)),
(ii)computing all derivatives of the B-splines of order $4$ by applying a second recurrence to these computed B-spline values (see de Boor (1972)),
(iii)multiplying the fourth-order B-spline values and their derivative by the appropriate B-spline coefficients, and summing, to yield the values of $s (x)$ and its derivatives.

The method of convex combinations is significantly faster than the recurrence based method. If higher derivatives of order

2

3

are not required, as much computation as possible is avoided.

4 References

Cox M G (1972) The numerical evaluation of B-splines J. Inst. Math. Appl. 10 134–149

Cox M G (1978) The numerical evaluation of a spline from its B-spline representation J. Inst. Math. Appl. 21 135–143

de Boor C (1972) On calculating with B-splines J. Approx. Theory 6 50–62

5 Arguments

1: $start$ – Nag_SplineVectorSort Input

On entry: indicates the completion state of the first phase of the algorithm.

$start = Nag_SplineVectorSort_Sorted$: The enclosing interval numbers ${ixloc}_{j}$ for the abscissae $x_{j}$ contained in x have not been determined, and you wish to use the sorted mode of vectorization.
$start = Nag_SplineVectorSort_Sorted_Indexed$: The enclosing interval numbers ${ixloc}_{j}$ have been determined and are provided in ixloc, however the required permutation and interval related information has not been determined and you wish to use the sorted mode of vectorization.
$start = Nag_SplineVectorSort_Sorted_Indexed_Perm$: You wish to use the sorted mode of vectorization, and the entire first phase has been completed, with the enclosing interval numbers supplied in ixloc, and the required permutation and interval related information provided in iwrk (from a previous call to e02bfc).
$start = Nag_SplineVectorSort_Unsorted$: The enclosing interval numbers ${ixloc}_{j}$ for the abscissae $x_{j}$ contained in x have not been determined, and you wish to use the unsorted mode of vectorization.
$start = Nag_SplineVectorSort_Unsorted_Indexed$: The enclosing interval numbers ${ixloc}_{j}$ for the abscissae $x_{j}$ contained in x have been supplied in ixloc, and you wish to use the unsorted mode of vectorization.

Constraint:

start = Nag_SplineVectorSort_Sorted

Nag_SplineVectorSort_Sorted_Indexed

Nag_SplineVectorSort_Sorted_Indexed_Perm

Nag_SplineVectorSort_Unsorted

Nag_SplineVectorSort_Unsorted_Indexed

Additional:

start = Nag_SplineVectorSort_Sorted

Nag_SplineVectorSort_Unsorted

should be used unless you are sure that the knot set is unchanged between calls.

2: $spline$ – Nag_Spline *

Pointer to structure of type Nag_Spline with the following members:

n – IntegerInput: On entry: $\bar{n} + 7$ , where $\bar{n}$ is the number of intervals of the spline (which is one greater than the number of interior knots, i.e., the knots strictly within the range $λ_{4}$ to $λ_{\bar{n} + 4}$ over which the spline is defined).

Constraint: $spline \to n \geq 8$ .

lamda – double *Input: On entry: a pointer to which memory of size $spline \to n$ must be allocated. $spline \to lamda [k - 1]$ must be set to the value of the $k$ th member of the complete set of knots, $λ_{k}$ , for $k = 1, 2, \dots, \bar{n} + 7$ .

Constraint: the $λ_{k}$ must be in nondecreasing order with $spline \to lamda [spline \to n - 4] > spline \to lamda [3]$ .

c – double *Input: On entry: a pointer to which memory of size $spline \to n - 4$ must be allocated. $spline \to c$ holds the coefficient $c_{i}$ of the B-spline $N_{i} (x)$ , for $i = 1, 2, \dots, \bar{n} + 3$ .

Under normal usage, the call to function e02bfc will follow at least one call to e01bac, e02bac or e02bec). In that case, the structure spline will have been set up correctly for input to e02bfc. If multiple sets of B-spline co-efficients are required for the same set of knots

λ

and the same set of abscissae

x

, multiple calls to e02bfc may be made with

spline \to c

pointing to different coefficient sets, with start set appropriately for efficiency.

3: $deriv$ – Nag_DerivType Input

On entry: determines the maximum order of derivatives required,

D_{ord}

, as well as the computational behaviour when absicssae correspond exactly to interior knots.

For abscissae satisfying

x_{j} = λ_{4}

x_{j} = λ_{\bar{n} + 4}

only right-handed or left-handed computation will be used respectively. For abscissae which do not coincide exactly with a knot, the handedness of the computation is immaterial.

$deriv = Nag_NoDerivs$: No derivatives required. $D_{ord} = 0$ . Only right-handed computation will be used at interior knots.
$deriv = Nag_LeftDerivs_1$ or $Nag_RightDerivs_1$: Only $s (x)$ and its first derivative are required. $D_{ord} = 1$ .
$deriv = Nag_LeftDerivs_2$ or $Nag_RightDerivs_2$: Only $s (x)$ and its first and second derivatives are required. $D_{ord} = 2$ .
$deriv = Nag_LeftDerivs_3$ or $Nag_RightDerivs_3$: $s (x)$ and its first, second and third derivatives are required. $D_{ord} = 3$ .

Constraint:

deriv = Nag_NoDerivs

Nag_LeftDerivs_1

Nag_RightDerivs_1

Nag_LeftDerivs_2

Nag_RightDerivs_2

Nag_LeftDerivs_3

Nag_RightDerivs_3

Additional: if left-handed computation of the spline

s

is required, a value of deriv must be chosen which computes at least the first derivative in a left-handed manner. As mentioned in Section 3, the handedness of the computation of

s

will only have an effect if at least

4

interior knots are identical.

4: $xord$ – Nag_Boolean Input

On entry: indicates whether x is supplied in a sufficiently ordered manner. If x is sufficiently ordered e02bfc will complete faster.

$xord = Nag_TRUE$: The abscissae in x are ordered at least by ascending interval, in that any two abscissae contained in the same interval are only separated by abscissae in the same interval. For example, $x_{j} < x_{j + 1}$ , for $j = 1, 2, \dots, nx - 1$ .
$xord = Nag_FALSE$: The abscissae in x are not sufficiently ordered.

5: $x [nx]$ – const double Input

On entry: the abscissae

x_{j}

, for

j = 1, 2, \dots, n_{x}

. If

start = Nag_SplineVectorSort_Sorted

Nag_SplineVectorSort_Unsorted

then evaluations will only be performed for these

x_{j}

satisfying

λ_{4} \leq x_{j} \leq λ_{\bar{n} + 4}

. Otherwise evaluation will be performed unless the corresponding element of ixloc contains an invalid interval number. Please note that if the

ixloc [j]

is a valid interval number then no check is made that

x [j]

actually lies in that interval.

Constraint: at least one abscissa must fall between

spline \to lamda [3]

and

spline \to lamda [spline \to n - 4]

6: $ixloc [nx]$ – Integer Input/Output

On entry: if

start = Nag_SplineVectorSort_Sorted_Indexed

Nag_SplineVectorSort_Sorted_Indexed_Perm

Nag_SplineVectorSort_Unsorted_Indexed

, if you wish

x_{j}

to be evaluated,

ixloc [j - 1]

must be the enclosing interval number

{ixloc}_{j}

of the abscissae

x_{j}

(see (1)). If you do not wish

x_{j}

to be evaluated, you may set the interval number to be either less than

4

or greater than

\bar{n} + 4

Otherwise, ixloc need not be set.

On exit: if

start = Nag_SplineVectorSort_Sorted_Indexed

Nag_SplineVectorSort_Sorted_Indexed_Perm

Nag_SplineVectorSort_Unsorted_Indexed

, ixloc is unchanged on exit.

Otherwise,

ixloc [j - 1]

, contains the enclosing interval number

{ixloc}_{j}

, for the abscissa supplied in

x [j - 1]

, for

j = 1, 2, \dots, n_{x}

. Evaluations will only be performed for abscissae

x_{j}

satisfying

λ_{4} \leq x_{j} \leq λ_{\bar{n} + 4}

. If evaluation is not performed

ixloc [j - 1]

is set to

0

x_{j} < λ_{4}

\bar{n} + 7

x_{j} > λ_{\bar{n} + 4}

Constraint: if

start = Nag_SplineVectorSort_Sorted_Indexed

Nag_SplineVectorSort_Sorted_Indexed_Perm

Nag_SplineVectorSort_Unsorted_Indexed

, at least one element of ixloc must be between

4

and

spline \to n - 3

7: $nx$ – Integer Input

On entry:

n_{x}

, the total number of abscissae contained in x, including any that will not be evaluated.

Constraint:

nx \geq 1

8: $s [\dim]$ – double Output

Note: the dimension, dim, of the array s must be at least

pds \times (D_{ord} + 1)

, see deriv for the definition of

D_{ord}

On exit: if

x_{j}

is valid,

S (j, d)

will contain the (

d - 1

)th derivative of

s (x)

, for

d = 1, 2, \dots, D_{ord} + 1

and

j = 1, 2, \dots, n_{x}

. In particular,

S (j, 1)

will contain the approximation of

s (x_{j})

for all legal values in x.

9: $pds$ – Integer Input

On entry: the stride separating row elements in the two-dimensional data stored in the array s.

Constraint:

pds \geq nx

, regardless of the acceptability of the elements of x.

10: $iwrk [liwrk]$ – Integer Input/Output

On entry: if

start = Nag_SplineVectorSort_Sorted_Indexed_Perm

, iwrk must be unchanged from a previous call to e02bfc with

start = Nag_SplineVectorSort_Sorted

Nag_SplineVectorSort_Sorted_Indexed

Otherwise, iwrk need not be set. Furthermore, iwrk may be NULL if

start = Nag_SplineVectorSort_Unsorted

Nag_SplineVectorSort_Unsorted_Indexed

On exit: if

start = Nag_SplineVectorSort_Unsorted

Nag_SplineVectorSort_Unsorted_Indexed

, iwrk is unchanged on exit.

Otherwise, iwrk contains the required permutation of elements of x, if any, and information related to the division of the abscissae

x_{j}

between the intervals derived from

spline \to lamda

11: $liwrk$ – Integer Input

On entry: the dimension of the array iwrk.

Constraint: if

start = Nag_SplineVectorSort_Sorted

Nag_SplineVectorSort_Sorted_Indexed

Nag_SplineVectorSort_Sorted_Indexed_Perm

liwrk \geq 3 + 3 \times nx

12: $fail$ – NagError * Input/Output

The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error Indicators and Warnings

NE_ABSCI_OUTSIDE_KNOT_INTVL: On entry, all elements of x had enclosing interval numbers in ixloc outside the domain allowed by the provided spline.
$⟨ value ⟩$ entries of x were indexed below the lower bound $⟨ value ⟩$ .
$⟨ value ⟩$ entries of x were indexed above the upper bound $⟨ value ⟩$ .
NE_ALLOC_FAIL: Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM: On entry, argument $⟨ value ⟩$ had an illegal value.
NE_INT: On entry, $nx = ⟨ value ⟩$ .
Constraint: $nx \geq 1$ .

On entry, $spline \to n = ⟨ value ⟩$ .
Constraint: $spline \to n \geq 8$ .
NE_INT_2: On entry, $liwrk = ⟨ value ⟩$ .
Constraint: $liwrk \geq 3 \times nx + 3 = ⟨ value ⟩$ .

On entry, $pds = ⟨ value ⟩$ .
Constraint: $pds \geq nx = ⟨ value ⟩$ .
NE_INT_CHANGED: On entry, $start = Nag_SplineVectorSort_Sorted_Indexed_Perm$ and nx is not consistent with the previous call to e02bfc.
On entry, $nx = ⟨ value ⟩$ .
Constraint: $nx = ⟨ value ⟩$ .
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.
NE_SPLINE_RANGE_INVALID: On entry, $spline \to lamda [3] = ⟨ value ⟩$ , $spline \to n = ⟨ value ⟩$ and $spline \to lamda [spline \to n - 4] = ⟨ value ⟩$ .
Constraint: $spline \to lamda [3] < spline \to lamda [spline \to n - 4]$ .
NW_SOME_SOLUTIONS: On entry, at least one element of x has an enclosing interval number in ixloc outside the set allowed by the provided spline. The spline has been evaluated for all x with enclosing interval numbers inside the allowable set.
$⟨ value ⟩$ entries of x were indexed below the lower bound $⟨ value ⟩$ .
$⟨ value ⟩$ entries of x were indexed above the upper bound $⟨ value ⟩$ .

7 Accuracy

The computed value of

s (x)

has negligible error in most practical situations. Specifically, this value has an absolute error bounded in modulus by

18 \times cmax \times machine precision

, where

cmax

is the largest in modulus of

c_{j}

c_{j} + 1

c_{j} + 2

and

c_{j} + 3

, and

j

is an integer such that

λ_{j} + 3 < x \leq λ_{j} + 4

. If

c_{j}

c_{j} + 1

c_{j} + 2

and

c_{j} + 3

are all of the same sign, then the computed value of

s (x)

has relative error bounded by

20 \times machine precision

. For full details see Cox (1978).

No complete error analysis is available for the computation of the derivatives of

s (x)

. However, for most practical purposes the absolute errors in the computed derivatives should be small. Note that this is in comparison to the derivatives of the spline, which may or may not be comparable to the derivatives of the function that has been approximated by the spline.

8 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

e02bfc 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.

9 Further Comments

If using the sorted mode of vectorization, the time required for the first phase to determine the enclosing intervals is approximately proportional to

O (n_{x} \log (\bar{n}))

. The time required to then generate the required permutations and interval information is

O (n_{x})

if x is ordered sufficiently, or at worst

O (n_{x} \min (n_{x}, \bar{n}) \log (\min (n_{x}, \bar{n})))

if x is not ordered. The time required by the second phase is then proportional to

O (n_{x})

If using the unsorted mode of vectorization, the time required is proportional to

O (n_{x} \log (\bar{n}))

if the enclosing interval numbers are not provided, or

O (n_{x})

if they are provided. However, the repeated calculation of various quantities will typically make this slower than the sorted mode when the ratio of abscissae to knots is high, or the abscissae are densely distributed over a relatively small subset of the intervals of the spline.

Note: the function does not test all the conditions on the knots given in the description of

spline \to lamda

in Section 5, since to do this would result in a computation time with a linear dependency upon

\bar{n}

instead of

\log (\bar{n})

. All the conditions are tested in e02bac and e02bec, however.

10 Example

This example fits a spline through a set of data points using e02bec and then evaluates the spline at a set of supplied abscissae.

e02bf: FL CL CPP AD PY MB

NAG CL Interfacee02bfc (dim1_​spline_​deriv_​vector)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG CL Interface
e02bfc (dim1_spline_deriv_vector)