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

【问题标题】：Results of parallel program with nested loops differ from serial program具有嵌套循环的并行程序的结果与串行程序不同
【发布时间】：2016-07-20 15:40:57
【问题描述】：

我想为这个单线程代码使用 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（一定有问题）。

我做错了什么？

【问题讨论】：

对所有 Fortran 问题使用标签 ptag:fortran]。在需要区分的地方添加版本标签。注意订阅fortran 和fortran95 的人数，您希望更多人看到您的问题，不是吗？
几乎可以说并行程序总是产生与其串行表亲不同的结果。具体说明结果有何不同，从差异中可以学到很多东西。
请尝试使用标题来说明您的问题，而不仅仅是主题。主题由标签很好地定义。
@VladimirF：你什么意思？你的意思是我也想通过观看次数来获得人气？不，我从不打算让更多人了解我的问题，我真的不在乎。我只是想我会从一些专家那里得到最正确的答案。我什至不知道获得观看次数有什么好处，除非我可以有更大的机会从一些专家那里得到答案。你能告诉我，这个网站是否禁止像我一样标记更多？如果是这样，我很乐意将其删除。谢谢。
@HighPerformanceMark：不同之处在于：通过单个程序我得到 12.1212（它是真实的），通过并行我得到 0.000（一定有问题）。

标签： fortran openmp gfortran fortran95

【解决方案1】：

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

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

尽管如此，下面的代码在我的机器上产生了“正确”的答案。我不确定它有多灵活，但它可以在这里工作。

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）。

【讨论】：

是的，你说的没错。这不是一个真实的案例，只是一个想象的案例，并且肯定 [DO I=2, 30000] 绝对不是必需的。我只是添加它以了解并行代码的工作原理，并且不需要 [DO K=1, 50000] 以查看差异CPU时间的代表。是的，我忘了把 [HSTEP1X] 放在 [LASTPRIVATE] 上。顺便说一句，代码也适用于我的机器。你的解释很棒。
哦，是的，我忘了评论 [J] 可以选择定义为 [PRIVATE]，因为它在循环内部，但 [Q2] 和 [P2] 都应保留为 [PRIVATE ].

【解决方案2】：

您是否尝试过使用

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

【讨论】：

我认为它不起作用，因为这里的主要问题是 [HSTEP1X] 在循环 [DO J=Q2,P2] 之前有自己的值（零）。因此，如果我们使用 [DEFAULT(PRIVATE)]，则 [HSTEP1X] 的值将保持为零。
你可能是对的 bob.bob.bob，gfortran 似乎比 ifort 难多了。我今天安装了它，并在差异中挣扎。无论如何，让我们忽略 DEFAULT ......这里似乎需要 REDUCTION ？
您能具体说明一下“更难”是什么意思吗？由于根据我的经验，我使用 Gfortran 并且从未遇到任何严重问题，即使上周我使用新的 Intel fortran 16 时，我也收到消息，好像它不支持格式为 .F95 及更高版本的文件。最后我应该将格式更改为 .F90 （但使用旧的英特尔 Fortran 13，因为新的仍然无法工作）然后它工作了....是的，你完全可以忽略 [DEFAULT(SHARED)]在这种情况下，因为没关系。这只是我喜欢完整编写它的习惯。
嗯，我一直在寻找 gfortran ewuivilent to intel's -diag-file and -vector-report to get ，所以更难的意思是“我无法在英特尔的 1/2 小时内完成，而不是几年。” OpenMP 似乎依赖 -O# 进行任何矢量化...因此矢量消息，我在 !$OMP DO SIMD REDUCTION(+:K) 上遇到错误