BS3IN
Computes a three-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 increasing.
XDATA must be increasing.
YDATA — Array of length NYDATA containing the data points in the y-direction. (Input)
YDATA must be increasing.
YDATA must be increasing.
ZDATA — Array of length NZDATA containing the data points in the z-direction. (Input)
ZDATA must be increasing.
ZDATA must be increasing.
FDATA — Array of size NXDATA by NYDATA by NZDATA containing the values to be interpolated. (Input)
FDATA (I, J, K) contains the value at (XDATA (I), YDATA(J), ZDATA(K)).
FDATA (I, J, K) contains the value at (XDATA (I), YDATA(J), ZDATA(K)).
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.
KZORD — Order of the spline in the z-direction. (Input)
KZORD must be less than or equal to NZDATA.
KZORD must be less than or equal to NZDATA.
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.
ZKNOT — Array of length NZDATA + KZORD containing the knot sequence in the z-direction. (Input)
ZKNOT must be nondecreasing.
ZKNOT must be nondecreasing.
BSCOEF — Array of length NXDATA * NYDATA * NZDATA containing the tensor-product B-spline coefficients. (Output)
BSCOEF is treated internally as a matrix of size NXDATA by NYDATA by NZDATA.
BSCOEF is treated internally as a matrix of size NXDATA by NYDATA by NZDATA.
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).
NZDATA — Number of data points in the z-direction. (Input)
Default: NZDATA = size (ZDATA,1).
Default: NZDATA = size (ZDATA,1).
LDF — 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).
MDF — Middle dimension of FDATA exactly as specified in the dimension statement of the calling program. (Input)
Default: MDF = size (FDATA,2).
Default: MDF = size (FDATA,2).
FORTRAN 90 Interface
Generic: CALL BS3IN (XDATA, YDATA, ZDATA, FDATA, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF [, …])
Specific: The specific interface names are S_BS3IN and D_BS3IN.
FORTRAN 77 Interface
Single: CALL BS3IN (NXDATA, XDATA, NYDATA, YDATA, NZDATA, ZDATA, FDATA, LDF, MDF, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF)
Double: The double precision name is DBS3IN.
Description
The routine BS3IN computes a tensor-product spline interpolant. The tensor-product spline interpolant to data {(xi, yj, zk, fijk)}, where 1 ≤ i ≤ Nx, 1 ≤ j ≤ Ny, and 1 ≤ k ≤ Nz has the form
where kx, ky, and kz are the orders of the splines (these numbers are passed to the subroutine in KXORD, KYORD, and KZORD, respectively). Likewise, tx, ty, and tz are the corresponding knot sequences (XKNOT, YKNOT, and ZKNOT). The algorithm requires that
Tensor-product spline interpolants can be computed quite efficiently by solving (repeatedly) three 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 pair ij we have Nz linear equations in the same number of unknowns as can be seen below:
The same interpolation matrix appears in all of the equations above:
Thus, we need only factor this matrix once and then apply it to the NxNy right-hand sides. Once this is done and we have computed hlij, then we must solve for the coefficients cnml using the relation
that is the bivariate tensor-product problem addressed by the IMSL routine BS2IN. The interested reader should consult the algorithm description in the two-dimensional routine if more detail is desired. The routine BS3IN is based on the routine SPLI2D by de Boor (1978, page 347).
Comments
1. Workspace may be explicitly provided, if desired, by use of B23IN/DB23IN. The reference is:
CALL B23IN (NXDATA, XDATA, NYDATA, YDATA, NZDAYA, ZDATA, FDATA, LDF, MDF, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF, WK, IWK)
The additional arguments are as follows:
WK — Work array of length
MAX((2 * KXORD − 1) * NXDATA, (2 * KYORD − 1) * NYDATA,
(2 * KZORD − 1) * NZDATA) +
MAX((3 * KXORD − 2) * NXDATA, (3 * KYORD − 2) * NYDATA,
(3 * KZORD − 2) * NZDATA)
+ NXDATA * NYDATA *NZDATA
+ 2 * MAX(NXDATA, NYDATA, NZDATA).
MAX((2 * KXORD − 1) * NXDATA, (2 * KYORD − 1) * NYDATA,
(2 * KZORD − 1) * NZDATA) +
MAX((3 * KXORD − 2) * NXDATA, (3 * KYORD − 2) * NYDATA,
(3 * KZORD − 2) * NZDATA)
+ NXDATA * NYDATA *NZDATA
+ 2 * MAX(NXDATA, NYDATA, NZDATA).
IWK — Work array of length MAX(NXDATA, NYDATA, NZDATA).
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 | 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 Ith 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 | 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. |
4 | 18 | The XDATA values must be strictly increasing. |
4 | 19 | The YDATA values must be strictly increasing. |
4 | 20 | The ZDATA values must be strictly increasing. |
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 × 2 grid are displayed.
USE BS3IN_INT
USE BSNAK_INT
USE UMACH_INT
USE BS3GD_INT
IMPLICIT NONE
! SPECIFICATIONS FOR PARAMETERS
INTEGER KXORD, KYORD, KZORD, LDF, MDF, NXDATA, NXKNOT, NXVEC,&
NYDATA, NYKNOT, NYVEC, NZDATA, NZKNOT, NZVEC
PARAMETER (KXORD=5, KYORD=2, KZORD=3, NXDATA=21, NXVEC=4,&
NYDATA=6, NYVEC=4, NZDATA=8, NZVEC=2, LDF=NXDATA,&
MDF=NYDATA, NXKNOT=NXDATA+KXORD, NYKNOT=NYDATA+KYORD,&
NZKNOT=NZDATA+KZORD)
!
INTEGER I, J, K, NOUT, NXCOEF, NYCOEF, NZCOEF
REAL BSCOEF(NXDATA,NYDATA,NZDATA), F,&
FDATA(LDF,MDF,NZDATA), FLOAT, VALUE(NXVEC,NYVEC,NZVEC)&
, X, XDATA(NXDATA), XKNOT(NXKNOT), XVEC(NXVEC), Y,&
YDATA(NYDATA), YKNOT(NYKNOT), YVEC(NYVEC), Z,&
ZDATA(NZDATA), ZKNOT(NZKNOT), ZVEC(NZVEC)
INTRINSIC FLOAT
! Define function.
F(X,Y,Z) = X*X*X + X*Y*Z
! Set up X-interpolation points
DO 10 I=1, NXDATA
XDATA(I) = FLOAT(I-11)/10.0
10 CONTINUE
! Set up Y-interpolation points
DO 20 I=1, NYDATA
YDATA(I) = FLOAT(I-1)/FLOAT(NYDATA-1)
20 CONTINUE
! Set up Z-interpolation points
DO 30 I=1, NZDATA
ZDATA(I) = FLOAT(I-1)/FLOAT(NZDATA-1)
30 CONTINUE
! Generate knots
CALL BSNAK (NXDATA, XDATA, KXORD, XKNOT)
CALL BSNAK (NYDATA, YDATA, KYORD, YKNOT)
CALL BSNAK (NZDATA, ZDATA, KZORD, ZKNOT)
! Generate FDATA
DO 50 K=1, NZDATA
DO 40 I=1, NYDATA
DO 40 J=1, NXDATA
FDATA(J,I,K) = F(XDATA(J),YDATA(I),ZDATA(K))
40 CONTINUE
50 CONTINUE
! Get output unit number
CALL UMACH (2, NOUT)
! Interpolate
CALL BS3IN (XDATA, YDATA, ZDATA, FDATA, KXORD, &
KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF)
!
NXCOEF = NXDATA
NYCOEF = NYDATA
NZCOEF = NZDATA
! Write heading
WRITE (NOUT,99999)
! Print over a grid of
! [-1.0,1.0] x [0.0,1.0] x [0.0,1.0]
! at 32 points.
DO 60 I=1, NXVEC
XVEC(I) = 2.0*(FLOAT(I-1)/3.0) - 1.0
60 CONTINUE
DO 70 I=1, NYVEC
YVEC(I) = FLOAT(I-1)/3.0
70 CONTINUE
DO 80 I=1, NZVEC
ZVEC(I) = FLOAT(I-1)
80 CONTINUE
! Call the evaluation routine.
CALL BS3GD (0, 0, 0, XVEC, YVEC, ZVEC,&
KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF, VALUE)
DO 110 I=1, NXVEC
DO 100 J=1, NYVEC
DO 90 K=1, NZVEC
WRITE (NOUT,'(4F13.4, F13.6)') XVEC(I), YVEC(K),&
ZVEC(K), VALUE(I,J,K),&
F(XVEC(I),YVEC(J),ZVEC(K))&
- VALUE(I,J,K)
90 CONTINUE
100 CONTINUE
110 CONTINUE
99999 FORMAT (10X, 'X', 11X, 'Y', 10X, 'Z', 10X, 'S(X,Y,Z)', 7X,&
'Error')
END
Output
X Y Z S(X,Y,Z) Error
-1.0000 0.0000 0.0000 -1.0000 0.000000
-1.0000 0.3333 1.0000 -1.0000 0.000000
-1.0000 0.0000 0.0000 -1.0000 0.000000
-1.0000 0.3333 1.0000 -1.3333 0.000000
-1.0000 0.0000 0.0000 -1.0000 0.000000
-1.0000 0.3333 1.0000 -1.6667 0.000000
-1.0000 0.0000 0.0000 -1.0000 0.000000
-1.0000 0.3333 1.0000 -2.0000 0.000000
-0.3333 0.0000 0.0000 -0.0370 0.000000
-0.3333 0.3333 1.0000 -0.0370 0.000000
-0.3333 0.0000 0.0000 -0.0370 0.000000
-0.3333 0.3333 1.0000 -0.1481 0.000000
-0.3333 0.0000 0.0000 -0.0370 0.000000
-0.3333 0.3333 1.0000 -0.2593 0.000000
-0.3333 0.0000 0.0000 -0.0370 0.000000
-0.3333 0.3333 1.0000 -0.3704 0.000000
0.3333 0.0000 0.0000 0.0370 0.000000
0.3333 0.3333 1.0000 0.0370 0.000000
0.3333 0.0000 0.0000 0.0370 0.000000
0.3333 0.3333 1.0000 0.1481 0.000000
0.3333 0.0000 0.0000 0.0370 0.000000
0.3333 0.3333 1.0000 0.2593 0.000000
0.3333 0.0000 0.0000 0.0370 0.000000
0.3333 0.3333 1.0000 0.3704 0.000000
1.0000 0.0000 0.0000 1.0000 0.000000
1.0000 0.3333 1.0000 1.0000 0.000000
1.0000 0.0000 0.0000 1.0000 0.000000
1.0000 0.3333 1.0000 1.3333 0.000000
1.0000 0.0000 0.0000 1.0000 0.000000
1.0000 0.3333 1.0000 1.6667 0.000000
1.0000 0.0000 0.0000 1.0000 0.000000
1.0000 0.3333 1.0000 2.0000 0.000000
