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.

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.