隐式 none 使程序不准确
Implicit none makes program inaccurate
我正在编写一个程序,该程序使用给定的子例程来计算球贝塞尔函数。我将给出 table 的子例程修改为只给出一个值的函数。然而,我意识到当我调用我的函数时我需要在我的程序中包含 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
否则它会给我一个错误的值或错误。下面我提供了一个可以正常工作的示例程序。
! n = 3, x = 2 should return ~ 6.07E-2
program hello
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
doubleprecision :: bessel, ans
WRITE(*,*)'Please enter n and x '
READ(*,*)N,X
ans = bessel(N,X)
print *, ans
end program
SUBROUTINE SPHJ(N,X,NM,SJ,DJ)
! =======================================================
! Purpose: Compute spherical Bessel functions jn(x) and
! their derivatives
! Input : x --- Argument of jn(x)
! n --- Order of jn(x) ( n = 0,1,??? )
! Output: SJ(n) --- jn(x)
! DJ(n) --- jn'(x)
! NM --- Highest order computed
! Routines called:
! MSTA1 and MSTA2 for computing the starting
! point for backward recurrence
! =======================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION SJ(0:N),DJ(0:N)
NM=N
IF (DABS(X).EQ.1.0D-100) THEN
DO 10 K=0,N
SJ(K)=0.0D0
10 DJ(K)=0.0D0
SJ(0)=1.0D0
DJ(1)=.3333333333333333D0
RETURN
ENDIF
SJ(0)=DSIN(X)/X
SJ(1)=(SJ(0)-DCOS(X))/X
IF (N.GE.2) THEN
SA=SJ(0)
SB=SJ(1)
M=MSTA1(X,200)
IF (M.LT.N) THEN
NM=M
ELSE
M=MSTA2(X,N,15)
ENDIF
F0=0.0D0
F1=1.0D0-100
DO 15 K=M,0,-1
F=(2.0D0*K+3.0D0)*F1/X-F0
IF (K.LE.NM) SJ(K)=F
F0=F1
15 F1=F
IF (DABS(SA).GT.DABS(SB)) CS=SA/F
IF (DABS(SA).LE.DABS(SB)) CS=SB/F0
DO 20 K=0,NM
20 SJ(K)=CS*SJ(K)
ENDIF
DJ(0)=(DCOS(X)-DSIN(X)/X)/X
DO 25 K=1,NM
25 DJ(K)=SJ(K-1)-(K+1.0D0)*SJ(K)/X
RETURN
END
INTEGER FUNCTION MSTA1(X,MP)
! ===================================================
! Purpose: Determine the starting point for backward
! recurrence such that the magnitude of
! Jn(x) at that point is about 10^(-MP)
! Input : x --- Argument of Jn(x)
! MP --- Value of magnitude
! Output: MSTA1 --- Starting point
! ===================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
N0=INT(1.1*A0)+1
F0=ENVJ(N0,A0)-MP
N1=N0+5
F1=ENVJ(N1,A0)-MP
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-MP
IF(ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA1=NN
RETURN
END
INTEGER FUNCTION MSTA2(X,N,MP)
! ===================================================
! Purpose: Determine the starting point for backward
! recurrence such that all Jn(x) has MP
! significant digits
! Input : x --- Argument of Jn(x)
! n --- Order of Jn(x)
! MP --- Significant digit
! Output: MSTA2 --- Starting point
! ===================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
HMP=0.5D0*MP
EJN=ENVJ(N,A0)
IF (EJN.LE.HMP) THEN
OBJ=MP
N0=INT(1.1*A0)
ELSE
OBJ=HMP+EJN
N0=N
ENDIF
F0=ENVJ(N0,A0)-OBJ
N1=N0+5
F1=ENVJ(N1,A0)-OBJ
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-OBJ
IF (ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA2=NN+10
RETURN
END
REAL*8 FUNCTION ENVJ(N,X)
DOUBLE PRECISION X
ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
RETURN
END
!end of file msphj.f90
doubleprecision function bessel(N,X)
implicit doubleprecision(A-Z)
DIMENSION SJ(0:250),DJ(0:250)
integer :: N
CALL SPHJ(N,X,N,SJ,DJ)
bessel = SJ(N)
end function
这里是一个无法运行的示例程序,使用相同的功能。
program hello
IMPLICIT none
doubleprecision :: bessel, ans
integer :: N, X
WRITE(*,*)'Please enter n and x '
READ(*,*)N,X
ans = bessel(N,X)
print *, ans
end program
我对 Fortran 比较陌生,不明白为什么我的程序不能运行。我感谢任何人都可以提供的帮助。
我想非工作示例程序使用与工作示例相同的 bessel
实现。
如果是这样,则N
和X
的类型(非工作主程序中的integer
)与[=中的相应参数之间存在类型冲突11=] 根据规范都是双精度的
implicit doubleprecision(A-Z)
bessel 中的所有内容默认为 doubleprecision
。您的主程序必须将 N
和 X
定义为 doubleprecision
.
正如我在上面的评论中所说,最好的解决方案是在任何地方都使用显式类型。
我正在编写一个程序,该程序使用给定的子例程来计算球贝塞尔函数。我将给出 table 的子例程修改为只给出一个值的函数。然而,我意识到当我调用我的函数时我需要在我的程序中包含 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
否则它会给我一个错误的值或错误。下面我提供了一个可以正常工作的示例程序。
! n = 3, x = 2 should return ~ 6.07E-2
program hello
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
doubleprecision :: bessel, ans
WRITE(*,*)'Please enter n and x '
READ(*,*)N,X
ans = bessel(N,X)
print *, ans
end program
SUBROUTINE SPHJ(N,X,NM,SJ,DJ)
! =======================================================
! Purpose: Compute spherical Bessel functions jn(x) and
! their derivatives
! Input : x --- Argument of jn(x)
! n --- Order of jn(x) ( n = 0,1,??? )
! Output: SJ(n) --- jn(x)
! DJ(n) --- jn'(x)
! NM --- Highest order computed
! Routines called:
! MSTA1 and MSTA2 for computing the starting
! point for backward recurrence
! =======================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION SJ(0:N),DJ(0:N)
NM=N
IF (DABS(X).EQ.1.0D-100) THEN
DO 10 K=0,N
SJ(K)=0.0D0
10 DJ(K)=0.0D0
SJ(0)=1.0D0
DJ(1)=.3333333333333333D0
RETURN
ENDIF
SJ(0)=DSIN(X)/X
SJ(1)=(SJ(0)-DCOS(X))/X
IF (N.GE.2) THEN
SA=SJ(0)
SB=SJ(1)
M=MSTA1(X,200)
IF (M.LT.N) THEN
NM=M
ELSE
M=MSTA2(X,N,15)
ENDIF
F0=0.0D0
F1=1.0D0-100
DO 15 K=M,0,-1
F=(2.0D0*K+3.0D0)*F1/X-F0
IF (K.LE.NM) SJ(K)=F
F0=F1
15 F1=F
IF (DABS(SA).GT.DABS(SB)) CS=SA/F
IF (DABS(SA).LE.DABS(SB)) CS=SB/F0
DO 20 K=0,NM
20 SJ(K)=CS*SJ(K)
ENDIF
DJ(0)=(DCOS(X)-DSIN(X)/X)/X
DO 25 K=1,NM
25 DJ(K)=SJ(K-1)-(K+1.0D0)*SJ(K)/X
RETURN
END
INTEGER FUNCTION MSTA1(X,MP)
! ===================================================
! Purpose: Determine the starting point for backward
! recurrence such that the magnitude of
! Jn(x) at that point is about 10^(-MP)
! Input : x --- Argument of Jn(x)
! MP --- Value of magnitude
! Output: MSTA1 --- Starting point
! ===================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
N0=INT(1.1*A0)+1
F0=ENVJ(N0,A0)-MP
N1=N0+5
F1=ENVJ(N1,A0)-MP
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-MP
IF(ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA1=NN
RETURN
END
INTEGER FUNCTION MSTA2(X,N,MP)
! ===================================================
! Purpose: Determine the starting point for backward
! recurrence such that all Jn(x) has MP
! significant digits
! Input : x --- Argument of Jn(x)
! n --- Order of Jn(x)
! MP --- Significant digit
! Output: MSTA2 --- Starting point
! ===================================================
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
A0=DABS(X)
HMP=0.5D0*MP
EJN=ENVJ(N,A0)
IF (EJN.LE.HMP) THEN
OBJ=MP
N0=INT(1.1*A0)
ELSE
OBJ=HMP+EJN
N0=N
ENDIF
F0=ENVJ(N0,A0)-OBJ
N1=N0+5
F1=ENVJ(N1,A0)-OBJ
DO 10 IT=1,20
NN=N1-(N1-N0)/(1.0D0-F0/F1)
F=ENVJ(NN,A0)-OBJ
IF (ABS(NN-N1).LT.1) GO TO 20
N0=N1
F0=F1
N1=NN
10 F1=F
20 MSTA2=NN+10
RETURN
END
REAL*8 FUNCTION ENVJ(N,X)
DOUBLE PRECISION X
ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
RETURN
END
!end of file msphj.f90
doubleprecision function bessel(N,X)
implicit doubleprecision(A-Z)
DIMENSION SJ(0:250),DJ(0:250)
integer :: N
CALL SPHJ(N,X,N,SJ,DJ)
bessel = SJ(N)
end function
这里是一个无法运行的示例程序,使用相同的功能。
program hello
IMPLICIT none
doubleprecision :: bessel, ans
integer :: N, X
WRITE(*,*)'Please enter n and x '
READ(*,*)N,X
ans = bessel(N,X)
print *, ans
end program
我对 Fortran 比较陌生,不明白为什么我的程序不能运行。我感谢任何人都可以提供的帮助。
我想非工作示例程序使用与工作示例相同的 bessel
实现。
如果是这样,则N
和X
的类型(非工作主程序中的integer
)与[=中的相应参数之间存在类型冲突11=] 根据规范都是双精度的
implicit doubleprecision(A-Z)
bessel 中的所有内容默认为 doubleprecision
。您的主程序必须将 N
和 X
定义为 doubleprecision
.
正如我在上面的评论中所说,最好的解决方案是在任何地方都使用显式类型。