BS2IN

Computes a two-dimensional tensor-product spline interpolant, returning the tensor-product B-spline coefficients.

Required Arguments

XDATA — Array of length NXDATA containing the data points in the X-direction. (Input)

XDATA must be strictly increasing.

XDATA must be strictly increasing.

YDATA — Array of length NYDATA containing the data points in the Y-direction. (Input)

YDATA must be strictly increasing.

YDATA must be strictly increasing.

FDATA — Array of size NXDATA by NYDATA containing the values to be interpolated. (Input)

FDATA (I, J) is the value at (XDATA (I), YDATA(J)).

FDATA (I, J) is the value at (XDATA (I), YDATA(J)).

KXORD — Order of the spline in the X-direction. (Input)

KXORD must be less than or equal to NXDATA.

KXORD must be less than or equal to NXDATA.

KYORD — Order of the spline in the Y-direction. (Input)

KYORD must be less than or equal to NYDATA.

KYORD must be less than or equal to NYDATA.

XKNOT — Array of length NXDATA + KXORD containing the knot sequence in the X-direction. (Input)

XKNOT must be nondecreasing.

XKNOT must be nondecreasing.

YKNOT — Array of length NYDATA + KYORD containing the knot sequence in the Y-direction. (Input)

YKNOT must be nondecreasing.

YKNOT must be nondecreasing.

BSCOEF — Array of length NXDATA * NYDATA containing the tensor-product B-spline coefficients. (Output)

BSCOEF is treated internally as a matrix of size NXDATA by NYDATA.

BSCOEF is treated internally as a matrix of size NXDATA by NYDATA.

Optional Arguments

NXDATA — Number of data points in the X-direction. (Input)

Default: NXDATA = size (XDATA,1).

Default: NXDATA = size (XDATA,1).

NYDATA — Number of data points in the Y-direction. (Input)

Default: NYDATA = size (YDATA,1).

Default: NYDATA = size (YDATA,1).

LDF — The leading dimension of FDATA exactly as specified in the dimension statement of the calling program. (Input)

Default: LDF = size (FDATA,1).

Default: LDF = size (FDATA,1).

FORTRAN 90 Interface

Generic: CALL BS2IN (XDATA, YDATA, FDATA, KXORD, KYORD, XKNOT, YKNOT,

BSCOEF [, …])

BSCOEF [, …])

Specific: The specific interface names are S_BS2IN and D_BS2IN.

FORTRAN 77 Interface

Single: CALL BS2IN (NXDATA, XDATA, NYDATA, YDATA, FDATA, LDF, KXORD, KYORD, XKNOT, YKNOT, BSCOEF)

Double: The double precision name is DBS2IN.

Description

The routine BS2IN computes a tensor product spline interpolant. The tensor product spline interpolant to data {(xi, yj, fij)}, where 1 ≤ i ≤ Nx and 1 ≤ j ≤ Ny, has the form

where kx and ky are the orders of the splines. (These numbers are passed to the subroutine in KXORD and KYORD, respectively.) Likewise, tx and ty are the corresponding knot sequences (XKNOT and YKNOT). The algorithm requires that

tx(kx) ≤ xi ≤ tx(Nx + 1) 1 ≤ i ≤ Nx

ty(ky) ≤ yj ≤ ty(Ny + 1) 1 ≤ j ≤ Ny

Tensor product spline interpolants in two dimensions can be computed quite efficiently by solving (repeatedly) two univariate interpolation problems. The computation is motivated by the following observations. It is necessary to solve the system of equations

Setting

we note that for each fixed i from 1 to Nx, we have Ny linear equations in the same number of unknowns as can be seen below:

The same matrix appears in all of the equations above:

Thus, we need only factor this matrix once and then apply this factorization to the Nx righthand sides. Once this is done and we have computed hmi, then we must solve for the coefficients cnm using the relation

for m from 1 to Ny, which again involves one factorization and Ny solutions to the different right-hand sides. The routine BS2IN is based on the routine SPLI2D by de Boor (1978, page 347).

Comments

1. Workspace may be explicitly provided, if desired, by use of B22IN/DB22IN. The reference is:

CALL B22IN (NXDATA, XDATA, NYDATA, YDATA, FDATA, LDF, KXORD, KYORD, XKNOT, YKNOT, BSCOEF, WK, IWK)

The additional arguments are as follows:

WK — Work array of length NXDATA * NYDATA + MAX((2 * KXORD −1) NXDATA, (2 * KYORD − 1) * NYDATA) + MAX((3 * KXORD − 2) * NXDATA, (3 * KYORD − 2) * NYDATA) + 2 * MAX(NXDATA, NYDATA).

IWK — Work array of length MAX(NXDATA, NYDATA).

2. Informational errors

Type | Code | Description |
---|---|---|

3 | 1 | Interpolation matrix is nearly singular. LU factorization failed. |

3 | 2 | Interpolation matrix is nearly singular. Iterative refinement failed. |

4 | 6 | The XDATA values must be strictly increasing. |

4 | 7 | The YDATA values must be strictly increasing. |

4 | 13 | Multiplicity of the knots cannot exceed the order of the spline. |

4 | 14 | The knots must be nondecreasing. |

4 | 15 | The I-th smallest element of the data point array must be greater than the I-th knot and less than the (I + K_ORD)-th knot. |

4 | 16 | The largest element of the data point array must be greater than the (N_DATA)-th knot and less than or equal to the (N_DATA + K_ORD)-th knot. |

4 | 17 | 7 The smallest element of the data point array must be greater than or equal to the first knot and less than the (K_ORD + 1)st knot. |

Example

In this example, a tensor product spline interpolant to a function f is computed. The values of the interpolant and the error on a 4 × 4 grid are displayed.

USE BS2IN_INT

USE BSNAK_INT

USE BS2VL_INT

USE UMACH_INT

IMPLICIT NONE

! SPECIFICATIONS FOR PARAMETERS

INTEGER KXORD, KYORD, LDF, NXDATA, NXKNOT, NXVEC, NYDATA,&

NYKNOT, NYVEC

PARAMETER (KXORD=5, KYORD=2, NXDATA=21, NXVEC=4, NYDATA=6,&

NYVEC=4, LDF=NXDATA, NXKNOT=NXDATA+KXORD,&

NYKNOT=NYDATA+KYORD)

!

INTEGER I, J, NOUT, NXCOEF, NYCOEF

REAL BSCOEF(NXDATA,NYDATA), F, FDATA(LDF,NYDATA), FLOAT,&

X, XDATA(NXDATA), XKNOT(NXKNOT), XVEC(NXVEC), Y,&

YDATA(NYDATA), YKNOT(NYKNOT), YVEC(NYVEC),VL

INTRINSIC FLOAT

! Define function

F(X,Y) = X*X*X + X*Y

! Set up interpolation points

DO 10 I=1, NXDATA

XDATA(I) = FLOAT(I-11)/10.0

10 CONTINUE

! Generate knot sequence

CALL BSNAK (NXDATA, XDATA, KXORD, XKNOT)

! Set up interpolation points

DO 20 I=1, NYDATA

YDATA(I) = FLOAT(I-1)/5.0

20 CONTINUE

! Generate knot sequence

CALL BSNAK (NYDATA, YDATA, KYORD, YKNOT)

! Generate FDATA

DO 40 I=1, NYDATA

DO 30 J=1, NXDATA

FDATA(J,I) = F(XDATA(J),YDATA(I))

30 CONTINUE

40 CONTINUE

! Interpolate

CALL BS2IN (XDATA, YDATA, FDATA, KXORD, KYORD, XKNOT, YKNOT,&

BSCOEF)

NXCOEF = NXDATA

NYCOEF = NYDATA

! Get output unit number

CALL UMACH (2, NOUT)

! Write heading

WRITE (NOUT,99999)

! Print over a grid of

! [0.0,1.0] x [0.0,1.0] at 16 points.

DO 50 I=1, NXVEC

XVEC(I) = FLOAT(I-1)/3.0

50 CONTINUE

DO 60 I=1, NYVEC

YVEC(I) = FLOAT(I-1)/3.0

60 CONTINUE

! Evaluate spline

DO 80 I=1, NXVEC

DO 70 J=1, NYVEC

VL = BS2VL (XVEC(I), YVEC(J), KXORD, KYORD, XKNOT,&

YKNOT, NXCOEF, NYCOEF, BSCOEF)

WRITE (NOUT,'(3F15.4,F15.6)') XVEC(I), YVEC(J),&

VL, (F(XVEC(I),YVEC(J))-VL)

70 CONTINUE

80 CONTINUE

99999 FORMAT (13X, 'X', 14X, 'Y', 10X, 'S(X,Y)', 9X, 'Error')

END

Output

X Y S(X,Y) Error

0.0000 0.0000 0.0000 0.000000

0.0000 0.3333 0.0000 0.000000

0.0000 0.6667 0.0000 0.000000

0.0000 1.0000 0.0000 0.000000

0.3333 0.0000 0.0370 0.000000

0.3333 0.3333 0.1481 0.000000

0.3333 0.6667 0.2593 0.000000

0.3333 1.0000 0.3704 0.000000

0.6667 0.0000 0.2963 0.000000

0.6667 0.3333 0.5185 0.000000

0.6667 0.6667 0.7407 0.000000

0.6667 1.0000 0.9630 0.000000

1.0000 0.0000 1.0000 0.000000

1.0000 0.3333 1.3333 0.000000

1.0000 0.6667 1.6667 0.000000

1.0000 1.0000 2.0000 0.000000