使用OpenMP并行化嵌套循环

Question

我有一个相对简单的循环,我正在使用蛮力方法计算粒子系统的净加速度.

我有一个工作的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

Answer 1

使用您当前的方法和数据结构,您将很难通过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

...,现在回到问题]

这里有两个问题:

1. 正如你已经知道的那样,你有一个数据竞争更新bodyArray(j)%acc; 多个线程将尝试更新相同的元素,并且没有这些更新的协调.垃圾结果.使用critical部分或命令语句序列化代码; 当你做对了,你也得到它与开始使用OpenMP之前一样慢.
2. 访问元素的模式bodyArray是缓存不友好的.我发现即使您在没有序列化计算的情况下处理数据竞争也不会感到惊讶,缓存不友好的影响是产生比原始代码更慢的代码.现代CPU的计算速度非常快,但内存系统很难为野兽提供动力,因此缓存效果可能非常大.试图同时在同一个rank-1数组上运行两个循环,这实际上是你的代码所做的,永远不会(？)以最大速度通过缓存移位数据.

我个人会尝试以下.我不打算保证这会更快,但是(我认为)比修复当前的方法更容易,并且像手套一样适合OpenMP. 我确实怀疑这是一个过于复杂的问题,但我还没有更好的想法.

首先,创建一个实数的二维数组,调用它forces,其中元素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将更多,等等.你可能只是简单地schedule(dynamic)在parallel指令中添加一个子句,你可能需要更加努力地平衡负载.

您还可以查看我的代码片段WRT缓存友好和适当的调整.你可能会发现,如果你按照我的建议做,你的原始代码会更好,那么减半计算量并不会节省太多时间.

另一种方法是将较低(或较高)三角形打包forces成秩-1阵列,并使用一些花哨的索引算法将2D (i,j)索引转换为该阵列的1D索引.这样可以节省存储空间,并且可以更容易实现缓存友好.