隐式 none 使程序不准确

Question

我正在编写一个程序，该程序使用给定的子例程来计算球形贝塞尔函数。我将给表的子例程修改为只给出一个值的函数。但是，我意识到当我调用我的函数时，我需要在我的程序中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相对较新，不明白为什么我的程序不起作用。我感谢任何人都可以提供的任何帮助。

Answer 1

我想非工作示例程序使用与工作示例相同的bessel实现。

如果是这样，则N和X的类型（integer在非工作主程序中）与bessel中的相应参数之间存在类型冲突，这些参数根据规范都是双精度

implicit doubleprecision(A-Z)

默认情况下，贝塞尔中的所有内容都是doubleprecision.您的主程序必须将N和X定义为doubleprecision。

正如我在上面的评论中所说，最好的解决方案是在任何地方使用显式类型。