積分法飛人的計算

Question

我已經有了現在的代碼，但它仍然不起作用。如果代碼是正確的，請幫助我如何編譯它。我曾試圖編譯它這樣：

gfortran trap.f -fopenmp

PROGRAM TRAP 
USE OMP_LIB 
DOUBLE PRECISION INTEG, TMPINT 
DOUBLE PRECISION A, B  
PARAMETER (A=3.0, B=7.0) 
INTEGER N  
PARAMETER (N=10) 
DOUBLE PRECISION H   
DOUBLE PRECISION X 
INTEGER I 
DOUBLE PRECISION F   
H = (B-A)/N 
INTEG = 0.0 
TMPINT = 0.0 
!$omp parallel firstprivate(X, TMPINT) shared(INTEG) 
!$omp do 
DO 10 I=1,N-1,1 
X=A+I*H 
TMPINT = TMPINT + F(X) 
10 CONTINUE 
!$omp end do 
!$omp critical 
INTEG = INTEG + TMPINT 
!$omp end critical 
!$omp end parallel 
NTEG = (INTEG+(F(A)+F(B))/2.0)*H 
PRINT *, "WITH N=", N, "INTEGRAL=", INTEG 
END 
FUNCTION F(X) 
DOUBLE PRECISION X 
F = X/(X + 1) * EXP(-X + 2) 
END 

編譯器提供了以下問題：

[http://i.stack.imgur.com/QPknv.png][1]

[http://i.stack.imgur.com/GYkmN.png][2]

Answer 1

你的程序有一個後綴爲.f，所以gfortran假設代碼是固定格式的，並且抱怨許多陳述是「不可分類的」。要解決此問題，請將文件名更改爲trap.f90並編譯爲gfortran -fopenmp trap.f90以採用自由格式。還有其他的問題：一個是功能F(X)返回類型不匹配在主程序聲明的類型，從而

FUNCTION F(X) 
implicit none      !<--- this is always recommended 
DOUBLE PRECISION X, F    !<--- add F here 
F = X/(X + 1) * EXP(-X + 2) 
END 

的另一個問題是NTEG可能是一個錯字修改F(X)需求的INTEG，所以它應修改爲

INTEG = (INTEG+(F(A)+F(B))/2.0)*H 

（如果我們有在主程序implicit none這是自動檢測的）。現在運行代碼，例如8個線程，給人

$ OMP_NUM_THREADS=8 ./a.out 
WITH N=   10 INTEGRAL= 0.28927708626319770 

而exact result爲0.28598 ...增加N值，我們可以確認該協議變得更好：

WITH N=   100 INTEGRAL= 0.28602065571967972 
WITH N=  1000 INTEGRAL= 0.28598803555916535 
WITH N=  10000 INTEGRAL= 0.28598770935198736 
WITH N=  100000 INTEGRAL= 0.28598770608991503 

---

BTW，它可能是更容易使用還原子句做同樣的事情，例如：

INTEG = 0.0 

!$omp parallel do reduction(+ : integ) private(x)         
DO I = 1, N-1                  
    X = A + I * H 
    INTEG = INTEG + F(X) 
ENDDO 
!$omp end parallel do                

INTEG = (INTEG+(F(A)+F(B))/2.0)*H 

Answer 2

您的代碼使用固定格式（.f）。因此，您必須按照固定格式的規則進行編碼：每行上的前六個字符具有特殊含義，除非您在第一個位置，一行續行（第六個位置）或語句標籤中指定註釋10 。

如果您相應地設置了代碼格式，編譯器會報告F(X)的返回值不匹配。由於您不使用implicit none，該類型由函數的第一個字母定義，並且F映射到（單精度）實數。所以你需要明確指定返回類型。

然後代碼如下：

 PROGRAM TRAP 
     USE OMP_LIB 
     DOUBLE PRECISION INTEG, TMPINT 
     DOUBLE PRECISION A, B  
     PARAMETER (A=3.0, B=7.0) 
     INTEGER N  
     PARAMETER (N=10) 
     DOUBLE PRECISION H   
     DOUBLE PRECISION X 
     INTEGER I 
     DOUBLE PRECISION F   
     H = (B-A)/N 
     INTEG = 0.0 
     TMPINT = 0.0 
c$omp parallel firstprivate(X, TMPINT) shared(INTEG) 
c$omp do 
     DO 10 I=1,N-1,1 
     X=A+I*H 
     TMPINT = TMPINT + F(X) 
10 CONTINUE 
c$omp end do 
c$omp critical 
     INTEG = INTEG + TMPINT 
c$omp end critical 
c$omp end parallel 
     INTEG = (INTEG+(F(A)+F(B))/2.0)*H 
     PRINT *, "WITH N=", N, "INTEGRAL=", INTEG 
     END 

     DOUBLE PRECISION FUNCTION F(X) 
     DOUBLE PRECISION X 
     F = X/(X + 1) * EXP(-X + 2) 
     END 

[請注意，我也固定NTAG =線進入INTEG=，因爲我相信這是預期。我沒有檢查代碼的有效性。 ]