Функция DLAMIC netlib для расчета неполных выходов гаммы с ошибкой

Question

Мне нужна процедура для расчета неполной гамма-функции. Конечно, я попробовал маршрут netlib и нашел функцию dgamic. Однако после составления следующей тестовой программы

program test_dgamic
  implicit none
  interface 
     double precision function dgamic(in1,in2)
       double precision, intent(in) :: in1,in2
     end function dgamic
  end interface 


  print *, 'dgamic:', dgamic(1.d0,1.d0)
end program test_dgamic

с gfortran версии 6.2.0 вот так

gfortran main.f90 -o main dgamic.f d9lgic.f d9lgit.f d9gmic.f d9gmit.f dlgams.f dlngam.f dgamma.f d9lgmc.f dcsevl.f dgamlm.f initds.f d1mach.f xerclr.f xermsg.f xerprn.f xersve.f xgetua.f i1mach.f j4save.f xerhlt.f xercnt.f fdump.f

и работает, я получаю следующее сообщение об ошибке slatec

***MESSAGE FROM ROUTINE INITDS IN LIBRARY SLATEC.
 ***POTENTIALLY RECOVERABLE ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  Chebyshev series too short for specified accuracy
 *  ERROR NUMBER = 1
 *   
 ***END OF MESSAGE

 ***JOB ABORT DUE TO UNRECOVERED ERROR.
0          ERROR MESSAGE SUMMARY
 LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
 SLATEC     INITDS     Chebyshev series too         1         1         1

Note: The following floating-point exceptions are signalling: IEEE_DIVIDE_BY_ZERO

Кто-нибудь понял, как этого избежать? Судя по всему, ошибка выглядит как недостаток дизайна.

Answer 1

Кажется, что проблема (опять же) связана с d1mach.f в Slatec, потому что нам нужно раскомментировать соответствующий раздел этого файла вручную. На практике удобнее использовать модифицированную версию d1mach.f, доступную на сайте BLAS (см. Эту страницу). Так что процедура может быть что-то вроде:

1) скачать slatec_src.tar.gz с оригинального сайта

2) скачать модифицированные (BLAS) версии d1mach.f и т. Д. И использовать их вместо версий Slatec

rm -f i1mach.f r1mach.f d1mach.f
wget http://www.netlib.org/blas/i1mach.f 
wget http://www.netlib.org/blas/r1mach.f
wget http://www.netlib.org/blas/d1mach.f 

3) и совместить, например, с тестовой программой

program main
    implicit none
    external dgamic
    double precision dgamic, a, x, y

    a = 1.0d0
    x = 1.0d0
    y = dgamic( a, x )

    print *, "a = ", a
    print *, "x = ", x
    print *, "y(slatec)        = ", y
    print *, "y(exact for a=1) = ", exp( -x )
end program 

который дает

 a =    1.0000000000000000     
 x =    1.0000000000000000     
 y(slatec)        =   0.36787944117144233     
 y(exact for a=1) =   0.36787944117144233

Для сравнения, если мы используем версию d1mach.f для Slatec, мы получаем ту же ошибку

 ***MESSAGE FROM ROUTINE INITDS IN LIBRARY SLATEC.
 ***POTENTIALLY RECOVERABLE ERROR, PROG ABORTED, TRACEBACK REQUESTED
 *  Chebyshev series too short for specified accuracy
 *  ERROR NUMBER = 1
 ...

потому что точность не установлена ​​в d1mach.f (нам нужно раскомментировать необходимый раздел, например, помеченный как "IBM PC", плюс некоторые модификации для старого Фортрана, что может быть утомительно...)