并行嵌套循环中的数据竞争

Question

我有一个三重嵌套循环，我想并行化，但是，我遇到了数据争用问题。我很确定我需要以某种方式使用缩减，但我不太知道如何使用。

这是有问题的循环：

#pragma omp parallel for simd collapse(3)
    for (uint64 u = 0; u < nu; ++u) {
        for (uint64 e = 0; e < ne; ++e) {
            for (uint64 v = 0; v < nv; ++v) {
                uAT[u][e] += _uT[u][e][v] * wA[e][v];
            }
        }
    }

有人可以向我解释一下，为什么这会导致数据竞争？我真的很想了解这一点，这样我将来就不会遇到这些问题。另外，这个循环可以并行吗？如果是这样，怎么办？

编辑：我怎么知道存在数据竞争？

这个循环应该完成的任务（并且它是串行完成的）是计算不连续伽辽金框架中函数的元素平均值。当我多次运行代码时，有时会得到不同的结果，尽管它应该总是产生相同的结果。产生的错误值总是小于应有的值，这就是为什么我假设某些值没有被添加。也许这张图可以更好地解释它：第三个单元格中的平均值显然是错误的（太小）。

Answer 1

关于多线程而不是SIMD的原始答案

通过使用该collapse(3)子句，迭代的整个迭代空间nu * ne * nv被分布在线程之间。这意味着对于u和的任意组合e，v循环也可以分布在多个线程之间。然后，它们可以并行访问同一元素uAT[u][e]，这就是数据竞争。

只要 比nu * ne您使用的 CPU 核心数量大得多，最简单的解决方案就是改用collapse(2). collapse由于（请参阅此处）的实现可能效率低下，您甚至可能希望完全放弃它，具体取决于nu是否足够大。

如果您确实需要所有三个循环的并行性来有效地使用您的硬件，您可以使用其中之一reduction(+: uAT[0:nu][0:ne])（将其添加到现有的编译指示中）或#pragma omp atomic update在uAT[u][e] += ...;. 应该对这些选项中哪一个更快进行基准测试。该reduction子句将使用更多的内存，因为每个线程都获得通过 寻址的整个内存的自己的私有副本uAT。另一方面，atomic update在最坏的情况下，可能会顺序排列并行工作的一部分，并提供比使用更差的性能collapse(2)。

编辑1：SIMD

我刚刚看到您还关心 SIMD 而不仅仅是多线程。我原来的答案是关于后者的。对于 SIMD，我首先看一下没有任何 OpenMP 指令的编译器输出（例如在Compiler Explorer上）。如果您的编译器（具有目标处理器架构的优化和信息）已经使用 SIMD 指令，则您可能不需要首先使用 OpenMP。

如果您仍然想使用 OpenMP 来矢量化循环，则忽略该collapse子句并将编译指示放在第二个循环前面将是我的第一个 ansatz。将其放在第一个循环前面collapse(2)也可能有效。即使是reduction应该可以工作，但在这种情况下似乎不必要的复杂，因为我希望有足够的并行性nu或nu * ne填充您的 SIMD 通道。我从未在 SIMD 上下文中使用过如下所述的数组缩减，因此我不太确定它会做什么（即为每个 SIMD 通道分配一个数组似乎不太现实），或者它是否是 OpenMP 的一部分标准（取决于标准的版本，请参见此处）。

我认为，按照您现在编写代码的方式，多线程数据竞争的描述在技术上仍然适用。但我不确定您的代码是否会导致（有效）矢量化，因此编译后的二进制文件可能不会出现数据竞争（它可能根本不会矢量化）。

编辑2：线程+SIMD

nu = 4我对、nv = 128和ne之间1024的循环嵌套的几个版本进行了基准测试524288（2^19）。我的基准测试是使用 google-benchmark （即 C++ 而不是 C，我认为这里不重要）和 gcc 11.3 （始终使用-march=native -mtune=native，即对于可移植性能这可能没有帮助）完成的。我确保并行初始化所有数据（首次接触策略）以避免不良的 NUMA 影响。我稍微修改了问题/代码以使用连续内存并手动执行多维索引。由于OP的代码没有显示数据类型，我使用了float.

我将在这里分享的三个最好的版本都是三个性能相对相似的版本。因此，根据硬件架构的不同，哪一个最好可能会有所不同。对我来说，受 Peter Cordes 评论启发的版本在这个答案下面表现最好：

#pragma omp parallel for
for (uint64_t e = 0UL; e < ne; ++e) {
  // In the future we will get `#pragma omp unroll partial(4)`.
  // The unrolling might actually not be necessary (or even a pessimization).
  // So maybe leave it to the compiler.
  #pragma GCC unroll 4
  for (uint64_t u = 0UL; u < nu; ++u) {
    float temp = 0.0f;
    #pragma omp simd reduction(+ : temp)
    for (uint64_t v = 0UL; v < nv; ++v) {
      temp += uT[(u * ne + e) * nv + v] * wA[e * nv + v];
    }
    uAT[u * ne + e] += temp;
  }
}

还可以有一个float temp[nu];并将展开的u循环放入循环内v，以更接近 Peter Cordes 的描述，但随后必须使用如上所述的数组缩减。这些数组减少始终导致我出现堆栈溢出，因此我选择了这个版本，该版本取决于nv足够小，wA仍然可以在迭代之间缓存u。

第二个版本的不同之处在于u循环位于外部：

#pragma omp parallel
for (uint64_t u = 0UL; u < nu; ++u) {
  #pragma omp for
  for (uint64_t e = 0UL; e < ne; ++e) {
    float temp = 0.0f;
    #pragma omp simd reduction(+ : temp)
    for (uint64_t v = 0UL; v < nv; ++v) {
      temp += uT[(u * ne + e) * nv + v] * wA[e * nv + v];
    }
    uAT[u * ne + e] += temp;
  }
}

第三个版本使用collapse(2)：

#pragma omp parallel for collapse(2)
for (uint64_t u = 0UL; u < nu; ++u) {
  for (uint64_t e = 0UL; e < ne; ++e) {
    float temp = 0.0f;
    #pragma omp simd reduction(+ : temp)
    for (uint64_t v = 0UL; v < nv; ++v) {
      temp += uT[(u * ne + e) * nv + v] * wA[e * nv + v];
    }
    uAT[u * ne + e] += temp;
  }
}

长话短说：

我认为从三个版本中都可以看出最重要的一点：

• #pragma omp parallel for simd对于大的简单循环很有用。如果您有循环嵌套，您可能应该将编译指示分开。
• 用于simd在连续迭代中访问连续元素的循环。
• 通过尊重前两点，您不需要可能昂贵的数组缩减。
• 通过使用临时（即寄存器）进行减少而不是写回内存，您可以更轻松地使用 OpenMP，并且可能在串行代码中具有更好的性能。由于浮点非关联性，如果不通过例如（gcc）允许，大多数编译器都无法完成此优化-ffast-math。
• 出于同样的原因，大多数编译器无法自行对归约进行向量化。
• 像使用或循环的顺序这样的细节是collapse较小的优化，只要您提供足够的并行性，这些优化就不那么重要。即，不要并行化if很小（正如问题作者在这个答案下面所写的那样）。euunu