Fortran 代码给出错误

Fortran code giving error

我从一本书中获取了这个 Fortran 程序,该程序基本上对某些数据进行了拟合优度测试并给出了输出。代码及其实际 result/output 如下:

real*4 x(50),xc(50,20),omega(50)
integer ir(50)
real*8 xx
c This code tests goodness of fit.
n=47
c The method of Bak, Nielsen, and Madsen is used.
data (x(i), i=1,47)/ 18, 22, 26, 16, 19, 21, 18, 22,
* 25, 31, 30, 34, 31, 25, 21, 24, 21, 28, 24, 26, 32,
* 33, 36, 39, 32, 33, 42, 44, 43, 48, 50, 56, 57, 59,
* 51, 49, 49, 57, 69, 72, 75, 76, 78, 73, 73, 75, 86/
do 999 icase=1,2
c Parameter icase =1 or 2 denotes SDE model 1 or 2.
xx=102038.
m=8
h=1.0
do 10 j=1,m+1
10 omega(j)=0.0
kk=4
akk=kk
h=h/akk
do 202 i=2,n
xs=x(i-1)
xe=x(i)
do 202 j=1,m
xk=xs
do 252 k=1,kk
call functs(icase,xk,f,g)
call random(xx,rand1,rand2)
252 xk=xk+h*f+sqrt(h)*g*rand1
xc(i,j)=xk
202 continue
do 402 i=2,n
irr=1
do 302 j=1,m
xe=x(i)
xcalc=xc(i,j)
if(xe.gt.xcalc) irr=irr+1
302 continue
402 ir(i)=irr
do 502 i=2,n
irr=ir(i)
omega(irr)=omega(irr)+1.0
502 continue
chi2=0.0
an=n
am=m
hlp=(an-1.0)/(am+1.0)
do 602 j=1,m+1
602 chi2=chi2+(omega(j)-hlp)**2/hlp
write(6,100) icase,chi2
100 format(5x,i7,5x,f9.2)
999 continue
stop
end
subroutine functs(icase,x,f,g)
th1=3510.0
th2=13500.0
f=th1/(x*x)
g=th2/(x*x)
if(icase.eq.1) goto 17
th1=.0361
th2=.6090
f=th1*x
g=sqrt(th2*x)
17 continue
return
end
subroutine random(xx,rand1,rand2)
real*8 xx,a,b,d,rng(2)
a=16807.
ib=2147483647
b=ib
do 55 i=1,2
id=a*xx/b
d=id
xx=a*xx-d*b
55 rng(i)=xx/b
pi=3.141592654
u1=rng(1)
u2=rng(2)
hlp=sqrt(-2.0*alog(u1))
rand1=hlp*cos(pi*2.0*u2)
rand2=hlp*sin(pi*2.0*u2)
return
end

这个程序的输出是:

1 18.57
2 4.09

然而,即使在使用了许多在线 Fortran 编译器之后,我也没有得到这些结果。它给出了非标准类型声明等错误

我需要帮助才能获得与上述相同的输出。

该代码是使用(旧的)Fortran 77 风格编写的,并添加了一些常见的扩展。由于它使用所谓的固定形式,因此源代码使用的列对于拥有正确的代码至关重要。特别是案例:

  • 注释由第一列的 c 字符定义
  • 续行由第六列的 * 定义
  • 标签必须使用前 5 列
  • 常规代码必须使用 7-72 列范围

正确缩进您的代码允许它在 GNU gfortran(使用 v.4.8.2 测试)和 Intel ifort(使用版本 15.0.2 测试)上使用 运行。要通知编译器您想为大多数编译器采用固定格式,您只需为源文件使用 .f 扩展名。否则你有合适的编译器选项。对于 gfortran,编译指定 -ffixed-form。下面提供了(最少)缩进的代码。

      real*4 x(50),xc(50,20),omega(50)
      integer ir(50)
      real*8 xx
c This code tests goodness of fit.
      n=47
c The method of Bak, Nielsen, and Madsen is used.
      data (x(i), i=1,47)/ 18, 22, 26, 16, 19, 21, 18, 22,
     * 25, 31, 30, 34, 31, 25, 21, 24, 21, 28, 24, 26, 32,
     * 33, 36, 39, 32, 33, 42, 44, 43, 48, 50, 56, 57, 59,
     * 51, 49, 49, 57, 69, 72, 75, 76, 78, 73, 73, 75, 86/
      do 999 icase=1,2
c Parameter icase =1 or 2 denotes SDE model 1 or 2.
      xx=102038.
      m=8
      h=1.0
      do 10 j=1,m+1
10    omega(j)=0.0
      kk=4
      akk=kk
      h=h/akk
      do 202 i=2,n
      xs=x(i-1)
      xe=x(i)
      do 202 j=1,m
      xk=xs
      do 252 k=1,kk
      call functs(icase,xk,f,g)
      call random(xx,rand1,rand2)
252   xk=xk+h*f+sqrt(h)*g*rand1
      xc(i,j)=xk
202   continue
      do 402 i=2,n
      irr=1
      do 302 j=1,m
      xe=x(i)
      xcalc=xc(i,j)
      if(xe.gt.xcalc) irr=irr+1
302   continue
402   ir(i)=irr
      do 502 i=2,n
      irr=ir(i)
      omega(irr)=omega(irr)+1.0
502   continue
      chi2=0.0
      an=n
      am=m
      hlp=(an-1.0)/(am+1.0)
      do 602 j=1,m+1
602   chi2=chi2+(omega(j)-hlp)**2/hlp
      write(6,100) icase,chi2
100   format(5x,i7,5x,f9.2)
999   continue
      stop
      end
      subroutine functs(icase,x,f,g)
      th1=3510.0
      th2=13500.0
      f=th1/(x*x)
      g=th2/(x*x)
      if(icase.eq.1) goto 17
      th1=.0361
      th2=.6090
      f=th1*x
      g=sqrt(th2*x)
17    continue
      return
      end
      subroutine random(xx,rand1,rand2)
      real*8 xx,a,b,d,rng(2)
      a=16807.
      ib=2147483647
      b=ib
      do 55 i=1,2
      id=a*xx/b
      d=id
      xx=a*xx-d*b
55    rng(i)=xx/b
      pi=3.141592654
      u1=rng(1)
      u2=rng(2)
      hlp=sqrt(-2.0*alog(u1))
      rand1=hlp*cos(pi*2.0*u2)
      rand2=hlp*sin(pi*2.0*u2)
      return
      end

如果您想使用在线资源进行编译,请确保正确复制粘贴代码(使用正确的缩进)并使用固定形式的选项。例如在下面的 shell 中使用 https://www.tutorialspoint.com/compile_fortran_online.php 编译类型:gfortran -ffixed-form *.f95 -o main

由于 Fortran 77 风格现在已经很老了,如果您要开始编写新代码,我个人建议转向自由形式的源代码并使用更新的 Fortran 功能。下面给出了使用现代风格可能重写的代码:

module my_kinds
    integer, parameter :: sp = selected_real_kind(9)
    integer, parameter :: dp = selected_real_kind(18)
end module my_kinds

program test_from_book
    use my_kinds
    real(sp) :: x(50),xc(50,20),omega(50)
    integer  :: ir(50)
    real(dp) :: xx
    ! This code tests goodness of fit.
    n=47
    ! The method of Bak, Nielsen, and Madsen is used.
    x = [ 18, 22, 26, 16, 19, 21, 18, 22, &
          25, 31, 30, 34, 31, 25, 21, 24, 21, 28, 24, 26, 32, &
          33, 36, 39, 32, 33, 42, 44, 43, 48, 50, 56, 57, 59, &
          51, 49, 49, 57, 69, 72, 75, 76, 78, 73, 73, 75, 86, &
          0 ,  0,  0]
    loop_999: do icase=1,2
        ! Parameter icase =1 or 2 denotes SDE model 1 or 2.
        xx=102038.
        m=8
        h=1.0
        do j=1,m+1
            omega(j)=0.0
        enddo
        kk=4
        akk=kk
        h=h/akk
        loop_202: do i=2,n
            xs=x(i-1)
            xe=x(i)
            do j=1,m
                xk=xs
                do k=1,kk
                    call functs(icase,xk,f,g)
                    call random(xx,rand1,rand2)
                    xk=xk+h*f+sqrt(h)*g*rand1
                enddo
                xc(i,j)=xk
            enddo
        enddo loop_202
        loop_402: do i=2,n
            irr=1
            do j=1,m
                xe=x(i)
                xcalc=xc(i,j)
                if(xe.gt.xcalc) irr=irr+1
            enddo
            ir(i)=irr
        enddo loop_402
        do i=2,n
            irr=ir(i)
            omega(irr)=omega(irr)+1.0
        enddo
        chi2=0.0
        an=n
        am=m
        hlp=(an-1.0)/(am+1.0)
        do j=1,m+1
            chi2=chi2+(omega(j)-hlp)**2/hlp
        enddo
        write(6,100) icase,chi2
        100 format(5x,i7,5x,f9.2)
    enddo loop_999
    stop
end

subroutine functs(icase,x,f,g)
    th1=3510.0
    th2=13500.0
    f=th1/(x*x)
    g=th2/(x*x)
    if(icase.ne.1) then
        th1=.0361
        th2=.6090
        f=th1*x
        g=sqrt(th2*x)
    endif
end

subroutine random(xx,rand1,rand2)
    use my_kinds
    real(dp) :: xx,a,b,d,rng(2)
    a=16807.
    ib=2147483647
    b=ib
    do i=1,2
        id=a*xx/b
        d=id
        xx=a*xx-d*b
        rng(i)=xx/b
    enddo
    pi=3.141592654
    u1=rng(1)
    u2=rng(2)
    hlp=sqrt(-2.0*alog(u1))
    rand1=hlp*cos(pi*2.0*u2)
    rand2=hlp*sin(pi*2.0*u2)
end