F90 compiler dependent code

Question

i try to run my f90 code. It's a converted version of an old f77 code.

When i try to compile it with different compiler (IFORT, GFORTRAN) I have two different results: run the program yourself with the two compiler and see, for istance with GNUplot, the plot:

plot &#39;orbitm.txt&#39; u 1:2

With the two compiler the output of the plot is VERY different !

I suppose that this different depend (also) in such a way by the COMMON command, i try to replace it with a MODULE but i find out some errors.

I put some change to my code as recommended by the comments:

     module data
      REAL*8 :: OME = 1.D0
      REAL*8 :: MU = 0.000954
     end module


      PROGRAM MAIN
      use data
      IMPLICIT NONE
      REAL*8  :: dist0 , dt , duepi , e , e0 , ermed , errh , H , k_max_r8
      REAL*8  :: ptau , ptau0 , px , px0 , py , py0 , t , tau , tau0 , n_step_r8
      REAL*8  :: t_per, x , x0 , y , y0, k_r8 , m_per
      INTEGER :: k , k_max  , n_step

      duepi = 8.d0*DATAN(1.D0)
      !duepi = 2.d0*3.1415926535897932d0
      t_per = duepi/OME
      n_step = 1000

!         do iX0=1,6
!          do iY0=-50,50

!       x0 = (0.5d0-MU)+(0.0001*iX0) !v&#224; da 0.449 a 0.599
!                y0 = (sqrt(3.d0)/2.d0)+(0.0001*iY0) ! va da 0.868 a 0.864

      OPEN (UNIT=11,FILE=&#39;orbitm.txt&#39;) 
      x0 = 0.47
      y0 = SQRT(3.D0)/2.D0
      tau0 = 0.d0
      px0 = OME*y0
      py0 = -OME*x0
      ptau0 = 1.d0

      x = x0
      y = y0
      tau = tau0
      px = px0
      py = py0
      ptau = ptau0
      n_step_r8 = real(n_step)
      dt = t_per/n_step_r8
      e0 = H(x,y,tau,px,py,ptau)

      k_max = 1000*n_step

      k = 0
      t = 0.d0
      errh = 0.d0

!---------
!  inizio loop integrazione
!--------
      ermed = 0.d0

      DO k = 1 , k_max
         k_r8 = real(k)
         CALL SYM4(x,y,tau,px,py,ptau,dt)
         e = H(x,y,tau,px,py,ptau)
         errh = ABS(e-e0)
         t = k_r8*dt
         IF ( MOD(k,n_step).EQ.0 ) THEN
            WRITE (11,&#39;(4g12.5)&#39;) x , y , px , py
         ENDIF
      ENDDO
      k_max_r8 = real(k_max)



      DO k = 1 , k_max
         CALL SYM4(x,y,tau,px,py,ptau,-dt)
         e = H(x,y,tau,px,py,ptau)
         errh = ABS(e-e0)
         t = t - dt
      ENDDO
!           write(*,*) ix0,ermed,errh
!              enddo ! iY0
!            enddo ! iX0

!           close(11)
      END


      REAL*8 FUNCTION H(X,Y,Tau,Px,Py,Ptau)
      use data
      IMPLICIT NONE
      REAL*8 :: c , Ptau , Px , Py , r1 , r2 , s , Tau , X , Y

      c = COS(OME*Tau)
      s = SIN(OME*Tau)

      r1 = SQRT((X+MU*c)**2+(Y-MU*s)**2)
      r2 = SQRT((X-(1.d0-MU)*c)**2+(Y+s*(1.d0-MU))**2)

      H = (Px*Px)/2.D0 + (Py*Py)/2.D0 + Ptau - (1.d0-MU)/r1 - MU/r2

      END


      SUBROUTINE F(X,Y,Tau,Fx,Fy,Ftau)
      use data
      IMPLICIT NONE
      REAL*8 :: c, Ftau , Fx , Fy , r1 , r2 , s , Tau , X , Y

      c = COS(OME*Tau)
      s = SIN(OME*Tau)

      r1 = SQRT((X+MU*c)**2+(Y-MU*s)**2)
      r2 = SQRT((X-(1.d0-MU)*c)**2+(Y+s*(1.d0-MU))**2)


      Fx = -((1.d0-MU)*(X+MU*c))/r1**3 - (MU*(X-c*(1.d0-MU)))/r2**3
      Fy = -((1.d0-MU)*(Y-MU*s))/r1**3 - (MU*(Y+s*(1.d0-MU)))/r2**3

      Ftau = -( (1.d0-MU)*OME*MU*(-s*(X+MU*c)-c*(Y-MU*s)) )/r1**3.d0 - ( MU*(1.d0-MU)*OME*(s*(X-(1.d0-MU)*c)+c*(Y+(1.d0-MU)*s)) )/r2**3.d0

      END



      SUBROUTINE SYM2(X,Y,Tau,Px,Py,Ptau,Dt)
      IMPLICIT NONE
      REAL*8 :: Dt , ftau , ftaunew , fx , fxnew , fy , fynew , Ptau
      REAL*8 :: ptaunew , Px , pxnew , Py , pynew , Tau , taunew , X
      REAL*8 :: xnew , Y , ynew

      CALL F(X,Y,Tau,fx,fy,ftau)


      xnew = X + Px*Dt + fx*(Dt**2.d0)/2.D0

      ynew = Y + Py*Dt + fy*(Dt**2.d0)/2.D0

      taunew = Tau + Dt

      CALL F(xnew,ynew,taunew,fxnew,fynew,ftaunew)
      pxnew = Px + Dt*(fx+fxnew)/2.D0
      pynew = Py + Dt*(fy+fynew)/2.D0
      ptaunew = Ptau + Dt*(ftau+ftaunew)/2.D0

      X = xnew
      Y = ynew
      Tau = taunew
      Px = pxnew
      Py = pynew
      Ptau = ptaunew

      END

      SUBROUTINE SYM4(X,Y,Tau,Px,Py,Ptau,Dt)
      IMPLICIT NONE
      REAL*8 :: alpha , beta , Dt , dt1 , dt2 , Ptau , Px , Py , sq2
      REAL*8 :: Tau , X , Y
      sq2 = 2.d0**(1.D0/3.D0)
      alpha = 1.D0/(2.d0-sq2)
      beta = sq2/(2.d0-sq2)
      dt1 = Dt*alpha
      dt2 = -Dt*beta
      CALL SYM2(X,Y,Tau,Px,Py,Ptau,dt1)
      CALL SYM2(X,Y,Tau,Px,Py,Ptau,dt2)
      CALL SYM2(X,Y,Tau,Px,Py,Ptau,dt1)
      END

Now the difference in compilation between IFORT AND GFORTRAN ARE SMALL but not equal to zero. Can I improve more the code with other prescription as, for instance: changing the calling or the function, split the subroutine or introduce others modules?

Thanks a lot !

Are there any differences when you run the F77 versions on both compilers? — cup, Feb 13 '14 at 17:25
Both files look like straight lines from `(0,0)` to about `(48000,13000)` to me. What kind of "VERY different" things am I supposed to see? — Kyle Kanos, Feb 13 '14 at 17:25
I'm certainly not going to run your code and pore through its entrails to spot the problem. What have you done to find it yourself ? If you really think it is down to the `common` block why not simply replace it with parameters in every scope, it only contains two variables ? — High Performance Mark, Feb 13 '14 at 17:26
the `common` here looks safe since it appears in the main program. Take the good old approach of adding debugging writes and track down where the actual problem is instead of making changes based on guesses. — agentp, Feb 13 '14 at 19:59
OK, i change (as Mark told) the `common` with module and declare the two variables. But the result is, again, the same: ifort and gfort give two different results. So the problem can be, as mentioned by M.S.B. the difference in numerical operations that made the two compilers. — Panichi Pattumeros PapaCastoro, Feb 14 '14 at 16:22
If you are chasing down subtle numerical issues, one thing i'd advise is explicitly specifying the precision of your literals, eg. `mu=0.001D0` — agentp, Feb 15 '14 at 14:56

score 1 · Answer 1 · edited May 23 '17 at 12:05

I suggest inserting intermediate output statements and seeing where the results first diverge. Possiblities: a bug that is revealed only in one compiler. You might see a sudden, discrete change. Try to home in on where the change first occurs as a clue to the problem. Or the differences might be the result of numerical instabilities in the algorithm. In this case the differences will start very small and grow. The order of numerical operations might not be exactly the same for the executable made by two compilers and this could make a difference in some calculations.

I find a common problem in porting legacy FORTRAN 77 code to modern compilers is that that many old FORTRAN 77 programs assumed that all local variables were saveed. This was and is not guaranteed by the language standards but was a common behavior. You can restore that behavior with a compiler option. See Are local variables in Fortran 77 static or stack dynamic?

F90 compiler dependent code

1 Answers1