积分法飞人的计算

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=，因为我相信这是预期。我没有检查代码的有效性。 ]