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.
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...)