Помогите пожалуйста.
MS Visual Studio 2005 + Intel Fortran 9.1 + IMSL FNL 5.0Никак не могу разобраться с IMSL-евской функцией PPVAL. Ведёт себя непредсказуемо.
Пример:
Хочу аппроксимировать функцию f(x)=dsin(Pi/2.d0*x) на промежутке [0,2] с помощью B-spline 4-го порядка. Для этого использую:
BSNAK - Computes the “not-a-knot” spline knot sequence.
BSINT - Computes the spline interpolant, returning the B-spline coefficients.
BSVAL - Evaluates a spline, given its B-spline representation.
В принципе, этого достаточно для построения сплайна. Но так как мне надо будет вычислять сплайн во многих промежуточных точках и важна скорость вычисления, то надо преобразовать сплайн в полином.
BSCPP - Converts a spline in B-spline representation to piecewise polynomial representation.
PPVAL - Evaluates a piecewise polynomial.
Код:
program spline
use BSINT_INT
use BSNAK_INT
use BSCPP_INT
use PPVAL_INT
use BSVAL_INT
implicit none
real(8) f
real(8) :: xi1=0.0d0, xi2=2.0d0
integer, parameter :: korder=4, ndata=20, ncoef=20, nknot=korder+ncoef
real(8) :: break(ncoef), ppcoef(korder,ncoef), bscoef(ncoef)
real(8), dimension(nknot) :: xknot
integer i,nppcf
real(8) x, xdata(ndata), fdata(ndata)
real(8) z1,z2,z3
real(8), parameter :: Pi=3.141592653589793238d0
f(x)=dsin(Pi/2.d0*x)
do i=1, ndata
x=xi1+(xi2-xi1)*dfloat(i-1)/dfloat(ndata-1)
xdata(i) = x
fdata(i) = f(x)
enddo
CALL D_BSNAK (NDATA, XDATA, KORDER, XKNOT)
CALL D_BSINT (NCOEF, XDATA, FDATA, KORDER, XKNOT, BSCOEF)
CALL D_BSCPP (KORDER, XKNOT, NCOEF, BSCOEF, NPPCF, BREAK, PPCOEF)
! Первый блок
I=2*ndata-3
x = dFLOAT(I)/dFLOAT(NDATA-1)
z1=f(x)
z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
z3=D_PPVAL(x,BREAK,PPCOEF)
print*,'x=',x
print*,z1,' ',z2,' ',z3
print*
! Второй блок
x = dFLOAT(2*ndata-3)/dFLOAT(NDATA-1)
z1=f(x)
z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
z3=D_PPVAL(x,BREAK,PPCOEF)
print*,'x=',x
print*,z1,' ',z2,' ',z3
print*
! Третий блок
I=2*ndata-3
x = dFLOAT(I)/dFLOAT(NDATA-1)
z1=f(x)
z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
z3=D_PPVAL(x,BREAK,PPCOEF)
print*,'x=',x
print*,z1,' ',z2,' ',z3
print*
pause 'FINISH: press ENTER.'
end program spline
Результат:
Код:
x=1.94736842105263
8.257934547233223E-002 8.258268483722089E-002 [b]0.000000000000000E+000[/b]
x=1.94736842105263
8.257934547233223E-002 8.258268483722089E-002 8.258268483722092E-002
x=1.94736842105263
8.257934547233223E-002 8.258268483722089E-002 8.258268483722092E-002
Вот те жирные нули в первой строке результата выглядят совершенно странно. Должно быть то же самое, что строкой ниже, то-есть второй столбик (рассчитанным по BSVAL) должен практически совпадать с третьим (рассчитанный по PPVAL). У мене IMSL 5.0. Интересно, что покажет IMSL 6.0.
Что интересно, при некоторых x функция PPVAL в каждом блоке считает правильно (например x=1.1234d0)! А при некоторых нет.
Обратите внимание, первый и третий блоки идентичны, а результат расчёта PPVAL в них разный!
При компиляции оптимизация отсутствовала. Было бы хорошо, если бы кто-то запустил этот код у себя и закинул сюда результат. Было бы хорошо знать, исправлен ли этот глюк в последующих версиях фортрана и ИМСЛя. И вообще следует разобраться, чей глюч: фортрана, ИМСЛя или мой.
PS.
Я уже поднял этот вопрос на ру-борде, но ребята там сильно отклонилиь от темы.
