使用OpenMP的critical和ordered指令

4

我对Fortran和OpenMP都比较新,但我正在努力适应。我有一段用于计算变异图的代码,现在想要并行化它。然而,我似乎遇到了竞争条件,因为一些结果会出现千分之一左右的偏差。

问题似乎是关于约简操作。使用OpenMP约简可以正常工作并给出正确的结果,但这并不是理想的解决方案,因为这些约简实际上发生在另一个子程序中(我将相关行复制到了OpenMP循环中进行测试)。因此,我将约减操作放入CRITICAL部分但没有成功。有趣的是,问题只会在实数情况下出现,而整数则不受影响。我已经考虑过加法的顺序是否有影响,但它们不应该产生这么大的误差。

为了确认,我将并行do中的所有内容都放在ORDERED块中,这当然会给出正确的结果(尽管没有任何加速)。我还尝试将所有内容放在CRITICAL部分内,但由于某种原因,这并没有给出正确的结果。我的理解是,OpenMP将在进入/退出CRITICAL部分时刷新共享变量,因此不应该存在任何缓存问题。

所以我的问题是:为什么在这种情况下CRITICAL部分无法正常工作?

以下是我的代码。除了nptmhmgam之外,所有共享变量都是只读的。

编辑:我试图通过将do循环替换为相同范围内的随机整数(即生成一对i,j在循环中;如果它们被“访问”,则产生新的随机数)来模拟多个线程引起的随机性,但出乎意料的是结果匹配了。然而,进一步检查发现我忘记给RNG提供种子,结果只是巧合而已。多尴尬啊!

TL;DR: 结果的不一致性是由浮点值的排序引起的。使用双精度可以解决问题。

!$OMP PARALLEL DEFAULT(none) SHARED(nd, x, y, z, nzlag, nylag, nxlag, &
!$OMP& dzlag, dylag, dxlag, nvarg, ivhead, ivtail, ivtype, vr, tmin, tmax, np, tm, hm, gam) num_threads(512)
!$OMP DO PRIVATE(i,j,zdis,ydis,xdis,izl,iyl,ixl,indx,vrh,vrt,vrhpr,vrtpr,variogram_type) !reduction(+:np, tm, hm, gam)
  DO i=1,nd        
!$OMP CRITICAL (main)
! Second loop over the data:
    DO j=1,nd

! The lag:
      zdis = z(j) - z(i)
      IF(zdis >= 0.0) THEN
        izl =  INT( zdis/dzlag+0.5)
      ELSE
        izl = -INT(-zdis/dzlag+0.5)
      END IF
 ! ---- SNIP ----

! Loop over all variograms for this lag:

      DO cur_variogram=1,nvarg
        variogram_type = ivtype(cur_variogram)

! Get the head and tail values:

        indx = i+(ivhead(cur_variogram)-1)*maxdim
        vrh   = vr(indx)
        indx = j+(ivtail(cur_variogram)-1)*maxdim
        vrt   = vr(indx)
        IF(vrh < tmin.OR.vrh >= tmax.OR. vrt < tmin.OR.vrt >= tmax) CYCLE

        ! ----- PROBLEM AREA -------
        np(ixl,iyl,izl,1)  = np(ixl,iyl,izl,1) + 1.   ! <-- This never fails
        tm(ixl,iyl,izl,1)  = tm(ixl,iyl,izl,1) + vrt  
        hm(ixl,iyl,izl,1)  = hm(ixl,iyl,izl,1) + vrh
        gam(ixl,iyl,izl,1) = gam(ixl,iyl,izl,1) + ((vrh-vrt)*(vrh-vrt))
        ! ----- END OF PROBLEM AREA -----

        !CALL updtvarg(ixl,iyl,izl,cur_variogram,variogram_type,vrt,vrh,vrtpr,vrhpr)
      END DO
    END DO
    !$OMP END CRITICAL (main)
  END DO
!$OMP END DO
!$OMP END PARALLEL

非常感谢您提前的帮助!

2
我建议尝试发布一个SSCCE,而不是一个不完整的代码片段。这将迫使您更好地分析实际引起问题的原因,并让其他人有可能在自己的环境中重现该问题。 - Massimiliano
@Massimiliano 通常我会同意你的看法,但是在这种情况下由于初始化和数据的问题,无法提供一个自包含的实例。 - Painless
@HighPerformanceMark 我报告的错误是绝对的。相对误差非常小(例如得到84.26539与84.26538之类的结果)。然而,绝对误差仍然似乎太大了,我无法忽略它。我还通过随机遍历do循环的数字范围来测试操作的顺序,并且产生了精确的结果。我也怀疑这背后是否有任何SIMD魔法。 - Painless
1
使用更大精度的实数类型,看看它会做什么。 - Vladimir F Героям слава
@VladimirF 是的,这显著提高了结果,所以感谢您的建议! - Painless
显示剩余3条评论
2个回答

3
如果你使用32位浮点数和算术,那么84.26539与84.26538之间的差异,即最低有效数字为1的差异,完全可以通过并行浮点运算的不确定性来解释。请记住,32位浮点数只有约7个十进制数字可供使用。
普通浮点运算不是严格联想的。对于实数(在数学而不是Fortran意义上)有(a+b)+c==a+(b+c),但是对于浮点数没有这样的规则。这在维基百科上关于浮点运算的文章中得到了很好的解释。
非确定性是由于在使用OpenMP时,操作的顺序被交给了运行时控制。跨线程值的求和(例如对+的缩减)将全局求和表达式的括号放置留给了运行时。甚至不能保证对同一个OpenMP程序的两次执行会产生完全相同的结果。
我怀疑即使在单个线程上运行OpenMP程序也可能会产生与等效的非OpenMP程序不同的结果。由于对OpenMP可执行文件可用线程数量的知识可能被推迟到运行时,因此编译器将需要创建并行化可执行文件,无论最终是否并行运行。

那是我的最初想法,但正如我在之前的评论中所说,我随机地遍历了两个do循环(即:我随机生成了两个1到nd之间的数字i和j;如果我之前“访问”过它们,则重新生成;继续直到所有nd * nd组合都被访问)。这产生了完全相同的结果。此外,OpenMP reduce语句(在发布的代码片段中已注释掉)每次也应该产生不同的结果,但实际上并没有。 - Painless
我不同意“OpenMP reduce语句...每次都应该产生不同的结果”,但我同意“OpenMP reduce语句可能不会每次产生相同的结果”。检查您的代码后,我没有得出其他结论或建议。 - High Performance Mark
有趣。我之前没有考虑到OpenMP中关联浮点顺序会如何影响结果。 - Z boson
嗯,似乎整个东西都被排序搞糊涂了,至少在使用OpenMP时是这样,但仍然无法解释我的实验结果,其中我随机迭代了索引。 - Painless
没关系,我是个傻瓜;我没有意识到init_seed()会自动生成一个种子,所以默认种子的结果是正确的。现在感觉有点尴尬 B-) - Painless

2

High Performance Mark提出了一个有趣的观点,与浮点数和结合律有关。这可以在C语言中轻松测试。

float a = -1.0E8f, b = 1.0E8f, c = 1.23456f;
printf("sum %f\n", a+b+c);   //output 1.234560
printf("sum %f\n", a+(b+c)); //output 0.000000

但我想指出,在OpenMP中有可能保持顺序。 我在这里讨论了这个问题:C++ OpenMP: Split for loop in even chunks static and join data at the end 编辑: 实际上,我混淆了可交换性和结合律。 如果您有一个是结合但不是可交换的运算符,则可以像我在上面的帖子中所做的那样使用OpenMP保留顺序。 但是,IEEE浮点数是可交换的但不是结合的,因此唯一能做的就是打破IEEE并使其成为可结合的。

这是一个相当不错的想法,如果可以的话,我一定会使用它! - Painless
@Painless,我犯了一个错误。我把结合律和交换律搞混了。我在那个链接中使用的技巧是针对一个不满足交换律但仍然满足结合律的运算符的。这种方法并不适用于反过来的情况(IEEE所要求的)。我已经更新了我的答案。 - Z boson

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接