DGAMIC netlib function to calculate incomplete Gamma exits with an error

Question

I need a procedure to calculate the incomplete Gamma function. Of course, I've tried the netlib route and found the dgamic function. However, after compiling the following test program

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

with gfortran version 6.2.0 like this

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

and running, I get the following slatec error message

***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

Has anyone got a clue how to avoid this? From the looks of the error, it looks like a design flaw.

Answer 1

It seems that the problem is (again) due to d1mach.f in Slatec, because we need to uncomment an appropriate section of that file manually. In practice, it is more convenient to use a modified version of d1mach.f available from the BLAS site (see this page). So the procedure may be something like:

1) download slatec_src.tar.gz from the original site

2) download modified (BLAS) versions of d1mach.f etc and use them instead of the Slatec versions

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) and comiple, e.g., with a test program

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 

which gives

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

For comparison, if we use the Slatec version of d1mach.f, we get the same error

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

because the precision is not set in d1mach.f (we need to uncomment a necessary section, e.g. labeled as "IBM PC", plus some modifications for legacy Fortran, which could be tedious...)