具有嵌套循环的并行程序的结果不同于串行程序

Question

我想将 OpenMP 用于此单线程代码：

PROGRAM SINGLE
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2

!Just intiial value
  DO I=1, 30000
     SUMGRM(I)=I*3        
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000    
  END DO

!Computing computer's running time (start) : for serial programming
 CALL CPU_TIME(TIME1)

 DO K=1, 50000
    DO I=2, 30000
       HSTEP1X=0.0    
         DO J=SUMGRM(I-1)+1, SUMGRM(I)-1
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X       
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X         
    END DO
 END DO

  PRINT *, 'Results  =', HSTEP1X
  PRINT *, '   '

!Computing computer's running time (finish) : for serial programming
 CALL CPU_TIME(TIME2)
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM SINGLE

如您所见，主要问题位于最内侧循环（J），这也是最外侧循环（I）的一个函数。我试过像这样并行化这个程序：

PROGRAM PARALLEL
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2,OMP_GET_WTIME
  INTEGER::Q2,P2

!Just intiial value
  DO I=1, 30000
     SUMGRM(I)=I*3        
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000  
  END DO

!Computing computer's running time (start) : for parallel programming
 TIME1= OMP_GET_WTIME()

 DO K=1, 50000
 !$OMP PARALLEL DO PRIVATE (HSTEP1X,Q2,P2)
    DO I=2, 30000
       HSTEP1X=0.0
       Q2=SUMGRM(I-1)+1
       P2=SUMGRM(I)-1
         DO J=Q2, P2
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X       
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X     
    END DO
 !$OMP END PARALLEL DO
 END DO

 PRINT *, 'Results  =', HSTEP1X
 PRINT *, '   '

!Computing computer's running time (finish) : for parallel programming
 TIME2= OMP_GET_WTIME()
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PARALLEL

我用的是gfortran with -O3 -fopenmp然后导出OMP_NUM_THREADS=...并行程序运行速度更快但是结果和单线程代码不一样。通过串行程序我得到了 12.1212 （这是正确的）并且通过并行我得到了 0.000 （一定有问题）。

我做错了什么？

Answer 1

首先我们可以注意到，默认情况下您可能会发现 j 和 hstep2x 都将在线程之间共享。我不认为这真的是你想要的，因为如果多个线程使用相同的迭代索引但试图在不同的范围内循环，它会导致一些非常奇怪的行为。

接下来让我们注意，您的串行代码实际上只是打印 i=30000 迭代的结果，因为 hstep1x 的值在每次迭代开始时重置为 0。因此，为了在 openmp 代码中获得 "correct" 答案，我们可以只专注于重现最终迭代——我认为这完全否定了在这里使用 openmp 的意义。我猜这只是一个简单的案例，您试图用它来代表您的实际问题——我认为您在制作这个时可能错过了一些实际问题。

然而，下面的代码在我的机器上产生了 "correct" 答案。我不确定它有多灵活，但它在这里有效。

PROGRAM PARALLEL
  INTEGER, DIMENSION(30000)::SUMGRM
  INTEGER, DIMENSION(90000)::GRI,H
  REAL*8::HSTEP1X,HSTEP2X
  REAL*8::TIME1,TIME2,OMP_GET_WTIME
  INTEGER::Q2,P2

!Just intiial value                                                                                                                                                                                                  
  DO I=1, 30000
     SUMGRM(I)=I*3
  END DO

  DO I=1, 90000
     GRI(I)=I
     H(I)=0.5*I/10000
  END DO

!Computing computer's running time (start) : for parallel programming                                                                                                                                                
 TIME1= OMP_GET_WTIME()

 DO K=1, 50000
!$OMP PARALLEL DO PRIVATE (Q2,P2,J,HSTEP2X) DEFAULT(SHARED) LASTPRIVATE(HSTEP1X)                                                                                                                                     
    DO I=2, 30000
       HSTEP1X=0.0
       Q2= SUMGRM(I-1)+1
       P2= SUMGRM(I)-1
         DO J=Q2,P2
            HSTEP2X=H(GRI(J))/0.99
            HSTEP1X=HSTEP1X+HSTEP2X
         END DO
       HSTEP2X=H(GRI(SUMGRM(I)))/0.99
       HSTEP1X=HSTEP1X+HSTEP2X
    END DO
!$OMP END PARALLEL DO                                                                                                                                                                                                
END DO

 PRINT *, 'Results  =', HSTEP1X
 PRINT *, '   '

!Computing computer's running time (finish) : for parallel programming                                                                                                                                               
 TIME2= OMP_GET_WTIME()
 PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PARALLEL

我在这里做了三件事：

确保 j 和 hstep2x 对每个线程都是私有的。
明确声明要共享的默认行为（此处不需要，但没关系）。
指定 hstep1x 为 lastprivate。这意味着在退出并行区域后，hstep1x 的值取自执行最后一次迭代的线程。（详见 here）。

Answer 2

您尝试过使用

!$OMP PARALLEL DO DEFAULT(PRIVATE) REDUCTION(+:HSTEP1X)

具有嵌套循环的并行程序的结果不同于串行程序

Results of parallel program with nested loops differ from serial program

fortran

openmp

gfortran

fortran95