【发布时间】：2014-02-13 17:05:58
【问题描述】：
我尝试运行我的 f90 代码。它是旧 f77 代码的转换版本。
当我尝试使用不同的编译器（IFORT、GFORTRAN）编译它时，我有两个不同的结果： 用两个编译器自己运行程序，看看，例如 GNUplot，情节：
plot 'orbitm.txt' u 1:2
使用两个编译器，绘图的输出非常不同！
我认为这种不同（也）依赖于COMMON 命令，我尝试用MODULE 替换它，但我发现了一些错误。
我按照 cmets 的建议对我的代码进行了一些更改：
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à 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='orbitm.txt')
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,'(4g12.5)') 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
现在 IFORT 和 GFORTRAN 之间的编译差异很小但不等于零。 是否可以通过其他方式对代码进行更多改进，例如：更改调用或函数、拆分子程序或引入其他模块？
非常感谢！
【问题讨论】：
-
在两个编译器上运行 F77 版本有什么不同吗？
-
在我看来，这两个文件都像是从
(0,0)到大约(48000,13000)的直线。我应该看到什么样的“非常不同”的东西？ -
我当然不会运行你的代码并仔细研究它的内部来发现问题。你自己做了什么来找到它？如果你真的认为是
common块，为什么不简单地用每个范围内的参数替换它，它只包含两个变量？ -
这里的
common看起来很安全，因为它出现在主程序中。采用添加调试写入的旧方法并跟踪实际问题所在，而不是根据猜测进行更改。 -
如果您要解决微妙的数字问题，我建议您明确指定文字的精度，例如。
mu=0.001D0
标签： fortran fortran90 gfortran intel-fortran fortran77