Problems with GAMMA functions?

[Angelo Graziosi <Angelo dot Graziosi at roma1 dot infn dot it>]

I have found this on Cygwin but, perhaps it is general.
=======================================================

Usually I have the habit to build, on Cygwin, CERNLIB with GFortran (but
also G77 and G95).

CERNLIB has its own tests, between which, 89 mathematical tests.

Before that the set of GAMMA functions were added in GFortran, the build
of CERNLIB completed succeffully all tests. The last I tried was with
GFortran 4.3.0-20070823-127745.

But at least with gfc 4.3.0-20070903-128061, the CERNLIB build gives two
failure:

   Test# 14 ( C302 ): *** failed ***
   Test# 66 ( E408 ): *** failed ***


These tests use GAMMAs!


I have syntesized the C302 test (gfc_gamma_test_c302m.F) below.

If one build it using the (default) internal gfc GAMMAs, it fails.

If, instead, one links with the EXTERNAL GAMMAs (from a Cernlib build made
with GFortran 4.3.0 .LE. 20070823), it is completed successfully.

In this case, one have to uncomment the line

     C      external gamma,dgamma


(The Cygwin build of CERNLIB with GFC-4.3.0-20070727 can be found here
http://www.webalice.it/angelo.graziosi/cygwin/cernlib/Cernlib.html).

Interesting is also to compare, in the above cases, the 'fort.10' file
generated by 'gfc_gamma_test_c302m'.



Cheers,
   Angelo.



gfc_gamma_test_c302m.F
===========================================================================
C
C gfortran -O0 -funroll-loops -fomit-frame-pointer -fno-second-underscore
C          -fno-automatic gfc_gamma_test_c302m.F -o gfc_gamma_test_c302m
C
C
C Uncommenting the 'external' you must build in this way:
C
C gfortran -O0 -funroll-loops -fomit-frame-pointer -fno-second-underscore
C          -fno-automatic gfc_gamma_test_c302m.F `cernlib mathlib`
C          -o gfc_gamma_test_c302m
C
      program main
      implicit none

      integer LIN,LOUT
      COMMON/IOLUNS/LIN,LOUT

      integer NTEST,NFAIL,IRC
      COMMON/GTSTAT/NTEST,NFAIL,IRC

      LOGICAL LTEST1
      COMMON /C302LT1/LTEST1

      LOUT = 10
      NTEST = 0
      NFAIL = 0

      call c302m

      end
C
C ===============================================================
C
      SUBROUTINE C302M
C     Program to test the MATHLIB routines GAMMA, DGAMMA and
C     QGAMMA (C302)

      LOGICAL LTEST, LTEST1,LTEST2
      COMMON /C302LT1/LTEST1

      COMMON/IOLUNS/LIN,LOUT
      COMMON/GTSTAT/NTEST,NFAIL,IRC

      CALL HEADER('C302',0)
      LTEST=.TRUE.
      LTEST1=.TRUE.
      LTEST2=.TRUE.
      CALL C302D
        LTEST=LTEST .AND. LTEST1
      IRC=ITEST('C302',LTEST)
      CALL PAGEND('C302')
      RETURN
      END
C ---------------------------------------------------------------
      SUBROUTINE C302D

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      REAL GAMMA
C
C      external gamma,dgamma
C
      CHARACTER*6 TFUNC(2)

      COMMON/IOLUNS/LIN,LOUT
      COMMON/GTSTAT/NTEST,NFAIL,IRC
C
      PARAMETER (HALF = 5D-1, PI=3.14159 26535 89793 24D0)
      DIMENSION C(0:20)
      LOGICAL LTEST1
      DIMENSION TOL(2)
      COMMON /C302LT1/LTEST1
C
C     Set maximum error allowed for test to be considered successful
C

      DATA TOL/1D-6, 5D-14/
      DATA TFUNC/'GAMMA','DGAMMA'/
      LTEST1=.TRUE.

      C(0)=1
      DO 2 N = 1,20
    2 C(N)=(2*N-1)*C(N-1)/2

      NF=2

      DO 1000 JF=1,NF
         ERRMAX= 0.0D0
        RERRMAX= 0.0E0
         WRITE(LOUT,'(/10X,''Test of C302 '',A)') TFUNC(JF)
         WRITE(LOUT,'(/9X,''X  '',7X,''Exact'',25X,''Calculated'',
     + 14X,''Rel. Error'')')

      DO 1 N =  1,20
                   X=N+HALF
                   T=C(ABS(N))*SQRT(PI)

      IF(JF.EQ.1) THEN
                        DR=GAMMA( SNGL(X) )
         IF(DR .NE. 0) ER=ABS(SNGL(DR-T)/SNGL(DR) )
         WRITE(LOUT,'(1X,F10.1,2E27.9,5X,E10.1)')
     +        SNGL(X),SNGL(T),SNGL(DR),SNGL(ER)

      ENDIF
      IF(JF.EQ.2) THEN
                       DR=DGAMMA(X)
         IF(DR .NE. 0) ER=ABS((DR-T)/DR)
         WRITE(LOUT,'(1X,F10.1,2E27.18,5X,E10.1)') X,T,DR,ER
      ENDIF

      ERRMAX= MAX( ERRMAX,ER )
    1 CONTINUE

      ETOL=TOL(JF)

      WRITE(LOUT,'(/''Largest relative error was'',E10.1)') ERRMAX
      LTEST1=LTEST1.AND.(ERRMAX.LE.ETOL)

      WRITE(LOUT,'(/''Testing error messages:'')')

      IF(JF.EQ.1) DR=GAMMA(-SNGL(HALF))
      IF(JF.EQ.2) DR=DGAMMA(-HALF)

1000  CONTINUE
      RETURN
      END
C
C ================== Aux routines ===============================
C
      SUBROUTINE HEADER(CODE,MODE)
C     This routine prints a page header with a testing routine name
C     message.
      CHARACTER*(*) CODE

      COMMON/IOLUNS/LIN,LOUT
      COMMON/GTSTAT/NTEST,NFAIL,IRC

      NTEST=NTEST+1
      IF(MODE.EQ.1) PRINT      1000, NTEST,CODE
                    WRITE(LOUT,1001) CODE
      RETURN
1000  FORMAT(' Test#',I3,' ( ',A,' ):     started')
1001  FORMAT('1',25X,30('*')/26X,'**   Testing Routine ',A,3X,'**'/
     +       26X,30('*'))
      END
C ---------------------------------------------------------------
      SUBROUTINE PAGEND(CODE)
C     This subroutine prints a page end message
      CHARACTER*(*) CODE

      COMMON/IOLUNS/LIN,LOUT
      COMMON/GTSTAT/NTEST,NFAIL,IRC

                    WRITE(LOUT,1001) CODE
      RETURN
1001  FORMAT(/26X,30('*')/26X,'**   End of Test of  ',A,3X,'**'/
     +       26X,30('*'))
      END
      FUNCTION ITEST(CODE,TEST)

      COMMON/IOLUNS/LIN,LOUT
      COMMON/GTSTAT/NTEST,NFAIL,IRC

      CHARACTER*(*) CODE
      LOGICAL TEST
      IF(TEST) THEN
        PRINT 1000,NTEST,CODE
        IF (LOUT .NE. 6) WRITE(LOUT,1000) NTEST,CODE
        ITEST=0
      ELSE
        PRINT 1001,NTEST,CODE
        IF (LOUT .NE. 6) WRITE(LOUT,1001) NTEST,CODE
        ITEST=1
      ENDIF
      NFAIL=NFAIL+ITEST
1000  FORMAT(' Test#',I3,' ( ',A,' ):     completed successfully')
1001  FORMAT(' Test#',I3,' ( ',A,' ): *** failed ***')
      RETURN
      END
===========================================================================