forums.silverfrost.com
Welcome to the Silverfrost forums

Author Message
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Wed Apr 21, 2021 10:27 pm    Post subject: 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: [code:1:aabccb7571]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
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Wed Apr 21, 2021 10:41 pm    Post subject: Source code The source code is located unter the following link: https://www.magentacloud.de/share/gnh315jfie
mecej4

Joined: 31 Oct 2006
Posts: 1680

 Posted: Wed Apr 21, 2021 11:53 pm    Post subject: 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?
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Thu Apr 22, 2021 10:50 am    Post subject: 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
mecej4

Joined: 31 Oct 2006
Posts: 1680

Posted: Thu Apr 22, 2021 11:07 am    Post subject:

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 .

 Code: 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)?
wahorger

Joined: 13 Oct 2014
Posts: 1021
Location: Morrison, CO, USA

 Posted: Thu Apr 22, 2021 5:48 pm    Post subject: Theo, is it just serendipity that the answers differ by a factor of 2? Bill
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Thu Apr 22, 2021 6:00 pm    Post subject: 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
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Fri Apr 23, 2021 2:38 pm    Post subject: 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
mecej4

Joined: 31 Oct 2006
Posts: 1680

 Posted: Fri Apr 23, 2021 4:51 pm    Post subject: 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.Last edited by mecej4 on Fri Apr 23, 2021 5:16 pm; edited 1 time in total
mecej4

Joined: 31 Oct 2006
Posts: 1680

Posted: Fri Apr 23, 2021 4:58 pm    Post subject: Function GAMCAR in modern Fortran

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.

 Code: !     +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-       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]

Last edited by mecej4 on Fri Apr 23, 2021 8:32 pm; edited 1 time in total
mecej4

Joined: 31 Oct 2006
Posts: 1680

Posted: Fri Apr 23, 2021 5:00 pm    Post subject:

[CONTINUED]
 Code: 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
Theo1002

Joined: 19 Apr 2021
Posts: 15
Location: Ingelheim am Rhein, Germany

 Posted: Sat Apr 24, 2021 10:32 am    Post subject: 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
mecej4

Joined: 31 Oct 2006
Posts: 1680

Posted: Sun Apr 25, 2021 1:14 pm    Post subject:

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:

 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:

 Code: 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.
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT + 1 Hour Page 1 of 1