глюк IMSL???

[godsdog]

Помогите пожалуйста.
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.
Я уже поднял этот вопрос на ру-борде, но ребята там сильно отклонилиь от темы.

[godsdog]

Проблема решена. Решение можно найти здесь и здесь.

[godsdog]

Коротко.
Следует явно указать опциональные параметры в PPVAL:

D_PPVAL(x,BREAK,PPCOEF,korder=korder,nintv=nppcf)