为什么OpenMP并行for循环不能正常工作？

Question

我想实现OpenMP来并行化我的代码。我从一个非常基本的例子开始，了解它是如何工作的，但我遗漏了一些东西......

所以，我的例子看起来像这样，没有并行化：

int main() {
  ...

  for (i = 0; i < n-1; i++) {
    u[i+1] = (1+h)*u[i]; // Euler
    v[i+1] = v[i]/(1-h); // implicit Euler
  }

  ...

  return 0;
}

我省略了“......”中的某些部分，因为它们不相关。它工作，如果我在文件上打印u[]和v[]数组，我得到预期的结果。

现在，如果我尝试通过添加以下内容来并行化它：

#include <omp.h>

int main() {
  ...

  omp_set_num_threads(2);

  #pragma omp parallel for
  for (i = 0; i < n-1; i++) {
    u[i+1] = (1+h)*u[i]; // Euler
    v[i+1] = v[i]/(1-h); // implicit Euler
  }

  ...

  return 0;
}

代码编译并运行程序，但u[]和v[]数组是半满的零。

如果我设置omp_set_num_threads( 4 )，我得到四分之三的零。 如果我设置omp_set_num_threads( 1 )，我会得到预期的结果。

所以它看起来只有第一个线程正在被执行，而不是其他的...

我究竟做错了什么？

Answer 1

OpenMP假定循环的每次迭代都独立于其他循环。当你写这个：

for (i = 0; i < n-1; i++) {
  u[i+1] = (1+h)*u[i]; // Euler
  v[i+1] = v[i]/(1-h); // implicit Euler
}

循环的迭代i正在修改迭代i+1。同时，迭代i+1可能同时发生。

除非您可以使迭代独立，否则这不是并行性的好用例。

而且，如果你考虑Euler的方法做了什么，显然不可能以这种方式并行处理你正在处理的代码。欧拉方法基于时间t+1处的信息计算时间t处的系统状态。因为你不知道在t+1知道什么是t，所以没有办法在Euler方法的迭代中并行化。

Answer 2

在并行化代码之前，必须确定它的并发性，即同时逻辑上发生的一组任务，然后找出一种方法使它们实际并行发生。

如上所述，由于其性质不存在并发性，因此这不是应用并行性的好例子。由于所谓的竞争条件，试图使用这样的并行性将导致错误的结果。

如果您只是想了解OpenMP的工作原理，请尝试提供一些示例，您可以清楚地识别出独立于概念的任务。我能想到的最简单的一种方法是通过积分计算曲线下面积。

Answer 3

欢迎来到并行（或“正常” - 并发）多个计算现实。

---

为什么？

处理循环的任何非顺序计划都会出现隐藏（未正确处理）数据泄露的问题 - { - access | -value}及时完整。

一个纯粹的[SERIAL]处理流程没有这样的危险，因为主要序列化的步骤间接引入（通过一个严格的顺序执行除了一步一个接一个作为序列）的顺序，其中没有机会同时“触摸”相同的内存位置两次或更多次。

一旦一个过程进入"just"-[CONCURRENT]或true-[PARALLEL]处理，这种“安心”就会无意中丢失。

突然间，有一个几乎随机的顺序（在“just”-[CONCURRENT]的情况下）或主要的“立即”奇点（避免任何“order”的原始含义 - 在真正的[PARALLEL]代码执行模式的情况下 - 像具有6DoF的机器人以真实[PARALLEL]的方式到达每个轨迹点，以纯粹的[SERIAL]方式并行驱动所有6个DoF轴，而不是一个接一个地驱动，而不是因为机器人手臂的3D轨迹将变得难以预测并且相互碰撞经常发生，所以现在有些 - 现在 - 其他 - 其他 - 后来 - 其余 - 它正以“正义”-[CONCURRENT]方式获得汽车装配线......）。

---

解：

使用称为原子操作的防御工具或主要方法 - 在可能的情况下设计（b）无锁定算法，或明确地发出信号并协调读写操作（当然，需要花费超时和降低的性能），为了保证价值不会被损坏成不一致的数字垃圾，如果保护步骤（确保所有“旧”写作在任何“下一步” - 前进之前安全“通过”以获得“正确” - 值）没有编码（如上所述）。

---

结语：

使用像OpenMP这样的工具来解决问题，它无法带来任何优势，这将导致花费时间和性能降低（因为需要处理所有与工具相关的开销，而在这种情况下，并行性的净效果几乎为零，算法不允许任何并行性的享受），所以最后一个人最终得到的方式更多。

了解OpenMP最佳实践的一个好点可能来自劳伦斯利弗莫尔国家实验室（确实非常称职）和类似的publications on using OpenMP.。

Answer 4

u[i+1] = (1+h)*u[i];                                                                                                                                              
v[i+1] = v[i]/(1-h);                                                                                                                                       

相当于

u[i] = pow((1+h), i)*u[0];
v[i] = v[0]*pow(1.0/(1-h), i);

因此，您可以像这样并行化代码

#pragma omp parallel for
for (int i = 0; i < n; i++) {
    u[i] = pow((1+h), i)*u[0];
    v[i] = v[0]*pow(1.0/(1-h), i);
}

如果你想减轻pow函数的成本，你可以每个线程执行一次，而不是像他的每次迭代一样（自t << n）。

#pragma omp parallel
{
    int nt = omp_get_num_threads();
    int t  = omp_get_thread_num();
    int s  = (t+0)*n/nt;
    int f  = (t+1)*n/nt;
    u[s]   = pow((1+h), s)*u[0];
    v[s]   = v[0]*pow(1.0/(1-h), s);
    for(int i=s; i<f-1; i++) {
        u[i+1] = (1+h)*u[i];
        v[i+1] = v[i]/(1-h);
    }
}

您还可以编写自己的针对整数幂优化的pow(double, int)函数。

请注意，我使用的关系实际上并非100％等效，因为浮点运算不是关联的。这通常不是问题，但这是人们应该注意的事情。