2017-05-07 202 views
-1

我已經從一本基本上運行善良的書中讀取了這個Fortran程序,以適合測試某些數據並給出輸出結果。代碼和它的實際結果/輸出給定爲下:這個程序的Fortran代碼給出錯誤

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編譯器後,我沒有得到這些結果。它會給出非標準類型聲明等錯誤。

我需要幫助以獲得與上述相同的輸出。

+3

請給出確切的錯誤,它們出現的行,並顯示你不想嘗試解決它。 – Carcigenicate

+0

請使用一些縮進,空格和空行來使您的代碼可讀。 Seriusly,這太可怕了。並告訴我們你得到的確切的錯誤消息以及創建這些消息的編譯器命令。 –

回答

3

該代碼是使用(舊)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 
+0

謝謝,現代風格的代碼運作良好! –