BSLS3
Computes a three-dimensional tensor-product spline approximant using least squares, 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 nondecreasing.
YDATA — Array of length NYDATA containing the data points in the y-direction. (Input)
YDATA must be nondecreasing.
ZDATA — Array of length NZDATA containing the data points in the z-direction. (Input)
ZDATA must be nondecreasing.
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)
KYORD — Order of the spline in the y-direction. (Input)
KZORD — Order of the spline in the z-direction. (Input)
XKNOT — Array of length KXORD + NXCOEF containing the knots in the x-direction. (Input)
XKNOT must be nondecreasing.
YKNOT — Array of length KYORD + NYCOEF containing the knots in the y-direction. (Input)
YKNOT must be nondecreasing.
ZKNOT — Array of length KZORD + NZCOEF containing the knots in the z-direction. (Input)
ZKNOT must be nondecreasing.
BSCOEF — Array of length NXCOEF*NYCOEF*NZCOEF that contains the tensor product
B-spline coefficients. (Output)
Optional Arguments
NXDATA — Number of data points in the x-direction. (Input)
NXDATA must be greater than or equal to NXCOEF.
Default: NXDATA = size (XDATA,1).
NYDATA — Number of data points in the y-direction. (Input)
NYDATA must be greater than or equal to NYCOEF.
Default: NYDATA = size (YDATA,1).
NZDATA — Number of data points in the z-direction. (Input)
NZDATA must be greater than or equal to NZCOEF.
Default: NZDATA = size (ZDATA,1).
LDFDAT — Leading dimension of FDATA exactly as specified in the dimension statement of the calling program. (Input)
Default: LDFDAT = size (FDATA,1).
MDFDAT — Second dimension of FDATA exactly as specified in the dimension statement of the calling program. (Input)
Default: MDFDAT = size (FDATA,2).
NXCOEF — Number of B-spline coefficients in the x-direction. (Input)
Default: NXCOEF = size (XKNOT,1) - KXORD.
NYCOEF — Number of B-spline coefficients in the y-direction. (Input)
Default: NYCOEF = size (YKNOT,1) - KYORD.
NZCOEF — Number of B-spline coefficients in the z-direction. (Input)
Default: NZCOEF = size (ZKNOT,1) - KZORD.
XWEIGH — Array of length NXDATA containing the positive weights of XDATA. (Input)
Default: XWEIGH = 1.0.
YWEIGH — Array of length NYDATA containing the positive weights of YDATA. (Input)
Default: YWEIGH = 1.0.
ZWEIGH — Array of length NZDATA containing the positive weights of ZDATA. (Input)
Default: ZWEIGH = 1.0.
FORTRAN 90 Interface
Generic: CALL BSLS3 (XDATA, YDATA, ZDATA, FDATA, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, BSCOEF [, …])
Specific: The specific interface names are S_BSLS3 and D_BSLS3.
FORTRAN 77 Interface
Single: CALL BSLS3 (NXDATA, XDATA, NYDATA, YDATA, NZDATA, ZDATA, FDATA, LDFDAT, MDFDAT, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, NXCOEF, NYCOEF, NZCOEF, XWEIGH, YWEIGH, ZWEIGH, BSCOEF)
Double: The double precision name is DBSLS3.
Description
The routine BSLS3 computes the coefficients of a tensor-product spline least-squares approximation to weighted tensor-product data. The input for this subroutine consists of data vectors to specify the tensor-product grid for the data, three vectors with the weights, the values of the surface on the grid, and the specification for the tensor-product spline. The grid is specified by the three vectors x = XDATA, y = YDATA, and z = ZDATA of length k = NXDATA, l = NYDATA , and m = NZDATA, respectively. A three-dimensional array f = FDATA contains the data values which are to be fit. The three vectors wx = XWEIGH, wy = YWEIGH, and wz = ZWEIGH contain the weights for the weighted least-squares problem. The information for the approximating tensor-product spline must also be provided. This information is contained in kx = KXORD, tx = XKNOT, and K = NXCOEF for the spline in the first variable, in ky = KYORD, ty = YKNOT and L = NYCOEF for the spline in the second variable, and in kz = KZORD, tz = ZKNOT and M = NZCOEF for the spline in the third variable.
The coefficients of the resulting tensor product spline are returned in c = BSCOEF, which is an K × L × M array. The procedure computes coefficients by solving the normal equations in tensor-product form as discussed in de Boor (1978, Chapter 17). The interested reader might also want to study the paper by E. Grosse (1980).
The final result produces coefficients c minimizing
where the function Bstu is the tensor-product of three B-splines of order kx, ky, and kz. Specifically, we have
The spline
can be evaluated at one point using BS3VL and its partial derivatives can be evaluated using BS3DR. If the values on a grid are desired then we recommend BS3GD.
Comments
1. Workspace may be explicitly provided, if desired, by use of B2LS3/DB2LS3. The reference is:
CALL B2LS3 (NXDATA, XDATA, NYDATA, NZDATA, ZDATA, YDATA, FDATA, LDFDAT, KXORD, KYORD, KZORD, XKNOT, YKNOT, ZKNOT, NXCOEF, NYCOEF, NZCOEF, XWEIGH, YWEIGH, ZWEIGH, BSCOEF, WK)
The additional argument is:
WK — Work array of length NYCOEF * (NZDATA + KYORD + NZCOEF) + NZDATA * (1 + NYDATA) + NXCOEF * (KXORD + NYDATA * NZDATA) + KZORD * NZCOEF + 3 * MAX0(KXORD, KYORD, KZORD).
2. Informational errors
|
Type |
Code |
Description |
|
3 |
13 |
There may be less than one digit of accuracy in the least squares fit. Try using higher precision if possible. |
|
4 |
7 |
Multiplicity of knots cannot exceed the order of the spline. |
|
4 |
8 |
The knots must be nondecreasing. |
|
4 |
9 |
All weights must be greater than zero. |
|
4 |
10 |
The data point abscissae must be nondecreasing. |
|
4 |
11 |
The smallest element of the data point array must be greater than or equal to the K_ORDth knot. |
|
4 |
12 |
The largest element of the data point array must be less than or equal to the (N_COEF + 1)st knot. |
Example
The data for this example arise from the function e(y - z) sin(x + y) + ɛ on the rectangle [0, 3] × [0, 2] × [0, 1]. Here, ɛ is a uniform random variable with range [−.5, .5]. We sample this function on a 4 × 3 × 2 grid and then try to recover it by using tensor-product cubic splines in all variables. We print out the values of the function e(y - z) sin(x + y) on a 4 × 3 × 2 grid and compare these values with the values of the tensor-product spline that was computed using the IMSL routine BSLS3.
USE BSLS3_INT
USE RNSET_INT
USE RNUNF_INT
USE UMACH_INT
USE BS3GD_INT
IMPLICIT NONE
INTEGER KXORD, KYORD, KZORD, LDFDAT, MDFDAT, NXCOEF, NXDATA,&
NXVAL, NYCOEF, NYDATA, NYVAL, NZCOEF, NZDATA, NZVAL
PARAMETER (KXORD=4, KYORD=4, KZORD=4, NXCOEF=8, NXDATA=15,&
NXVAL=4, NYCOEF=8, NYDATA=15, NYVAL=3, NZCOEF=8,&
NZDATA=15, NZVAL=2, LDFDAT=NXDATA, MDFDAT=NYDATA)
!
INTEGER I, J, K, NOUT
REAL BSCOEF(NXCOEF,NYCOEF,NZCOEF), EXP, F,&
FDATA(NXDATA,NYDATA,NZDATA), FLOAT, RNOISE,&
SIN, SPXYZ(NXVAL,NYVAL,NZVAL), X, XDATA(NXDATA),&
XKNOT(NXCOEF+KXORD), XVAL(NXVAL), XWEIGH(NXDATA), Y,&
YDATA(NYDATA), YKNOT(NYCOEF+KYORD), YVAL(NYVAL),&
YWEIGH(NYDATA), Z, ZDATA(NZDATA),&
ZKNOT(NZCOEF+KZORD), ZVAL(NZVAL), ZWEIGH(NZDATA)
INTRINSIC EXP, FLOAT, SIN
! Define a function
F(X,Y,Z) = EXP(Y-Z)*SIN(X+Y)
!
CALL RNSET (1234579)
CALL UMACH (2, NOUT)
! Set up knot sequences
! X-knots
DO 10 I=1, NXCOEF - KXORD + 2
XKNOT(I+KXORD-1) = 3.0*(FLOAT(I-1)/FLOAT(NXCOEF-KXORD+1))
10 CONTINUE
DO 20 I=1, KXORD - 1
XKNOT(I) = XKNOT(KXORD)
XKNOT(I+NXCOEF+1) = XKNOT(NXCOEF+1)
20 CONTINUE
! Y-knots
DO 30 I=1, NYCOEF - KYORD + 2
YKNOT(I+KYORD-1) = 2.0*(FLOAT(I-1)/FLOAT(NYCOEF-KYORD+1))
30 CONTINUE
DO 40 I=1, KYORD - 1
YKNOT(I) = YKNOT(KYORD)
YKNOT(I+NYCOEF+1) = YKNOT(NYCOEF+1)
40 CONTINUE
! Z-knots
DO 50 I=1, NZCOEF - KZORD + 2
ZKNOT(I+KZORD-1) = 1.0*(FLOAT(I-1)/FLOAT(NZCOEF-KZORD+1))
50 CONTINUE
DO 60 I=1, KZORD - 1
ZKNOT(I) = ZKNOT(KZORD)
ZKNOT(I+NZCOEF+1) = ZKNOT(NZCOEF+1)
60 CONTINUE
! Set up X-grid.
DO 70 I=1, NXDATA
XDATA(I) = 3.0*(FLOAT(I-1)/FLOAT(NXDATA-1))
70 CONTINUE
! Set up Y-grid.
DO 80 I=1, NYDATA
YDATA(I) = 2.0*(FLOAT(I-1)/FLOAT(NYDATA-1))
80 CONTINUE
! Set up Z-grid
DO 90 I=1, NZDATA
ZDATA(I) = 1.0*(FLOAT(I-1)/FLOAT(NZDATA-1))
90 CONTINUE
! Evaluate the function on the grid
! and add noise.
DO 100 I=1, NXDATA
DO 100 J=1, NYDATA
DO 100 K=1, NZDATA
RNOISE = RNUNF()
RNOISE = RNOISE - 0.5
FDATA(I,J,K) = F(XDATA(I),YDATA(J),ZDATA(K)) + RNOISE
100 CONTINUE
! Use default weights equal to 1.0
!
! Compute least-squares
CALL BSLS3 (XDATA, YDATA, ZDATA, FDATA, KXORD, KYORD, KZORD, XKNOT, &
YKNOT, ZKNOT, BSCOEF)
! Set up grid for evaluation.
DO 110 I=1, NXVAL
XVAL(I) = FLOAT(I-1)
110 CONTINUE
DO 120 I=1, NYVAL
YVAL(I) = FLOAT(I-1)
120 CONTINUE
DO 130 I=1, NZVAL
ZVAL(I) = FLOAT(I-1)
130 CONTINUE
! Evaluate on the grid.
CALL BS3GD (0, 0, 0, XVAL, YVAL, ZVAL, KXORD, KYORD, KZORD, XKNOT, &
YKNOT, ZKNOT, BSCOEF, SPXYZ)
! Print results.
WRITE (NOUT,99998)
DO 140 I=1, NXVAL
DO 140 J=1, NYVAL
DO 140 K=1, NZVAL
WRITE (NOUT,99999) XVAL(I), YVAL(J), ZVAL(K),&
F(XVAL(I),YVAL(J),ZVAL(K)),&
SPXYZ(I,J,K), F(XVAL(I),YVAL(J),ZVAL(K)&
) - SPXYZ(I,J,K)
140 CONTINUE
99998 FORMAT (8X, 'X', 9X, 'Y', 9X, 'Z', 6X, 'F(X,Y,Z)', 3X,&
'S(X,Y,Z)', 4X, 'Error')
99999 FORMAT (' ', 3F10.3, 3F11.4)
END
X Y Z F(X,Y,Z) S(X,Y,Z) Error
0.000 0.000 0.000 0.0000 0.1987 -0.1987
0.000 0.000 1.000 0.0000 0.1447 -0.1447
0.000 1.000 0.000 2.2874 2.2854 0.0019
0.000 1.000 1.000 0.8415 1.0557 -0.2142
0.000 2.000 0.000 6.7188 6.4704 0.2484
0.000 2.000 1.000 2.4717 2.2054 0.2664
1.000 0.000 0.000 0.8415 0.8779 -0.0365
1.000 0.000 1.000 0.3096 0.2571 0.0524
1.000 1.000 0.000 2.4717 2.4015 0.0703
1.000 1.000 1.000 0.9093 0.8995 0.0098
1.000 2.000 0.000 1.0427 1.1330 -0.0902
1.000 2.000 1.000 0.3836 0.4951 -0.1115
2.000 0.000 0.000 0.9093 0.8269 0.0824
2.000 0.000 1.000 0.3345 0.3258 0.0087
2.000 1.000 0.000 0.3836 0.3564 0.0272
2.000 1.000 1.000 0.1411 0.1905 -0.0494
2.000 2.000 0.000 -5.5921 -5.5362 -0.0559
2.000 2.000 1.000 -2.0572 -1.9659 -0.0913
3.000 0.000 0.000 0.1411 0.4841 -0.3430
3.000 0.000 1.000 0.0519 -0.4257 0.4776
3.000 1.000 0.000 -2.0572 -1.9710 -0.0862
3.000 1.000 1.000 -0.7568 -0.8479 0.0911
3.000 2.000 0.000 -7.0855 -7.0957 0.0101
3.000 2.000 1.000 -2.6066 -2.1650 -0.4416
