PDA

View Full Version : Unexpected behavior from BSINT



mecej4
01-26-2014, 01:54 PM
Consider this example, which is almost the same as the manual example for BSINT, but with slightly different input data points, a different function to approximate and a different set of points at which the spline is to be evaluated.

USE bsint_int
USE bsnak_int
USE bsval_int

IMPLICIT NONE
INTEGER korder, ndata, nknot
PARAMETER (korder=3, ndata=21, nknot=ndata+korder)
!
INTEGER i, ncoef
REAL bscoef(ndata), bt, f, fdata(ndata), sqrt, x, xdata(ndata), &
xknot(nknot), xt
INTRINSIC sqrt
! Define function
f(x) = sqrt(1-x)
! Set up interpolation points
DO i = 1, ndata
xt = (i-11.0)/10.0
xdata(i) = xt
fdata(i) = f(xdata(i))
END DO
WRITE (*, '(1x,A,2x,1pE14.7,/)') 'XDATA(1) :', xdata(1)
! Generate knot sequence
CALL bsnak(ndata, xdata, korder, xknot)
! Interpolate
CALL bsint(ndata, xdata, fdata, korder, xknot, bscoef)
!
WRITE (*, 110)
! Print on a finer grid
ncoef = ndata
DO i = 1, 4
xt = 2.0*((i-1.0)/3.0) - 1.0
! Evaluate spline
bt = bsval(xt, korder, xknot, ncoef, bscoef)
WRITE (*, 100) xt, bt, f(xt) - bt
END DO
100 FORMAT (' ', F11.7, 10X, F8.4, 12X, F11.6)
110 FORMAT (/, 6X, 'X', 19X, 'S(X)', 18X, 'Error', /)
END PROGRAM
When I compiled and ran this with the Intel Fortran compiler, I obtained the interpolated value as 0.0 when default compiler optimization was allowed, and 1.4141 when optimization was disabled. How to explain this?

Two factors contribute to this failure:

1) The optimized code for the statement
xt = (i-11.0)/10.0gives xt = -9.9999994E-01 when i = 1, whereas the unoptimized code gives xt = -1.0000000E+00. Considering that we are using 4-byte REALs, this is reasonable. Unfortunately, when we try to obtain the spline value with xt = -1.0 later, with the optimized EXE BSVAL decides that this value is outside the interval covered by the spline.

2) It is an un/under-documented property of BSVAL that it returns 0.0 when asked to evaluate the B-spline at a point outside the interval covered. No error message is given, nor could I find this property described in the documentation. To catch this error, one would have to program more defensively, and explicitly call the error service routines ERSET, IERCD, etc., after calling BSVAL to find if an error occurred.

Something similar to this happens with the IMSL manual example for BS3IN. The example code given by Pate in the thread 3-dimension interpolation (http://forums.roguewave.com/showthread.php?553-3-dimension-interpolation) fails in the same fashion when compiled with /Ot or higher optimization level with the Intel compiler.

The user could, of course, add checks to the last loop in the program above

if(xt < xknot(1))then
xt = xknot(1)
elseif(xt > xknot(ndata+korder))then
xt = xknot(ndata+korder)
endif
but the user could be made aware of the need for such measures in the documentation.