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.
YDATA — Array of length NYDATA containing the data points in the Y-direction. (Input)
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)).
KXORD — Order of the spline in the X-direction. (Input)
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.
XKNOT — Array of length NXDATA + KXORD containing the knot sequence in the X-direction. (Input)
XKNOT must be nondecreasing.
YKNOT — Array of length NYDATA + KYORD containing the knot sequence in the Y-direction. (Input)
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.
Optional Arguments
NXDATA — Number of data points in the X-direction. (Input)
Default: NXDATA = size (XDATA,1).
NYDATA — Number of data points in the Y-direction. (Input)
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).
FORTRAN 90 Interface
Generic: CALL BS2IN (XDATA, YDATA, FDATA, KXORD, KYORD, XKNOT, YKNOT,
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
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
