Different results for 32-Bit and 64-Bit versions

Using Silverfrost ftn95 for reactivating older Fortran programs which run in the 80s on an IBM 3032 Mainframe generally works fine, when replacing IBM specific commands/conventions. Currently, however, I'm stuck with a program calculating matrix elements. What puzzles me is the fact, that the 32Bit-Version delivers different results compared with the results of the 64Bit-Version (I'm using Plato for compilation and SDBG version 8.70.00). Whereas the 32Bit-version produces correct results, the ones calculated with the 64Bit version are obvioulsy wrong. Trying to isolate the problem, kept me busy the last days, but is not so easy due to 20 subroutines and many function calls. Since I didn't succeed to isolate the root cause of the differences between the results of the 32- and 64-Bit versions up to now, I would like to ask the following questions: Did anybody experience similar deviations between results of 32- and 64-Bit versions? Can certain compiler options help to isolate the problem? Please note, that I have no experience with command line use of SDBG, but use the Plato environment. @mecej4: Thank you for the response and information contained therein concerning possible reasons for the different results. As suggested I add the part of source code here, which yields the different results. Update: Meanwhile I could isolate the problem to origin from the subroutine wavefun. When calling the Whittaker function, variable GAMBA differs (as result from calling GAMCAR (Gamma function). The input used for WAVEFUN is as follows: -1.500D+00 -1 1.840D+02 3.286D-01 3.002D-01 5.000D+03 0 0.000D+0 3.2858593764968097D-01 3.0015179886716270D-01

The subroutine WAVEFUN (with many writes to file 14 for test purposes) is listed below:

C     ******************************************************************
C     ******************************************************************
C      SUBROUTINE WAVFUN(E,KAPPA,Z,U1,U2,R0,ITEST,DELTA)                 
C     ******************************************************************
C     ******************************************************************
C     THIS PROGRAM COMPUTES THE RADIAL KONTINUUM WAVEFUNCTIONS OF THE   
C     DIRAC EQUATION FOR EACH Z WITH RESPECT OF AN EXTENDED NUCLEUS.    
C     U1 AND U2 ARE THE RADIAL KONTINUUMWAVEFUNCTIONS.                  
C     SEE FOR EXAMPLE EQ. 5.4' OF REF 2.                                
C     THE RADIUS OF THE NUCLEUS IS GIVEN BY R=R0*A**(1./3.)             
C     UNITS ARE H/(2.*PI) = C = M = 1.                                  
C     RF IS R IN FM.                                                    
C     REF 1.  B. MUELLER, J. RAFELSKI AND W. GREINER                    
C     INSTITUT FUER THEORETISCHE PHYSIK - FRANKFURT                     
C     PREPRINT (1973)                                                   
C     REF 2.  E. M. ROSE, RELATIVISTISCHE ELEKTRONENTHEORIE.            
C     REF 3.  B. MUELLER, J. RAFELSKI, AND W. GREINER                   
C             Z. PHYSIK 257,183 (1972)                                  
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 PHI1,PHI2,FAK3,FAK4,GAM,CDSQRT,ABC,FAK1,FAK2,DCONJG
      COMMON /ZELLE/A1,A2
      OPEN (5,file='D:\FORTRAN\MEDDAT\TST\FT05WF.DAT',status='READONLY')      
      OPEN (6,file='D:\FORTRAN\MEDDAT\TST\FT06WF.DAT',status='OLD')
      OPEN (14,file='D:\FORTRAN\MEDDAT\TST\FT14WF.DAT',status='OLD')
      READ (5,*) E,KAPPA,Z,U1,U2,R0,ITEST,DELTA,A1,A2
      WRITE (14,1680) E,KAPPA,Z,U1,U2,R0,ITEST,DELTA
1680  FORMAT (1X,'Entering WAVEFUN, E,KAPPA,Z,U1,U2,R0,ITEST,DELTA',
     &1X,1PD10.3,I3/1X,1P4D15.3,I6,1PD10.3)          
      COMWAV=386.1592D0                                                 
      ALPHA=1./137.03602D0                                              
      PI=3.1415926535897D0                                              
      ETA=0.D0                                                          
      ZA=Z*ALPHA                                                        
      S=DFLOAT(IABS(KAPPA))                                             
      ABC=DCMPLX(S*S-ZA*ZA,0.D0)                                        
      GAM=CDSQRT(ABC)
C	2021 April 20, Write for test purpose only
      Write (14,1670) GAM,S,Z
1670  Format (1x,'Wavef GAM,S,Z',1P4D20.10)      
      P=DSQRT(E*E-1.)                                                   
      Y=ZA*E/P                                                          
      IF (E.LT.0) GO TO 20                                              
C     SEE EQ. 4.1 AND 4.2 OF REF1                                       
      FAK1=DCMPLX(DSQRT(E+1.),0.D0)                                     
      FAK2=FAK1                                                         
      FAK3=DCMPLX(0.D0,DSQRT(E-1.D0))                                   
      FAK4=-FAK3                                                        
      GO TO 30                                                          
   20 CONTINUE                                                          
      FAK1=DCMPLX(DSQRT(-E-1.),0.D0)                                    
      FAK2=FAK1                                                         
      FAK3=-DCMPLX(0.D0,DSQRT(-E+1.D0))                                 
      FAK4=-FAK3                                                        
   30 CONTINUE                                                          
C     DETERMINATION OF THE NORMALIZATION INSIDE THE NUCLEUS.            
      R=R0/COMWAV                                                       
      RS=R
C     2021, April 20 - Write only for test purposes      
      WRITE (14,1660) Kappa,ETA,Delta
1660  FORMAT (1x,'In Wavefun before Phi:     ',1X,I3,1P2D20.10)                                                                    
C     ------------------------------------------------------------------
      CALL PHI(R,Y,P,Z,R0,E,GAM,KAPPA,K,M,PHI1,ETA,DELTA)               
C     ------------------------------------------------------------------
C     2021, April 20 - Write only for test purposes      
      WRITE (14,1661) Kappa,ETA,Delta
1661  FORMAT (1x,'In Wavefun after Phi:      ',1X,I3,1P2D20.10)   
      PHI2=DCONJG(PHI1)                                                 
      U1=DREAL(FAK1*PHI1+FAK2*PHI2)                                     
      U2=DREAL(FAK3*PHI1+FAK4*PHI2)                                     
C     THE WAVEFUNCTIONS ARE READY.                                      
C     POSSIBLE TEST OF THE RIGHT ASYMPTOTIC.                            
C     SEE EQ. 3.7A AND 3.7B OF REF 3.                                   
C     T HAS TO BE 1 FOR LARGE R.                                        
      IF (ITEST.EQ.0) STOP 111                                           
      T=0.                                                              
      WRITE (6,501) Z,E,KAPPA                                           
  501 FORMAT (' Z =',F7.1,'   E =',E16.8,'    KAPPA =',I5/)             
      T1=U1                                                             
      T2=U2                                                             
      IF (E.LT.0.) GO TO 60                                             
      T1=T1/DSQRT((E+1.)/PI/P)                                          
      T2=T2/DSQRT((E-1.)/PI/P)                                          
      GO TO 70                                                          
   60 CONTINUE                                                          
      T1=T1/DSQRT((-E-1.)/PI/P)                                         
      T2=T2/DSQRT((-E+1.)/PI/P)                                         
   70 CONTINUE                                                          
      T1=T1*T1                                                          
      T2=T2*T2                                                          
      T=T1+T2                                                           
      RF=R*COMWAV                                                       
      WRITE(6,500) RF,U1,U2,T,DELTA                                     
500   FORMAT (' ASYMPTOTIK ', 4D16.7, 5X,D16.7)                         
      STOP                                                            
      END                                                               
C     ******************************************************************
C     ******************************************************************
      SUBROUTINE PHI (R,Y,P,Z,R0,E,GAM,KAPPA,K,M,PHI1,ETA,DELTA)        
C     ******************************************************************
C     ******************************************************************
C     THIS PROGRAM EVALUATES THE FUNCTION PHI1.                         
C     SEE EQU. 5.5 AND 6.5 OF REF 1.                                    
      IMPLICIT REAL*8 (A-H,O-Z)
C     2021, April 20, WHITT removed in next line                                         
      COMPLEX*16 EPLUS,EMIN,PHI1,ALP,X,A,B,C,D,XI,XM              
      COMPLEX*16 CDEXP,GAM,F1,F2,F3,CDSQRT,WI1,WI2,S1,S2                      
C	  2021-04-20 Common /Witz/ extended by S1,S2                                                     
      COMMON /WITZ/WI1,WI2,S1,S2
      PI=3.1415926535897D0
      EPLUS = 0.D0
      EMIN  = 0.D0
C     2021, April 20 - Write only for test purposes      
      WRITE (14,1662) Kappa,ETA,Delta
1662  FORMAT (1x,'In Wavefun, in Phi(begin): ',1X,I3,1P2D20.10)                                                     
      IF (DABS(DIMAG(GAM)).GT.1.D-6) GO TO 40                           
      I=1                                                               
C     ..................................................................
      EPLUS=CDSQRT(ALP(I,GAM,Y,E,KAPPA))                                
C     ..................................................................
      I=-1                                                              
C     ..................................................................
      EMIN=CDSQRT(ALP(I,GAM,Y,E,KAPPA))                                 
C     ..................................................................
   40 CONTINUE                                                          
      IF (DABS(ETA).GT.1.D-30) GO TO 50                                 
C     ------------------------------------------------------------------
      CALL TANETA(R0,E,GAM,KAPPA,Y,Z,ETA,EPLUS,EMIN,SINETA,COSETA)      
C     ------------------------------------------------------------------
   50 CONTINUE
C     2021, April 20 - Write only for test purposes      
      WRITE (14,1663) Kappa,ETA,Delta
1663  FORMAT (1x,'In Wavefun, Phi before Norm',1X,I3,1P2D20.10)                                                              
C     ------------------------------------------------------------------
      CALL NORM(P,GAM,Y,ETA,EPLUS,EMIN,ANORM,E,KAPPA,Z,SINETA,   
     &COSETA,DELTA)                                                           
C     ------------------------------------------------------------------
C     2021, April 20 - Write only for test purposes      
      WRITE (14,1664) Kappa,ETA,Delta
1664  FORMAT (1x,'In Wavfun, Phi after Norm: ',1X,I3,1P2D20.10)    
      XI=0.5+DCMPLX(0.D0,Y)                                             
      XI=-XI                                                            
      XM=GAM                                                            
      A=COSETA*EPLUS                                                    
      B=SINETA*EMIN                                                     
      IF (DABS(DIMAG(GAM)).LT.1.D-6) GO TO 30                           
      ALPHA=1.D0/137.03602D0                                            
      ZA=Z*ALPHA                                                        
      GS=DABS(DFLOAT(KAPPA))                                            
      GG=DSQRT(ZA*ZA-GS*GS)                                             
      F1=CDEXP(DCMPLX(0.D0,ETA))                                        
      F2=CDEXP(DCMPLX(0.D0,-ETA)-PI*GG)                                 
      F3=(GAM-DCMPLX(0.D0,Y))/(DFLOAT(KAPPA)+DCMPLX(0.D0,Y)/E)          
   30 CONTINUE                                                          
      X=DCMPLX(0.D0,2.*P*R)                                             
      C=WI1                                                             
      D=WI2                                                             
      IF (DABS(DIMAG(GAM)).GT.1.D-6) GO TO 20                           
      PHI1=(C*A+D*B)*ANORM/CDSQRT(X)                                    
      RETURN                                                            
   20 CONTINUE                                                          
      PHI1=(F1*C+F2*F3*D)*ANORM/CDSQRT(X)                               
      RETURN                                                            
      END                                                               
C     ******************************************************************
C     ******************************************************************
      SUBROUTINE NORM (P,GAM,Y,ETA,EPLUS,EMIN,ANORM,E,KAPPA,Z,SINETA,
     &COSETA,DELTA) 
C     ******************************************************************
C     ******************************************************************
C     THIS PROGRAM EVALUATES THE NORMALISATION FACTOR OF THE            
C     WAVEFUNCTIONS.                                                    
C     SEE EQ. 5.9 AND 6.7 OF REF 1.                                     
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 GAMCAR,EPLUS,EMIN,A,B,C,D,CDEXP,AA,BB,CC,GAM           
      COMMON /TEST/A                                                    
      PI=3.1415926535897D0                                              
      F=DEXP(0.5*PI*Y)                                                  
      G=2.*DSQRT(PI*P)                                                  
      IF (DABS(DIMAG(GAM)).GT.1.D-6) GO TO 10                           
      A=2.*GAM+1.                                                       
C     ..................................................................
      A=GAMCAR(A)                                                       
C     ..................................................................
      B=GAM+1.+DCMPLX(0.D0,Y)                                           
C     ..................................................................
      B=GAMCAR(B)                                                       
C     ..................................................................
      C=-2.*GAM+1.                                                      
C     ..................................................................
      C=GAMCAR(C)                                                       
C     ..................................................................
      D=-GAM+1.+DCMPLX(0.D0,Y)                                          
C     ..................................................................
      D=GAMCAR(D)                                                       
C     ..................................................................
      A=COSETA*EPLUS*A/B+SINETA*EMIN*C/D                                
C     ..................................................................
      H=CDABS(A)                                                        
C     ..................................................................
      H=1./H                                                            
      ANORM=F/G*H                                                       
      GOTO 11                                                           
   10 CONTINUE                                                          
      A=2.*GAM+1.                                                       
C     ..................................................................
      A=GAMCAR(A)                                                       
C     ..................................................................
      B=GAM+1.+DCMPLX(0.D0,Y)                                           
C     ..................................................................
      B=GAMCAR(B)                                                       
C     ..................................................................
      C=-2.*GAM+1.                                                      
C     ..................................................................
      C=GAMCAR(C)                                                       
C     ..................................................................
      D=-GAM+1.+DCMPLX(0.D0,Y)                                          
C     ..................................................................
      D=GAMCAR(D)                                                       
C     ..................................................................
      AA=DCMPLX(0.D0,ETA)                                               
      AA=CDEXP(AA)                                                      
      ALPHA=1.D0/137.03602D0                                            
      ZA=Z*ALPHA                                                        
      GS=DABS(DFLOAT(KAPPA))                                            
      GG=DSQRT(ZA*ZA-GS*GS)                                             
      BB=DCMPLX(0.D0,-ETA)-PI*GG                                        
      BB=CDEXP(BB)                                                      
      CC=(GAM-DCMPLX(0.D0,Y))/(DFLOAT(KAPPA)+DCMPLX(0.D0,Y)/E)          
      A=AA*A/B+BB*CC*C/D                                                
      H=1.D0/CDABS(A)                                                   
      ANORM=F/G*H                                                       
11    AIM=DIMAG(A)                                                      
      ARE=DREAL(A)                                                      
      DELTA=DATAN(AIM/ARE)+Y*(DLOG(DABS(Y))-1.)-PI/4.                   
      IF(AIM.LT.0.AND.ARE.LT.0) DELTA=DELTA-PI                          
      IF(AIM.GT.0.AND.ARE.LT.0) DELTA=DELTA+PI                          
      RETURN                                                            
      END                                                               
C     ******************************************************************
C     ******************************************************************
      SUBROUTINE TANETA(R0,E,GAM,KAPPA,Y,Z,ETA,EPLUS,EMIN,SINETA,COSETA)
C     ******************************************************************
C     ******************************************************************
C     THIS PROGRAM EVALUATES THE MATCHING PHASE ETA WHICH DESCRIBES     
C     THE EFFECT OF THE EXTENDED NUCLEUS.                               
C     SEE EQ 5.11 AND 6.9 OF REF 1.                                     
C     REF 10. BRONSTEIN- SEMENDJAJEW, TASCHENBUCH DER MATHEMATIK,       
C             VERLAG HARRI DEUTSCH, FRANKFURT AM MAIN 1968              
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 A,BBB,X,EPLUS,EMIN,GAM,WI1,WI2,S1,S2                                 
      COMMON /ZELLE/A1,A2
C	  2021-04-20 Common /Witz/ added                                                  
      COMMON /WITZ/WI1,WI2,S1,S2                                                     
C      PI=3.1415926535897D0
      COMWAV=386.1592D0                                                 
      N=1                                                               
      P=DSQRT(E*E-1.)                                                   
      RS=R0/COMWAV                                                      
      A=DCMPLX(0.D0,2.*P*RS)                                            
      C=1.D0                                                              
      IF (E.LT.0.) C=-C                                                 
      I=1                                                               
      T=DSQRT((E-1.)/(E+1.))                                            
C     ..................................................................
      X=BBB(I,GAM,Y,E,KAPPA,EPLUS,EMIN,A,Z)                             
C     ..................................................................
C	2021 April 20, Write for test purpose only
      Write (14,1613) I,KAPPA,GAM,Y,E,EPLUS,EMIN,A,X,Z
1613  Format (1x,'I,KAPPA,GAM,Y,E,EPLUS,EMIN,A,,X,Z',2I3,1P4D20.10/
     &2(15x,1P5D20.10/))
      F1=DREAL(X)                                                       
      F2=DIMAG(X)
      WRITE (14,6120) F1,F2
6120  FORMAT (1x,'F1 and F2/TANETA:',1P2D20.10)                                                             
      I=-1                                                              
C     ..................................................................
      X=BBB(I,GAM,Y,E,KAPPA,EPLUS,EMIN,A,Z)                             
C     ..................................................................
C	2021 April 20, Write for test purpose only
      Write (14,1614) I,KAPPA,GAM,Y,E,EPLUS,EMIN,A,X,Z
1614  Format (1x,'I,KAPPA,GAM,Y,E,EPLUS,EMIN,A,X,Z',2I3,1P4D20.10/
     &2(15x,1P5D20.10/))
      F3=DREAL(X)                                                       
      F4=DIMAG(X)
C     2021, April 20: In next line 1.D-6 changed to 1.D-10                                                             
      IF (DABS(DIMAG(GAM)).GT.1.D-6) GO TO 10                           
      TEN=-(A2*F1+C*T*A1*F2)/(A2*F3+C*T*A1*F4)                          
C     IF (DABS(TEN).GT.1.D10) GO TO 20                                   
      ETA=DATAN(TEN)                                                    
      SINETA=DSIN(ETA)                                                  
      COSETA=DCOS(ETA)                                                  
C     SEE FOR EXAMPLE P. 156 OF REF 10.
      Write (14,1612) ETA,A1,A2,F1,F2,F3,F4,C,T,TEN      
      RETURN                                                            
   10 CONTINUE                                                          
      TEN=(A2*F1+C*T*A1*F2)/(A2*F4-C*T*A1*F3)                           
      IF (DABS(TEN).GT.1.D10) GO TO 20                                  
      ETA=DATAN(TEN)                                                    
      SINETA=DSIN(ETA)                                                  
      COSETA=DCOS(ETA)
      Write (14,1611) ETA,A1,A2,F1,F2,F3,F4,C,T 
      WRITE (14,1621)  ETA,A1,A2,F1,F2,F3,F4,C,T
1621  FORMAT (3(1X,1P3D30.20))                                                           
      RETURN                                                            
   20 CONTINUE
C	2021 April 20, Write for test purpose only
      Write (14,1610) ETA,A1,A2,F1,F2,F3,F4,C,T,TEN
1610  Format (1x,'A20:ETA A1 A2 F1-4 C T TEN',1P5D20.10/15x,1P5D20.10)
1611  Format (1x,'B20:ETA A1 A2 F1-4 C T TEN',1P5D20.10/15x,1P5D20.10)
1612  Format (1x,'B10:ETA A1 A2 F1-4 C T TEN',1P5D20.10/15x,1P5D20.10)
      F=TEN/DABS(TEN)                                                   
      COSETA=0.D0                                                         
      SINETA=F                                                          
      RETURN                                                            
      END                                                               
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION BBB(I,GAM,Y,E,KAPPA,EPLUS,EMIN,X,Z)                      
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     SEE EQ 5.12 AND 6.9 OF REF 1                                      
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16  WHITT,XI,XM,X,BBB,S,EPLUS,EMIN,CDSQRT,S1,S2,GAM,S4    
      COMPLEX*16 WI1,WI2 
C	  2021-04-20 Common /Witz/ extended by S1,S2                                                     
      COMMON /WITZ/WI1,WI2,S1,S2                                              
      PI=3.1415926535897D0 
      XI=DCMPLX(0.D0,Y)+0.5                                             
      XI=-XI                                                            
      IF (DABS(DIMAG(GAM)).GT.1.D-6) GO TO 10                           
      XM=DFLOAT(I)*GAM                                                  
C     ..................................................................
      S=WHITT(XI,XM,X)                                                  
C     ..................................................................
      IF (I.EQ.1) WI1=S                                                 
      IF (I.EQ.-1) WI2=S                                                
      XI=EPLUS                                                          
      IF (I.LT.0) XI=EMIN                                               
      BBB=XI*S/CDSQRT(X)
C	  2021, April 20 Next 2 lines only for test purposes      
      WRITE (14,1600) XI,S,X,BBB
1600  FORMAT (1x,'In BBB bef.10,XI,S,X,BBB:',2X,1p4d20.10/17X,1p4d20.10)                                                      
      GO TO 40                                                          
   10 CONTINUE                                                          
      IF (I.EQ.-1) GOTO 15                                              
      XM=GAM
C	  2021, April 20 Next 2 lines only for test purposes       
      WRITE (14,1602) I,XI,XM,X
1602  FORMAT (1x,'In BBB WHIT-S1,I,XI,XM,X:',2X,I3,1p4d22.14/2(31X,
     &1p4d22.14/))         
C     ..................................................................
      S1=WHITT(XI,XM,X)                                                 
C     ..................................................................
      WRITE (14,*) S1
      XM=-XM
C	  2021, April 20 Next 2 lines only for test purposes       
      WRITE (14,1603) I,XI,XM,X
1603  FORMAT (1x,'In BBB WHIT-S2,I,XI,XM,X:',2X,I3,1p4d22.14/2(31X,
     &1p4d22.14/))                                                                   
C     ..................................................................
      S2=WHITT(XI,XM,X)                                                 
C     ..................................................................
      WRITE (14,*) S2
      WI1=S1                                                            
      WI2=S2                                                            
15    CONTINUE                                                          
      ALPHA=1.D0/137.03602D0                                            
      ZA=Z*ALPHA                                                        
      GS=DABS(DFLOAT(KAPPA))                                            
      GG=DSQRT(ZA*ZA-GS*GS)                                             
      S3=DFLOAT(I)*DEXP(-PI*GG)                                         
      S4=(GAM-DCMPLX(0.D0,Y))/(DFLOAT(KAPPA)+DCMPLX(0.D0,Y)/E)          
      BBB=(S1+S3*S4*S2)/CDSQRT(X)
C	  2021, April 20 Next 2 lines only for test purposes       
      WRITE (14,1601) X,S1,S2,S3,S4,BBB
1601  FORMAT (1x,'In BBB bef.40,X,S1-4,BBB:',2X,1p4d20.10/2(17X,
     &1p4d20.10/))                                              
   40 RETURN                                                            
      END                                                               
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION ALP(I,GAM,Y,E,KAPPA)                                     
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     SEE EQU. 5.4A AND 5.4B OF REF1.                                   
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 ALP,A,B,C,D,CDEXP,GAM                                  
      PI=3.1415926535897D0                                              
      A=DCMPLX(DFLOAT(KAPPA),0.D0)+DCMPLX(0.D0,Y)/E                     
      A=1./A                                                            
      G=-PI*DREAL(GAM)                                                  
      IF (I.LT.0) G=-G                                                  
      B=DCMPLX(0.D0,G)                                                  
      C=CDEXP(B)                                                        
      D=-(GAM+DCMPLX(0.D0,Y))                                           
      IF (I.LT.0) D=GAM-DCMPLX(0.D0,Y)                                  
      ALP=D*A*C                                                         
      RETURN                                                            
      END                                                               
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION WHITT(XI,XM,RM)                                          
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     THIS PROGRAM EVALUATES THE WHITTAKER FUNCTION FOR                 
C     COMPLEX ARGUMENTS.                                                
C     SEE EQ. 3.11 OF REF 1.                                            
C     THIS PROGRAMM HAS BEEN TESTED WITH THE HELP OF THE RECURRENCE     
C     RELATIONS (13.4.28) AND (13.4.29) OF REF. 6                       
C     THE ACCURACY IS BETTER THAN 7 PLACES FOR ABS (Z) LE 20.           
C     REF 6   ABRAMOWITZ, SEGUN - HANDBOOK OF MATHEMATICAL FUNCTIONS.   
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 CHYPX,GAMCAR,WHITT                                     
      COMPLEX*16 S,A,B,RM,XI,XM,CDLOG,CDEXP                             
      COMPLEX*16 GAMBA,GAMB,GAMA                                        
      COMMON/ADLER/GAMBA,GAMB,GAMA                                      
      A=0.5-XI+XM                                                       
      B=1.+2.*XM                                                        
      IF (CDABS(RM).LT.20.D0) GO TO 10                                  
      GAMBA=0.D0                                                        
      IF (CDABS(GAMBA).GT.1.D-35) GO TO 10
      WRITE (14,1630)A,B
1630  FORMAT (1X,'WHITT: A,B before call',1X,1P4D22.14)
C     ..................................................................
      GAMBA=GAMCAR(B-A)                                                 
      GAMB=GAMCAR(B)                                                    
      GAMA=GAMCAR(A)
      WRITE (14,1620)GAMBA,GAMB,GAMA,A,B
1620  FORMAT (1X,'WHITT: GAMBA,GAMB,GAMA,A,B',1X,1P4D22.14/
     &2(16X,1P4D22.14/))                                                          
C     ..................................................................
   10 CONTINUE                                                          
C     ..................................................................
      S=CHYPX(A,B,RM)                                                   
C     ..................................................................
      A=(XM+0.5)*CDLOG(RM)-0.5*RM
      WRITE (14,1600) XM,RM,A,S
1600  FORMAT (1X,'In WHITT: XM,RM,A,S ',1X,1P4D22.14/22X,1P4D22.14)                                             
      WHITT=CDEXP(A)*S                                                  
      RETURN                                                            
      END                                                               
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION GAMCAR(CMX)                                              
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     THIS SUBPROGRAM EVALUATES THE GAMMA FUNCTION FOR ANY              
C     COMPLEX ARGUMENT BY USE OF AN ASYMPTOTIC EXPANSION.               
C     THIS PROGRAM HAS BEEN CHECKED FOR 1. LE RE(Z) LE 2. AND           
C     0.1 LE IM(Z) LE 10.0.                                             
C     ERROR IS LESS THAN 10**-6                                         
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 GAMCAR                                                 
      COMPLEX*16 CMX,COMX,GAMX,E1,E11,TX,GAM1,CDLOG,CDEXP,SUMME,FAK       
C      PI=3.1415926535897D0                                              
      SUMME=CMX                                                           
      FAK=(1.,0.)                                                       
      IF (DREAL(CMX).LE.0.) GO TO 5                                     
   12 CONTINUE                                                          
      COMX=CMX-1.                                                       
      E1=COMX+.5                                                        
      TX=COMX+5.5                                                       
      E11=E1*CDLOG(TX)                                                  
      GAM1=CDEXP(E11-TX)*2.50662827465D0                                
      GAMX=GAM1*(1.+76.18009173D0/(COMX+1.)-86.50532033D0/(COMX+2.)     
     &+24.01409822D0/(COMX+3.)-1.231739516D0/(COMX+4.)+.12085003D-2/(COM
     &X+5.)-.536382D-5/(COMX+6.))                                       
      GAMCAR=GAMX                                                       
      GAMCAR=GAMCAR*FAK                                                 
      GO TO 2                                                           
    5 CONTINUE                                                          
      J=0                                                               
    6 CMX=CMX+1.                                                        
      J=J+1                                                             
      T=DREAL(CMX)                                                      
      IF (T.LE.0) GO TO 6                                               
      DO 11 K=1,J                                                       
      FAK=FAK*1./(SUMME+DFLOAT(K)-1.)                                     
   11 CONTINUE                                                          
      GO TO 12                                                          
    2 CMX=SUMME                                                           
      RETURN                                                            
      END                                                               
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION CHYPX(A,B,Z)                                             
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     THIS SUBPROGRAM EVALUATES THE CONFUENT HYPERGEOMETRIC FUNCTION    
C     FOR A COMPLEX ARGUMENT AND COMPLEX PARAMETERS.                    
C     NUP IS THE NUMBER OF TERMS.                                       
C     EPS DETERMINES THE ACCURACY.                                      
C     EPS IS THE PART OF THE LAST TERM ON THE WHOLE SUM.                
C     THIS PROGRAM HAS BEEN TESTED FOR  -0.9 LE RE(A) LE 0.9 AND FOR    
C     0.1 LE RE(B) LE 0.9 AND FOR 0.1 LE RE(Z) LE 10.                   
C     THE ACCURACY IS BETTER THAN 6 PLACES.                             
C     THIS PROGRAM NEEDS OVER COMMON THE INFORMATION WHAT IS GAMMA(B-A),
C     GAMMA(B) AND GAMMA(A).                                            
C     FOR THE SAKE OF ACCURACY THE POWER SERIES EXPANSION IS CALCULATED 
C     IN DOUBLE PRECISION ACCURACY FOR SPECIAL VALUES OF Z.             
C     THE ACCURACY OF THE ASYMPTOTIC IS ONLY BETTER THAN 2 PLACES.      
C     THE ASYMPTOTIC HAS TO BE CHANGED TO DOUBLE PRECISION ACCURACY.    
C     NEW VERSION   MAI 1976.                                           
C     FOR A = (3,1.5) AND B = (2,0.5) THE ACCURACY WAS FOUND TO BE      
C     BETTER THAN 12 PLACES FOR CDABS(Z).LT.20 AND 8 PLACES FOR         
C     20.LT.CDABS(Z).LT.50.                                             
C     THE ACCURACY IN THE ASYMPTOTIC REGION IS DETERMINED BY THE        
C     GAMMA FUNCTION.                                                   
      IMPLICIT REAL*8 (A-H,O-Z)
C     2021, April 20, GAMCAR removed in next line                                               
      COMPLEX*16 CHYPX                                           
      COMPLEX*16 A,B,Z,S,TRM,AT,BT,ONE,EN,S1,S2,T1,T2,CT,DT,CDLOG,CDEXP 
      COMPLEX*16 GAMBA,GAMA,GAMB                                        
      COMPLEX*16 G1,G2,G3,G4                                            
      COMMON/ADLER/GAMBA,GAMB,GAMA                                      
      DATA EPS,ONE/1.D-12,(1.,0.)/ 
      WRITE (14,1600) A,B,Z,GAMBA,GAMB,GAMA
1600  FORMAT (1X,'In Chypx, A,B,Z,GAMBA,GAMB,GAMA ',1X,
     &1P4D22.14/2(18X,1P4D22.14/))                                           
      IF (CDABS(Z).GE.20.) GO TO 50                                     
C     EVALUATION OF A POWER SERIES EXPANSION FOR ABS(Z) .LE. 20.        
      NUP=100                                                           
      S=ONE                                                             
      TRM=ONE                                                           
      AT=A-ONE                                                          
      BT=B-ONE                                                          
      EN=(0.,0.)                                                        
      DO 10 N=1,NUP                                                     
      AT=AT+ONE                                                         
      BT=BT+ONE                                                         
      EN=EN+ONE                                                         
      TRM=TRM/EN*AT/BT*Z                                                
      S=S+TRM                                                           
      IF (CDABS(TRM).LE.EPS*CDABS(S)) GO TO 20                          
   10 CONTINUE                                                          
      WRITE (6,100) A,B,Z                                               
  100 FORMAT(1H1,'CHYPX DOES NOT CONVERGE  '/,1H0,3(2D15.6,5X))         
   20 CHYPX=S
      WRITE (14,1610) S
1610  FORMAT (1X,'In CHYPX, S: ',1PD24.16)                                                                 
      RETURN                                                            
C     EVALUATION OF AN ASYMPTOTIC EXPANSION FOR ABS(Z) .GT. 20.         
   50 CONTINUE                                                          
      EPSA=1.D-6                                                       
      II=1                                                              
      S1=ONE                                                            
      S2=ONE                                                            
      AT=A-ONE                                                          
      BT=1.+A-B-ONE                                                     
      CT=B-A-ONE                                                        
      DT=1.-A-ONE                                                       
      G1=A                                                              
      G2=1.+A-B                                                         
      G3=B-A                                                            
      G4=1.-A                                                           
      IMAX=15                                                           
      DO 70 I=1,IMAX                                                    
      IF (I.EQ.1) GO TO 40                                              
      AT=AT*(G1+I-1)                                                    
      BT=BT*(G2+I-1)                                                    
      CT=CT*(G3+I-1)                                                    
      DT=DT*(G4+I-1)                                                    
   40 CONTINUE                                                          
      IF (I.GT.1) GO TO 30                                              
      AT=AT+ONE                                                         
      BT=BT+ONE                                                         
      CT=CT+ONE                                                         
      DT=DT+ONE                                                         
   30 CONTINUE                                                          
      T1=AT*BT*(-Z)**(-I)                                               
      T2=CT*DT*Z**(-I)                                                  
      II=I*II                                                           
      T1=T1/DFLOAT(II)                                                  
      T2=T2/DFLOAT(II)                                                  
      S1=S1+T1                                                          
      S2=S2+T2                                                          
C     W=CDABS(T1)+CDABS(T2)                                             
C     WW=CDABS(S1)+CDABS(S2)                                            
C     IF (W.LE.EPS*WW.AND.I.GT.3) GO TO 7                               
      IF(CDABS(T1).LE.EPSA*CDABS(S1).AND.CDABS(T2).LE.EPSA*CDABS(S2)    
     & .AND. I.GT.4) GOTO 7                                             
   70 CONTINUE                                                          
      WRITE (6,110) T1,S1,T2,S2                                         
  110 FORMAT (' NO CONVERGENCE IN CHYPX ',1P4D16.8/,1P4D16.8)           
      WRITE (6,111) A,B,Z                                               
111   FORMAT (1P6D16.8)                                                 
    7 CONTINUE                                                          
      S=CDLOG(Z)                                                        
      E=DIMAG(Z)                                                        
      IF((DABS(E)-0.D0).LT.1.D-20) GO TO 90                                              
      E=E/DABS(E)*3.1415926535897D0                                     
      GO TO 80                                                          
   90 E=3.1415926535897D0                                               
   80 CONTINUE                                                          
      T1=CDEXP(A*(DCMPLX(0.D00,E)-S))/GAMBA                             
      T2=CDEXP(Z+(A-B)*S)/GAMA                                          
      S=(S1*T1+S2*T2)*GAMB                                              
      GO TO 20                                                          
      END

Probably there is a simple reason for the different results, but I don't recognize it yet.

Best regards Theo

The source code is located unter the following link: https://www.magentacloud.de/share/gnh315jfie

Theo, you posted only the source file and one data file to MagentaCloud. The file OPEN statements are inconsistent with the usage of the files. I corrected this by changing 'OLD' to 'REPLACE' in the OPEN statements for units 6 and 14.

As you may have noticed, there is a limit of about 100 lines per post, so long program sources may lose major portions when posted inline. For short code extracts, please select the source lines with the mouse and click on the Code button. The lines will then be formatted nicely.

The function GAMCAR takes one argument, CMX. This variable is both read and updated in the function. This is usually a cause of trouble, because the normal expected behaviour of a function is to return a function value without changing any of its arguments. In your program, this function GAMCAR is also referenced with an expression as an argument (on line 419). With most modern compilers, you are not allowed to place a value into an expression -- you can only place values into variables. I corrected this by changing the name of the dummy argument of GAMCAR to CMXARG, and setting CMX=CMXARG at the beginning of the function.

After these changes, I compiled and ran the program. The program printed 'STOP: 111' and stopped, since the input data provides '0' for ITEST.

The results in the output file generated, FT14WF.DAT, was nearly the same with two different compilers: FTN95 8.71, 64-bit, and Lahey LF71, 32-bit. There were small differences in the least significant digits, but that is not unusual with floating point calculations.

I then compared the output of 32-bit FTN95 and 64-bit FTN95. Again, the results agreed to at least 12 significant digits.

What next? What did you do to produce the different results that you did not feel were correct?

mecej4, first of all thank you for your second timely response and reviewing the source code! I added the two files delivering different results when compiling the original source code with 32 Bit, resp. 64 Bit compiler (FT14WF-32BIT.DAT and FT14WF-64BIT.DAT). The differences in GAMBA are as follows (first printout):

(1.47646338812876D-02,-1.63783687568201D-02) 32-Bit (7.38231694064379D-03,-8.18918437841007D-03) 64-Bit

This answers your last question. By error propagation from this subroutine I got faulty results in the matrix elements.

I was'nt aware that passing arguments forth and back via the subroutine parameter list is a problem for current compilers, since that was common use (and allowed) in IBM- Fortran (about 40 years ago..). So thanks a lot for this valuable hint! Same with line 419: only argument, no expression in function call.

After changing this according to your recommendation, I got agreeing results - within 11-12 significant digits (c.f. file FT14WF-64Bit-new.DAT) between FTN95 32- and 64- Bit runs. Stop 111 is a regular RETURN in the subroutine. So the FTN95 32Bit compiler obviously tolerates the mentioned coding deficiencies, wheras the 64Bit compiler doesn't.

I guess I wouldn't have figured it out without your help ! What next? I will search the complete source code concerning modification of passed arguments in subroutine calls (hours of enjoyment;-) plus possible expressions in function argument lists, correct and than check the final result. Thank you very much for your help, mecej4!

Kind regards Theo

Theo, I am gratified to see you catching up so quickly, and I have good news for you. Silverfrost FTN95 is very good at catching bugs in Fortran programs at compile time (i.e., before even running the program), but there are many bugs that can only be caught at run time, and FTN95 provides facilities to catch those as well.

Try the /debug, /check and /undef options. When your code contains more Fortran 95 /2003 features, you may also find /checkmate useful.

Here is a buggy but short program that may interest you. Try running it with different options, 32 and 64 bit, and you may find some of the outcomes quite puzzling and surprising. In this connection, see also https://medium.com/@patkerr/should-the-value-of-pi-change-f18f41599919 .

program OneIsTwo
implicit none
integer One, Two
One = 1
Two = 2
call Sub(One, Two)
Print *, 'One = ',One, 'Two = ',Two
call Sub(1, 2)
Print *, 'Later, One = ',1, ' and Two = ',2
end program

subroutine sub(A,B)
implicit none
integer A,B,C
C=A
A=B
B=C
return
end subroutine

What would the old IBM computer output for this program (after reformatting the source to Fortran 4 or Fortran 77 conventions to make it compatible)?

Theo, is it just serendipity that the answers differ by a factor of 2?

Bill

Well the first result is as expected in good old IBM Fortran (One = 2 and Two =1) and the second yields an error - as expected.

Using checkmate indeed delivers more information: Runtime error from program:d:\fortran\me-dive\fixedformat1.exe Run-time Error *** Error 14, Attempt to alter an actual argument that is a constant, an expression, an alias, an INTENT(IN) argument, or a DO variable

SUB - in file fixedformat1.for at line 16 [+009f]

ONEISTWO - in file fixedformat1.for at line 8 [+01af] I have used the /debug, /check and /undef parameters and eliminated tons of undefined variables in the FT77 pgms (IBM Fortran used to initialize every variable nicely with 0), but that didn't give me a hint to the problem with the argument list in the subroutine calls.

Thanks again for your help - I have truncated the source code in magenta cloud to the relevant functions.

Kind regards Theo

Good point Bill (eagle-eye)! I was so busy chasing the occurence, where the differences between the two compiler versions occur, that I didn't even notice that the results differ by a factor 2.

The factor 2 indeed seems to be serendipity, since (to my knowledge) there is no recursion relation for the Gamma functions in the form Gamma(B-A)=Gamma(B)/N or Gamma(A)/N [in case GAMCAR(B-A) would be interpreted either as GAMCAR(B) or GAMCAR(A)] for complex arguments A,B and integer N. The difference origins from the faulty argument passing at FUNCTION GAMCAR(CMX) via CMX as described by mecej4. It seems that this leads to returning a faulty function value by a factor 2.

Best regards Theo

Theo and Bill,

That factor of 2 that you noticed is just happenstance. When a Fortran program is in violation of the standard, most often the compiler and RTL may do anything, including

1. Refuse to compile the program.
2. Produce an EXE that runs a bit, gives an error message and stops.
3. Produce an EXE that hangs up, and you may have to use the OS tools to abort it.
4. Produce incorrect results. If they are only slightly incorrect, the user may not even notice that something is amiss.
5. Produce correct results.
6. Start WW-III.

and so on.

The Fortran standard simply uses the phrase 'behavior is undefined' to cover all these possibilities.

Theo's IBM 3033 code contains a surprisingly good and efficient function subprogram for computing the Gamma function of a complex argument within certain ranges of the argument. I found it a striking instance where the facilities of modern Fortran can enable the programmer to express the algorithm in a compact, efficient and readable form, in contrast to the 'spaghetti' style of the original.

!     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      complex(dble_complex) function gamcar (cmx)
      USE kind_param, ONLY:  dble_complex, double
!     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!     THIS SUBPROGRAM EVALUATES THE GAMMA FUNCTION FOR ANY
!     COMPLEX ARGUMENT BY USE OF AN ASYMPTOTIC EXPANSION.
!     THIS PROGRAM HAS BEEN CHECKED FOR 1. LE RE(Z) LE 2. AND
!     0.1 LE IM(Z) LE 10.0.
!     ERROR IS LESS THAN 10**-6
!**************************************************************
      implicit none
      integer :: j, k
      real(double) :: t
      complex(dble_complex) :: cmx, comx, tx, gam1
!-----------------------------------------------
      k = nint(dble(cmx))
      comx = cmx + k
      tx   = comx + 5.5
      gam1 = exp((comx + 0.5)*log(tx) - tx)*2.50662827465D0
      gamcar = (1. + 76.18009173D0/(comx + 1.) - 86.50532033D0/(comx + 2.) &
                   + 24.01409822D0/(comx + 3.) - 1.231739516D0/(comx + 4.) &
                   + 0.12085003D-2/(comx + 5.) -   0.536382D-5/(comx + 6.) &
               )*gam1/product([(cmx+j, j=0,k)])
      return
      end function gamcar

Contrast this with the original (with potentially illegal modification of function argument), continued in the next post in the thread (forum limit)

[CONTINUED]

[CONTINUED]

C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      FUNCTION GAMCAR(CMX)                                              
C     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C     THIS SUBPROGRAM EVALUATES THE GAMMA FUNCTION FOR ANY              
C     COMPLEX ARGUMENT BY USE OF AN ASYMPTOTIC EXPANSION.               
C     THIS PROGRAM HAS BEEN CHECKED FOR 1. LE RE(Z) LE 2. AND           
C     0.1 LE IM(Z) LE 10.0.                                             
C     ERROR IS LESS THAN 10**-6                                         
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      COMPLEX*16 GAMCAR                                                 
      COMPLEX*16 CMX,COMX,GAMX,E1,E11,TX,GAM1,SUMME,FAK       
      SUMME=CMX                                                           
      FAK=(1.,0.)                                                       
      IF (DBLE(CMX).LE.0.) GO TO 5                                     
   12 CONTINUE                                                          
      COMX=CMX-1.                                                       
      E1=COMX+.5                                                        
      TX=COMX+5.5                                                       
      E11=E1*LOG(TX)                                                  
      GAM1=EXP(E11-TX)*2.50662827465D0                                
      GAMX=GAM1*(1.+76.18009173D0/(COMX+1.)-86.50532033D0/(COMX+2.)     
     &+24.01409822D0/(COMX+3.)-1.231739516D0/(COMX+4.)+.12085003D-2/(COM
     &X+5.)-.536382D-5/(COMX+6.))                                       
      GAMCAR=GAMX                                                       
      GAMCAR=GAMCAR*FAK                                                 
      GO TO 2                                                           
    5 CONTINUE                                                          
      J=0                                                               
    6 CMX=CMX+1.                                                        
      J=J+1                                                             
      T=DBLE(CMX)                                                      
      IF (T.LE.0) GO TO 6                                               
      DO 11 K=1,J                                                       
      FAK=FAK*1./(SUMME+K-1.)                                     
   11 CONTINUE                                                          
      GO TO 12                                                          
    2 CMX=SUMME
      RETURN                                                            
      END

mecej4, thank you for explaining what happens when a Fortran program is in violation of standards.

The specific code for evaluating the Gamma function is a module which was not coded by myself but by one of the members of the Theoretical Physics group in Frankfurt (larger programs were jointly coded by the members of the Atomic Physics group). At that time I would have code it simlarly within FT77 Standard.

Your recoding is fascinating short and elegant ! (maybe I should have a deeper look into FTN95 coding standards). I will definitly test and use this.

Best regards Theo

The IMSL library has the generic function GAMMA -- not surprising --, so I compared your results with what IMSL gives. If you have IMSL or have access to a computer with IMSL, you can test this code:

PROGRAM TGAM
USE GAMMA_INT
IMPLICIT NONE
integer, parameter :: double = kind(0d0)
!
complex(double) :: ab,gab
complex sab,gimsl
ab = (0d0,2.69747d0)
sab = ab
gab = gamcar(ab)
gimsl = gamma(sab)
print 10,ab,gab,gimsl
10 format(/,' Argument : ',2ES13.5,/,' GamCar   : ',2ES13.5,/, &
            ' IMSL     : ',2ES13.5)

CONTAINS
!     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
      complex(double) function gamcar (cmx)
!     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
!     THIS SUBPROGRAM EVALUATES THE GAMMA FUNCTION FOR ANY
!     COMPLEX ARGUMENT BY USE OF AN ASYMPTOTIC EXPANSION.
!     THIS PROGRAM HAS BEEN CHECKED FOR 1. LE RE(Z) LE 2. AND
!     0.1 LE IM(Z) LE 10.0.
!     ERROR IS LESS THAN 10**-6
!**************************************************************
      implicit none
      integer :: j, k
      real(double) :: t, G0 = 2.50662827465D0
      real(double) :: G(6) = [76.18009173D0, -86.50532033D0, 24.01409822D0, &
                              -1.231739516D0, 0.12085003D-2, -0.536382D-5]
      complex(double) :: cmx, comx, tx, gam1
!-----------------------------------------------
      k = nint(dble(cmx))
      comx = cmx + k
      tx   = comx + 5.5
      gam1 = exp((comx + 0.5)*log(tx) - tx)*G0
      gamcar = (1. + sum([(G(j)/(comx + j), j = 1, 6)])) &
              *gam1/product([(cmx+j, j=0,k)])
      return
      end function gamcar
end program

Output:

T:\LANG\THEO\F90>f95 %f90flags% gamcar.f90 %link_fnl%

T:\LANG\THEO\F90>gamcar

 Argument :   0.00000E+00  2.69747E+00
 GamCar   :   1.47647E-02 -1.63785E-02
 IMSL     :   1.47647E-02 -1.63785E-02

For the argument that I used, at least, it appears that COMPLEX*16 may not be needed.

The NSWC Fortran subroutine library contains the routine CGAMMA, which appears to be quite suitable as a replacement for your GAMCAR:

https://wp.csiro.au/alanmiller/NSWC/cgamma.f90

For NSWC Fortran routines + FTN95 see:

https://gitlab.com/silverfrost/nswc

Thanks, Kenneth. I downloaded the DLLs and the test program source from the Silverfrost Gitlab page. I built the 32-bit exe, and it ran fine, displaying your examples. The 64-bit exe was also built, but when run it displayed 'The procedure entry point DGAMMA@ could not be located in the dynamic link library ...NSWC64.DLL.

Mecej4, Thanks for the feedback. Yes, I can replicate that error for the 64 bit test program using the the current FTN version, with the dlls on the Gitlab page that were prepared using FTN v8.97.2.

If I create a new 64 bit dll with the current FTN version all appears to run correctly.

Here are the dlls I created this evening: https://www.dropbox.com/scl/fi/0sloctwtvuctmfk7is4g9/nswc-FTN95-Ver.9.14.0.zip?rlkey=a1tbps9umydbje5u4cjrn5k66&st=5cc3nq9c&dl=0

Thanks for the new DLLs, Kenneth. Everything is fine now.

Not entirely fine - if you use Plato.

The Plato help says 'When there is no project open, you can view the Project Explorer window in order to add references that are used when compiling and linking the current file.'

Thus a simple way to compile and run the example file from within Plato is:

1. Open the file nswctest.f95 in Plato

2. In the VIEW tab, ensure that the following are selected: Project Explorer, Standard Tool Bar, and Build Tool Bar

3. Select either the WIN32 or X64 compiler in the Standard Tool Bar.

4. In the Project Explorer window, right click on References, and select Add References.
   This will open a dialogue which allows you to select:

   nswc32.dll for the 32 bit compiler WIN32

   or

   nswc64.dll for the 64 bit compiler X64

5. Use the Compile button in the Build Tool Bar to compile nswctest.f95

6. Use the Build button in the Build Tool Bar to compile nswctest.f95 and link the object file with the dll selected in the References window.

7. Use the Start button in the Build Tool Bar to execute the program.

Now this definitely worked back in 2022, while the present Plato does not allow me to add the required DLL as a reference. Something has changed, or I have forgotten a step.

Sorry about this. There is a regression in the current Plato v6.0.4.

Unfortunately the pending full release will use v6.0.4.

Users who need this feature can send me a PM for a link to v6.0.5 whilst noting that the next full release will revert to v6.0.4.