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.

@H. P. Mark , actually tried to wrap all of that in a module, but all those C comments broke my compilation. I'm trying to write a test case while cradling an outraged baby in other arm -- i guess that puts me into a severely defocused fortran programmer class ;) -- and couldn't be bothered to look through gfortran info while was writing the original post. — Mali Remorker, Mar 19 '17 at 12:42
You can put everything in a module and still use the fixed form source. Of course you cannot cut and paste it into a free form source .f90 file. But there is no problem just adding `module` and `end module` into the original source file. — Vladimir F Героям слава, Mar 19 '17 at 14:09

roygvib · Accepted Answer · 2017-03-19T14:49:08.480

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

I think it would be very useful if d1mach.f in the original Slatec site includes a "stop" statement rather than silent return (e.g., a line like 'stop "d1mach needs manual update"'). This way, we notice the need to modify it automatically... — roygvib, Mar 19 '17 at 15:08

DGAMIC netlib function to calculate incomplete Gamma exits with an error

1 Answers1