|
forums.silverfrost.com Welcome to the Silverfrost forums
|
View previous topic :: View next topic |
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 |
|
Back to top |
|
|
Theo1002
Joined: 19 Apr 2021 Posts: 15 Location: Ingelheim am Rhein, Germany
|
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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? |
|
Back to top |
|
|
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 |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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)? |
|
Back to top |
|
|
wahorger
Joined: 13 Oct 2014 Posts: 1217 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 |
|
Back to top |
|
|
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 |
|
Back to top |
|
|
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 |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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 |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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 |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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 |
|
|
Back to top |
|
|
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 |
|
Back to top |
|
|
mecej4
Joined: 31 Oct 2006 Posts: 1885
|
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. |
|
Back to top |
|
|
|
|
You cannot post new topics in this forum You cannot reply to topics in this forum You cannot edit your posts in this forum You cannot delete your posts in this forum You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|