使用OpenMP并行化嵌套循环

Question

使用OpenMP并行化嵌套循环

5

我有一个相对简单的循环，使用暴力方法计算粒子系统的净加速度。

我拥有一个工作正常的OpenMP循环，它循环遍历每个粒子，并将其与其他每个粒子进行比较，具有n^2的复杂度：

   !$omp parallel do private(i) shared(bodyArray, n) default(none)
   do i = 1, n
        !acc is real(real64), dimension(3)
        bodyArray(i)%acc = bodyArray(i)%calcNetAcc(i, bodyArray)
   end do

这个程序已经可以很好地工作。

现在我想通过仅计算每个物体上的力一次来减少计算时间，利用了F(a->b) = -F(b->a)的事实，将需要计算的相互作用数量减半（n^2 / 2）。我在以下循环中实现了此操作：

call clearAcceleration(bodyArray) !zero out acceleration

    !$omp parallel do private(i, j) shared(bodyArray, n) default(none)
       do i = 1, n
          do j = i, n

             if ( i /= j .and. j > i) then 
                bodyArray(i)%acc = bodyArray(i)%acc + bodyArray(i)%accTo(bodyArray(j))
                bodyArray(j)%acc = bodyArray(j)%acc - bodyArray(i)%acc
             end if

          end do
       end do

但我在并行化这个循环时遇到了很多困难，一直得到垃圾结果。我认为问题出在这一行代码：

bodyArray(j)%acc = bodyArray(j)%acc - bodyArray(i)%acc

并且各种 'j' 的写入没有正确地加起来。我尝试使用原子语句，但是数组变量不允许这样做。然后我尝试了关键部分，但这增加了大约20倍的时间，仍然没有给出正确的结果。我还尝试添加有序语句，但是结果都是 NaN。有没有简单的方法可以使用OpenMP使此循环正常工作？

以下是代码，它略微提高了速度，但不是我想要的 ~2 倍速度。

!$omp parallel do private(i, j) shared(bodyArray, forces, n) default(none) schedule(guided)
      do i = 1, n
         do j = 1, i-1
               forces(j, i)%vec = bodyArray(i)%accTo(bodyArray(j))
               forces(i, j)%vec = -forces(j, i)%vec
         end do
      end do

      !$omp parallel do private(i, j) shared(bodyArray, n, forces) schedule(static)
      do i = 1, n
         do j = 1, n
            bodyArray(i)%acc = bodyArray(i)%acc + forces(j, i)%vec
         end do
      end do

- Exascale

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- High Performance Mark · Accepted Answer

根据您目前的方法和数据结构，使用OpenMP将很难获得良好的加速效果。建议考虑循环嵌套。

!$omp parallel do private(i, j) shared(bodyArray, n) default(none)
   do i = 1, n
      do j = i, n

         if ( i /= j .and. j > i) then 
            bodyArray(i)%acc = bodyArray(i)%acc + bodyArray(i)%accTo(bodyArray(j))
            bodyArray(j)%acc = bodyArray(j)%acc - bodyArray(i)%acc
         end if

      end do
   end do

[实际上，在您考虑之前，按照以下方式进行修订...]

!$omp parallel do private(i, j) shared(bodyArray, n) default(none)
   do i = 1, n
      do j = i+1, n

            bodyArray(i)%acc = bodyArray(i)%acc + bodyArray(i)%accTo(bodyArray(j))
            bodyArray(j)%acc = bodyArray(j)%acc - bodyArray(i)%acc

      end do
   end do

...,现在回到问题上

这里有两个问题：

正如你已经意识到的那样，你正在更新 bodyArray(j)%acc，存在数据竞争的问题；多个线程将尝试更新相同的元素，并且这些更新没有协调。结果不可靠。使用 critical sections 或排序语句可以使代码串行化；当你正确地使用 OpenMP 时，你也会发现它与未使用 OpenMP 时一样慢。
对 bodyArray 元素的访问模式是不友好的缓存。如果您没有序列化计算而解决了数据竞争，我不会感到惊讶，缓存不友好性的影响是产生比原始代码更慢的代码。现代 CPU 在计算方面非常快，但内存系统很难满足其需求，因此缓存效应可能是巨大的。同时运行两个循环来同时遍历相同的 rank-1 数组，在本质上就是你的代码所执行的操作，永远 (?) 不会以最大速度通过缓存移动数据。

个人建议尝试以下方法。我不能保证这会更快，但比修复当前方法并完美适配 OpenMP 更容易（我认为）。我有一个困扰我已久的疑惑，即这可能会使事情过于复杂，但是我还没有更好的想法。

首先，创建一个名为 forces 的实数 2D 数组，其中元素 force(i,j) 是元素 i 对 j 施加的力。然后，使用类似以下的代码（未经测试，如果您想追踪此行，请自行负责）：

forces = 0.0  ! Parallelise this if you want to
!$omp parallel do private(i, j) shared(forces, n) default(none)
   do i = 1, n
      do j = 1, i-1
            forces(i,j) = bodyArray(i)%accTo(bodyArray(j)) ! if I understand correctly
      end do
   end do

然后对每个粒子的力进行求和（以下是正确的内容，但我还没有仔细检查过）。

!$omp parallel do private(i) shared(bodyArray,forces, n) default(none)
do i = 1, n
      bodyArray(i)%acc = sum(forces(:,i))
end do

如我之前所写，计算速度非常快，如果你有足够的内存，通常值得用一些空间换取时间。

现在你可能遇到了一个关于forces循环嵌套中负载平衡的问题。大多数OpenMP实现默认情况下会执行静态工作分配（这不是标准要求，但似乎最常见，请检查你的文档）。因此，线程1将获得要处理的前n/num_threads行，但这些是你正在计算的三角形顶部微小的行。线程2将获得更多的工作，线程3仍然会获得更多的工作，以此类推。你可以尝试在parallel指令中添加一个schedule(dynamic)子句来解决这个问题，或者你可能需要更加努力地平衡负载。

你还可以查看我的代码片段，关于缓存友好性进行调整。如果你按照我的建议去做，你可能会发现原来的代码更好，减少计算量并不能节省太多时间。

另一种方法是将forces的下（或上）三角形打包成一个秩为1的数组，并使用一些高级索引算法将二维(i,j)索引转换为该数组中的一维索引。这将节省存储空间，并且可能更容易使缓存友好。