计算科学 - OpenMP：对相互作用力的并行计算 - 吾爱随笔录

OpenMP：对相互作用力的并行计算

计算科学计算物理学并行计算 C

2021-12-13 04:44:42

我最近开始学习 OpenMP。尽管我已经有了一些直觉，但我仍然对如何在某些对计算物理学家非常有用的情况下进行操作有一些疑问。

我的问题：在分子模拟中寻找中心力是典型的，其中对于粒子和。此外，我们可以考虑牛顿第三定律，然后对于一个有个相互作用粒子的系统，我们只需要运行一半的粒子，因为另一个只是简单地通过改变符号获得。在 C/C++中： $F=F(r)$ $r = |\vec{r}_i -\vec{r}_j|$ $i$ $j$ $N$ $^{(1)}$

 for(int i=0;i<N;i++){
    for(int j=i+1;j<N;j++) {
            F = force(r);
            f[i] += F;
            f[j] -= F;
        }
     }

这是一个简单的例子，需要对给定的向量进行求和和减法，即使不涉及中心力，类似的例子也很常见。但是，我不确定如何有效地并行化上述循环。一种选择可能是计算所有内容，而不是只计算一半：

#pragma omp for private(F), reduction(+:f[:N])
  for(int i=0;i<N;i++){
    for(int j=0;j<N;j++) {
            if(i!=j){
                F = force(x);
                f[i] += F;
                f[j] -= F;
            }
        }
     }

但这远非理想，因为 F 的计算可能非常昂贵，并且会进行两次。我看到并行化第一个代码的主要问题是减少，因为它必须用于减法reduction(-:f[:N])和加法，reduction(+:f[:N])以及随之而来的冲突。

我的替代方法是通过使用两个中间数组，一个用于求和组件，另一个用于减法，让我们调用它们sums_f和subs_f：

double *sums_f, *subs_f;
// Reserve memory for sums_f and subs_f and set them to zero
#pragma omp parallel
{
    #pragma omp for reduction(+: sums_f[:N], -:subs_f[:N]) 
    for(int i=0;i<N;i++){
        for(int j=i+1;j<N;j++) {
                F = force(r);
                sums_f[i] += F;
                subs_f[j] -= F;
            }
        }
    }
 
    #pragma omp for
    for(int i=0;i<N;i++)
        f[i] += sums_f[i] - subs_f[i];
 }

但我认为对于大的内存而言，这可能会非常昂贵。 $N$

总而言之，我看不出有一种替代方法可以使它们处理效率高，而不需要几个变量或额外的步骤来降低并行化的效率。您对此类案件的处理有什么建议吗？

注：[1] 更详细的解释参见：How to compute force in multi-particle MD

更新：

20210513.1：重要提示：实际上动态分配的数组不适用于reduction. 这可能排除了它在这种情况下的使用。

20210514.1：对于 3D 对象，我使用将每个粒子的力数据存储为行乘 3 列的 3D 维数组。但是，这意味着必须执行 3x2步骤。结果是失去了并行化的时间增益。更好的解决方案是将力数据存储为一维数组。 $N$ atomic

2个回答

不一定是答案，只是在这里更容易格式化代码。该#pragma omp atomic指令确保在并行循环中避免竞争条件。它类似于critical但仅适用于简单的表达式，例如+=，并且速度更快。我相信这是实现您的想法的最简单方法。

#pragma omp parallel for
    for (int i = 0; i < N; i++)
    {
        for (int j = i + 1; j < N; j++)
        {
            F = force(r);
#pragma omp atomic
            f[i] += F;
#pragma omp atomic
            f[j] -= F;
        }
    }

编辑：现在更聪明了

如果您想force(r)更均匀地分布工作，请定义

k = j (j + 1) / 2 - i - 1

$k = j(j+1)/2 - i - 1$

j = ⌈ \sqrt{(2 k + 2.25)} - 0.5 ⌉

$j = \lceil\sqrt{(2k+2.25)}-0.5\rceil$

i = j (j + 1) / 2 - k - 1

$i = j(j+1)/2 - k - 1$

不要让我证明它创造了一个独特的配对——这是在它起作用之前搞砸和尝试的结果。我已经对其进行了测试N=45000（超出您的需要long）。这应该避免必须显式创建(i,j)对。

std::pair<int, int> ij(int k)
{
    int j = std::ceil(sqrt(2 * k + 2.25) - 0.5);
    int tri_j = j * (j + 1) / 2;
    int i = tri_j - k - 1;
    return std::make_pair(i, j);
}

//number of unique pairs (i,j)
int M = N*(N-1)/2;

#pragma omp parallel for
for (int k = 0; k < M; k++)
{
    std::pair<int,int> ij = calc_ij(k);
    int i = ij.first;
    int j = ij.second;

    //presumably the bottleneck
    double F = force(r);
    
//practically instantaneous
#pragma omp atomic
    f[i] += F;
#pragma omp atomic
    f[j] -= F;
}

我很想看看线性索引是否会提高执行时间。

所以我决定实现这个，当然我发现了一堆我在之前的答案中没有预见到的因素，我删除了这些因素，因为它包含一些专利废话。

顺序码

      for (int ip=0; ip<N; ip++) {
        for (int jp=ip+1; jp<N; jp++) {
          struct force f = force_calc(points[ip],points[jp]);
          add_force( forces+ip,f );
          sub_force( forces+jp,f );
        }
      }

add/sub 例程要考虑力箭头的不同方向。

通过循环完整的 N^2 并行并忽略额外的工作

    for (int step=0; step<STEPS; step++) {
      for (int ip=0; ip<N; ip++) {
        double sumx=0., sumy=0., sumf=0.;
#pragma omp parallel for reduction(+:sumx,sumy,sumf)
        for (int jp=0; jp<N; jp++) {
          if (ip==jp) continue;
          struct force f = force_calc(points[ip],points[jp]);
          sumx += f.x; sumy += f.y; sumf += f.f;
        } // end parallel jp loop
        struct force sumforce = {sumx,sumy,sumf};
        add_force( forces+ip, sumforce );
      } // end ip loop

请注意，我们只能并行化内部循环。如果您尝试并行化两者并添加折叠子句，您将得到不正确的结果。

并行代码 2，在三角形上循环，但在sub_force例程中使用原子：

#pragma omp parallel for schedule(guided,4)
      for (int ip=0; ip<N; ip++) {
        for (int jp=ip+1; jp<N; jp++) {
          struct force f = force_calc(points[ip],points[jp]);
          add_force( forces+ip,f );
          sub_force( forces+jp,f );
        }
      }

现在是时间安排。我正在使用总共 56 个内核的双插槽 Cascade Lake。

================ #threads = 1 ================
               Sequential: 1.167970e+01; total force: 7.382722e+08
       Full loop Parallel: 3.132034e+00; total force: 7.382722e+08, speedup= 3.73
   Atomic update parallel: 1.136083e+01; total force: 7.382722e+08, speedup= 1.03
================ #threads = 18 ================
               Sequential: 1.166757e+01; total force: 7.589678e+08
       Full loop Parallel: 1.910710e+00; total force: 7.589678e+08, speedup= 6.11
   Atomic update parallel: 1.733956e+00; total force: 7.588650e+08, speedup= 6.73
================ #threads = 37 ================
               Sequential: 1.167301e+01; total force: 7.378324e+08
       Full loop Parallel: 4.362817e+00; total force: 7.378324e+08, speedup= 2.68
   Atomic update parallel: 1.087576e+00; total force: 7.376554e+08, speedup=10.73
================ #threads = 56 ================
               Sequential: 1.168075e+01; total force: 7.428077e+08
       Full loop Parallel: 9.035676e+00; total force: 7.428077e+08, speedup= 1.29
   Atomic update parallel: 9.829427e-01; total force: 7.425870e+08, speedup=11.88

第一个大惊喜：在没有并行性的情况下，双循环快了 4 倍（不是因为双倍工作而慢了 2 倍）。这是因为三角形版本会在整个内存中跳跃，这意味着更糟糕的缓存和 TLB 行为。

惊喜二：双循环越来越差。我猜测阵列太短（10k 个粒子）而无法获得良好的效率。

惊喜三：使用原子更新实际上是高效的。我们得到了不错的加速，尽管远非线性。ip我正在使用指导计划，因为每次迭代的工作都不是恒定的。也许如果力计算更昂贵，这将更好地扩展。

编辑，但为了获得最佳解决方案，我们使所有内容都原子化，并折叠两个循环。现在还有很多工作要让所有核心都满意。

#pragma omp parallel for collapse(2)
      for (int ip=0; ip<N; ip++) {
        for (int jp=0; jp<N; jp++) {
          if (ip==jp) continue;
          struct force f = force_calc(points[ip],points[jp]);
          add_force( forces+ip, f );
        } // end parallel jp loop
      } // end ip loop

其中add例程使用#omp atomic更新。

定时：

================ #threads = 1 ================
               Sequential: 2.029093e+01;
         Full loop atomic: 3.114875e+01; speedup= 0.65
================ #threads = 18 ================
         Full loop atomic: 1.739044e+00; speedup=11.67
================ #threads = 37 ================
         Full loop atomic: 8.486451e-01; speedup=23.95
================ #threads = 56 ================
         Full loop atomic: 5.630970e-01; speedup=36.09

这让我有点意外。现代硬件必须真正知道如何处理这些原子的东西。

其它你可能感兴趣的问题

上一篇拍摄方法的首选应用下一篇用 numpy 形成一个特定的（平均）块矩阵