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.
YDATA — Array of
length NYDATA
containing the data points in the y-direction. (Input)
YDATA must
be increasing.
ZDATA — Array of
length NZDATA
containing the data points in the z-direction. (Input)
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)).
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.
KZORD — Order of
the spline in the z-direction. (Input)
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.
YKNOT — Array of
length NYDATA +
KYORD containing
the knot sequence in the y-direction. (Input)
YKNOT must be
nondecreasing.
ZKNOT — Array of
length NZDATA +
KZORD containing
the knot sequence in the z-direction. (Input)
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.
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).
NZDATA — Number
of data points in the z-direction. (Input)
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).
MDF — Middle
dimension of FDATA exactly as
specified in the dimension statement of the calling program.
(Input)
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).
IWK — Work array of length MAX(NXDATA, NYDATA, NZDATA).
2. Informational errors
Type Code
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
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |