用于计算不完整 Gamma 的 DGAMIC netlib 函数退出时出现错误

Question

我需要一个程序来计算不完整的 Gamma 函数。当然，我试过 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

似乎问题（再次）是由于 Slatec 中的 d1mach.f，因为我们需要手动取消注释该文件的适当部分。实际上，使用 BLAS 站点提供的 d1mach.f 的修改版本更方便（参见 this 页面）。所以程序可能是这样的：

1) 从 original site

下载 slatec_src.tar.gz

2) 下载 d1mach.f 等的修改 (BLAS) 版本并使用它们代替 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"，再加上对遗留 Fortran 的一些修改，这可能很乏味...）

用于计算不完整 Gamma 的 DGAMIC netlib 函数退出时出现错误

DGAMIC netlib function to calculate incomplete Gamma exits with an error

fortran

gfortran

netlib