Fortran 中 OpenMP 并行 do 循环的优化

Question

背景

我正在使用 Fortran90 和 OpenMP 模拟分子动力学中 N 个带电粒子的运动。施加到每个离子 i 的力的解析表达式是已知的，并且是离子 i 和其他离子 (r_x、r_y、r_z) 位置的函数。我使用并行的 2 嵌套 do 循环计算每对离子之间的库仑相互作用。我可以确定循环结束时每个离子的加速度（a2_x、a2_y、a2_z）（然后用 velocity-Verlet 更新速度和位置）。

方法

我在我的程序中使用以下代码来计算应用于每个离子的库仑力。我计算下一个时间步的加速度 (a2_x)，从当前时间步的位置 (r_x) 开始。这是一个 3D 问题，我写了所有的行，但其中大部分对于 x、y 和 z 都是一样的，所以在第一次阅读时你可以只考虑 _x 变量来了解它是如何工作的。

我在 C 线程上并行化我的循环，ia 和 ib 是用于将 N 离子分成 C 部分的数组。例如C=4线程和N=16离子（Se编辑下面的评论）

integer, parameter :: ia(C) = [1,5,9,13]
integer, parameter :: ib(C) = [4,8,12,16]

然后库仑计算如下

!$omp parallel default(none) &
!$omp private(im, i,j,rji,r2inv) &
!$omp firstprivate(r_x,r_y,r_z, N, ia, ib) &
!$omp shared(a2_x, a2_y, a2_z)
    im = omp_get_thread_num() + 1 ! How much threads
  ! Coulomb forces between each ion pair
  ! Compute the Coulomb force applied to ion i
    do i = ia(im,1), ib(im,1) ! loop over threads
        do j = 1, N ! loop over all ions
            rji(1)  = r_x(j) - r_x(i) ! distance between the ion i and j over x
            rji(2)  = r_y(j) - r_y(i) ! over y
            rji(3)  = r_z(j) - r_z(i) ! over z
          ! then compute the inverse square root of distance between the current ion i and the neighbor j
            r2inv   = 1.d0/dsqrt(rji(1)*rji(1) + rji(2)*rji(2) + rji(3)*rji(3) + softening)
            r2inv   = r2inv * r2inv * r2inv * alpha(1) ! alpha is 1/4.pi.eps0
          ! computation of the accelerations
            a2_x(i) = a2_x(i) - rji(1)*r2inv
            a2_y(i) = a2_y(i) - rji(2)*r2inv
            a2_z(i) = a2_z(i) - rji(3)*r2inv
        enddo
    enddo
    !$omp end parallel

问题

我正在尝试优化程序中这个耗时的部分。操作数量相当多，扩展速度很快，有N个。你能告诉我你对这个程序的看法吗？我有一些具体问题。

1. 我被告知我应该将位置r_x、r_y和r_z作为private变量，这对我来说似乎违反直觉，因为我想使用先前定义的离子位置进入该循环，所以我使用firstprivate。那正确吗 ？

2. 我不确定关于其他变量的并行化是否是最佳的。 rji 和 r2inv 不应该共享吗？因为为了计算离子 i 和 j 之间的距离，我“超越”了线程，你明白我的意思了吗？我需要分布在两个不同线程上的离子之间的信息。

3. 我首先分裂离子的方式是否最佳？

4. 我分别为每个离子遍历所有离子，这将在计算离子 i 和 i 之间的距离时导致除以零。为了防止这种情况，我定义了一个软化变量，它的值非常小，因此它不完全为零。我这样做是为了避免 if i==i 耗时。

5. 平方根也可能很耗时？

   如需任何其他详细信息，请随时询问。

   编辑（备注）

   1. 我的电脑有一个 10 核 CPU Xeon W2155、32 Go RAM。我打算渲染大约 1000 个离子，同时考虑 4000 个，这需要很多时间。

   2. 我在其他可能消耗一些 CPU 时间的子例程中有这个 Coulomb 子例程。例如，一个可能很耗时的例程专门为每个离子生成随机数，具体取决于它们是否已经被激发，并应用正确的效果，无论它们是否吸收光子。所以这是很多 RNG，如果每个离子。

      编辑（命题测试）

      3. 将 !$omp do 与 schedule(dynamic,1)、schedule(guided) 或 schedule(nonmonotonic:dynamic) 和/或 collapse(2) 结合使用并没有缩短运行时间。它至少延长了三倍。建议我的模拟中的元素数量 (N) 太低，看不到明显的改进。如果我尝试渲染更多数量的元素（4096、8192 ...），我会尝试这些选项。

      4. 使用 !$omp do 而不是在内核之间自制离子分布确实显示了运行时间方面的等效性。实施起来更容易，我会保留它。

      5. 用**(-1/2) 替换逆dsqrt 表明在运行时间方面是等效的。

      6. 延迟平方根并将其与 r2inv 的三次方相结合也是等效的。所以我把整个系列的操作换成了**(-1.5)。

      7. 与rji(1)*r2inv 的想法相同，我之前执行rji*r2inv 并且仅在下一行中使用结果。

Answer 1

1. 一般来说，并行区只需要读取的变量可以是shared。但是，在某些情况下，每个线程都有 firstprivate 副本可以提供更好的性能（副本可以在每个内核的本地缓存中），特别是对于重复读取的变量。
2. 绝对不是！如果你这样做，这些变量就会出现竞争条件
3. 看起来不错，但使用 !$OMP DO 指令而不是手动将工作分配到不同的线程通常更简单（并且在最坏的情况下也是有效的）

       !$OMP DO
       do i = 1, N ! loop over all ions
           do j = 1, N ! loop over all ions
   

   4. 为什么不呢，前提是你可以选择一个不会改变你的模拟的softening值（这是你必须针对if解决方案进行测试的东西）
   5. 不知何故，但在某些时候你无法避免求幂。我会延迟 sqrt 和这样的划分：

   r2inv   = (rji(1)*rji(1) + rji(2)*rji(2) + rji(3)*rji(3) + softening)
   r2inv   = r2inv**(-1.5) * alpha(1) ! alpha is 1/4.pi.eps0
   

   将工作除以 2

   这些力是对称的，对于给定的 (i,j) 对只能计算一次。这也很自然地避免了 i==j 的情况和软化值。但是，迭代之间的工作负载是高度不平衡的，并且需要 dynamic 子句。这实际上是一个案例，手动将迭代分发到线程可以更有效；）...

   !$omp parallel default(none) &
   !$omp private(im, i,j,rji,r2inv) &
   !$omp firstprivate(r_x,r_y,r_z, N, ia, ib) &
   !$omp shared(a2_x, a2_y, a2_z)
     ! Coulomb forces between each ion pair
     ! Compute the Coulomb force applied to ion i
       !$omp do schedule(dynamic,1)
       do i = 1, N-1 ! loop over all ions
           do j = i+1, N ! loop over some ions
               rji(1)  = r_x(j) - r_x(i) ! distance between the ion i and j over x
               rji(2)  = r_y(j) - r_y(i) ! over y
               rji(3)  = r_z(j) - r_z(i) ! over z
             ! then compute the inverse square root of distance between the current ion i and the neighbor j
               r2inv   = (rji(1)*rji(1) + rji(2)*rji(2) + rji(3)*rji(3))
               r2inv   = r2inv**(-1.5) * alpha(1) ! alpha is 1/4.pi.eps0
             ! computation of the accelerations
               rji(:) = rji(:)*r2inv
               a2_x(i) = a2_x(i) - rji(1)
               a2_y(i) = a2_y(i) - rji(2)
               a2_z(i) = a2_z(i) - rji(3)
               a2_x(j) = a2_x(j) + rji(1)
               a2_y(j) = a2_y(j) + rji(2)
               a2_z(j) = a2_z(j) + rji(3)
           enddo
       enddo
       !$omp end do
       !$omp end parallel
   

   或者，可以使用 guided 子句，在迭代中进行一些更改以使第一个迭代中的工作负载较低：

       !$omp do schedule(guided)
       do i = 2, N ! loop over all ions
           do j = 1, i-1 ! loop over some ions