fortran - OpenMP 中对变量求和给出错误答案的线程

Question

为了练习并行化 do 循环，我在 Fortran 中进行以下积分

$\integral{0}{1} \frac{4}{1+x^{2}} = \pi$

以下是我实现的代码:

      program mpintegrate

      integer i,nmax,nthreads,OMP_GET_NUM_THREADS
      real xn,dx,value
      real X(100000)

      nthreads = 4
      nmax = 100000
      xn = 0.0
      dx = 1.0/nmax
      value = 0.0

      do i=1,nmax
         X(i) = xn
         xn = xn + dx
      enddo

      call OMP_SET_NUM_THREADS(nthreads)

!$OMP Parallel 

!$OMP Do Schedule(Static) Private(i,X)

      do i=1,nmax
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo

!$OMP End DO NoWait

!$OMP End Parallel

      print *, value

      end

我编译程序没有问题

gfortran -fopenmp -o mpintegrate mpintegrate.f

问题出在我执行程序时。当我按原样运行程序时，我得到的值范围为 (1,4)。但是，当我使用 omp do 循环取消注释 print 语句时，最终值大约是它应该的值 pi。

为什么value中的答案不正确？

Answer 1

这里的一个问题是 X 不需要是私有(private)的(并且需要在平行线上指定，而不是在 do 线上)；每个人都需要看到它，并且为每个线程拥有单独的副本是没有意义的。更糟糕的是，您从此处访问私有(private)副本获得的结果是未定义的，因为一旦您进入私有(private)区域，该私有(private)变量就尚未初始化。您可以使用 firstprivate 而不是 private，后者使用并行区域之前的内容为您初始化它，但这里最简单/最好的只是 shared.

让结束不等待也没有多大意义，因为结束并行无论如何都必须等待每个人都完成。

但是，话虽这么说，您仍然存在一个相当重大(且经典)的正确性问题。如果您在循环中更明确一点，那么这里发生的事情就会更清楚(为了清楚起见，放弃时间表，因为问题不依赖于所选的时间表):

!$OMP Parallel do Private(i) Default(none) Shared(value,X,dx,nmax)

      do i=1,nmax
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo

!$OMP End Parallel Do
print *, value

重复运行会给出不同的值:

$ ./foo
   1.6643878
$ ./foo
   1.5004054
$ ./foo
   1.2746993

问题是所有线程都写入同一个共享变量值。这是错误的 - 每个人都同时写入，结果是乱码，因为一个线程可以计算它自己的贡献，准备将其添加到 value 中，就在它即将这样做时，另一个线程可以这样做它写入值，然后很快就会被破坏。并发写入同一个共享变量是一个经典race condition ，一个标准的错误系列，在 OpenMP 等共享内存编程中特别常见。

除了错误之外，它还很慢。由于内存系统中的争用，多个线程争夺相同的几个字节的内存(内存足够接近以落入同一高速缓存行)可能会非常慢。即使它们不是完全相同的变量(如本例所示)，这种内存争用 - False Sharing如果它们恰好是相邻变量，则可以显着减慢速度。取出显式线程号设置，并使用环境变量:

$ export OMP_NUM_THREADS=1
$ time ./foo
   3.1407621

real    0m0.003s
user    0m0.001s
sys 0m0.001s

$ export OMP_NUM_THREADS=2
$ time ./foo
   3.1224852

real    0m0.007s
user    0m0.012s
sys 0m0.000s

$ export OMP_NUM_THREADS=8
$ time ./foo
   1.1651508

real    0m0.008s
user    0m0.042s
sys 0m0.000s

因此，使用更多线程运行时，速度会慢近 3 倍(而且越来越错误)。

那么我们可以做些什么来解决这个问题呢？我们可以做的一件事是使用 atomic 指令确保每个人的添加不会互相覆盖:

!$OMP Parallel do Schedule(Static) Private(i) Default(none) Shared(X,dx, value, nmax)
      do i=1,nmax
!$OMP atomic
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo
!$OMP end parallel do

这解决了正确性问题:

$ export OMP_NUM_THREADS=8
$ ./foo
   3.1407621

但对速度问题没有任何作用:

$ export OMP_NUM_THREADS=1
$ time ./foo
   3.1407621

real    0m0.004s
user    0m0.001s
sys 0m0.002s

$ export OMP_NUM_THREADS=2
$ time ./foo
   3.1407738

real    0m0.014s
user    0m0.023s
sys 0m0.001s

(请注意，使用不同数量的线程，您会得到略有不同的答案。这是由于最终总和的计算顺序与串行情况不同。对于单精度实数，由于不同的原因，差异会出现在第 7 位数字中。操作顺序很难避免，这里我们进行了 100,000 次操作。)

那么我们还能做什么呢？一种方法是让每个人记录自己的部分总和，然后在完成后将它们全部加起来:

!... 

      integer, parameter :: nthreads = 4
      integer, parameter :: space=8
      integer :: threadno
      real, dimension(nthreads*space) :: partials
!...

      partials=0

!...

!$OMP Parallel Private(value,i,threadno) Default(none) Shared(X,dx, partials)
      value = 0
      threadno = omp_get_thread_num()
!$OMP DO
      do i=1,nmax
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo
!$OMP END DO
      partials((threadno+1)*space) = value
!$OMP end parallel

      value = sum(partials)
      print *, value
      end

这是有效的 - 我们得到了正确的答案，如果你研究线程的数量，你会发现它非常快 - 我们已经在部分和数组中间隔开条目以避免错误共享(并且它是false，这一次，因为每个人都写入数组中的不同条目 - 不会覆盖)。

尽管如此，仅仅为了在线程之间获得正确的总和，这是一项愚蠢的工作量!有一种更简单的方法可以做到这一点 - OpenMP 有一个 reduction构造自动执行此操作(并且比上面的手工版本更有效:)

!$OMP Parallel do reduction(+:value) Private(i) Default(none) Shared(X,dx)
      do i=1,nmax
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo
!$OMP end parallel do
      print *, value 

现在程序可以正常运行，速度很快，而且代码也相当简单。最终的代码，在更现代的 Fortran 语言中，看起来像这样:

      program mpintegrate
      use omp_lib
      integer, parameter :: nmax = 100000
      real :: xn,dx,value
      real :: X(nmax)
      integer :: i

      integer, parameter :: nthreads = 4

      xn = 0.0
      dx = 1.0/nmax
      value = 0.0
      partials=0

      do i=1,nmax
         X(i) = xn
         xn = xn + dx
      enddo

      call omp_set_num_threads(nthreads)

!$OMP Parallel do reduction(+:value) Private(i) Default(none) Shared(X,dx)
      do i=1,nmax
         value = value + dx*(4.0/(1+X(i)*X(i)))
      enddo
!$OMP end parallel do
      print *, value

      end