OpenMP和C++:私有变量

OpenMP and C++: private variables

本文关键字：变量 C++ OpenMP      更新时间：2023-10-16

我对OpenMP和c++很陌生，也许正因为如此，我遇到了一些非常基本的问题。

我试图在所有变量都是私有的情况下进行静态调度（以防万一，以验证获得的结果与非并行的结果相同）。

当我看到像bodies这样的变量时，问题就出现了，我不知道它们来自哪里，因为它们以前没有定义过。

是否可以将所有出现的变量（如bodies）定义为私有变量？怎么能做到

  std::vector<phys_vector> forces(bodies.size());
  size_t i, j; double dist, f, alpha;

  #pragma omp parallel for schedule(static) private(i, j, dist, f, alpha)
  for (i=0; i<bodies.size(); ++i) {
    for (j = i+1; j<bodies.size(); ++j) {
      dist = distance(bodies[i], bodies[j]);
      if (dist > param.min_distance()) {
        f = attraction(bodies[i], bodies[j], param.gravity(), dist);
        alpha = angle(bodies[i],bodies[j]);
        phys_vector deltaf{ f * cos(alpha) , f * sin(alpha) };
        forces[i] += deltaf;
        forces[j] -= deltaf;
      }
    }
  }
  return forces;
}

PS：对于当前代码，执行结果与非并行执行不同。

NoseKnowsAll已正确识别您的问题。

我想进一步解释一下为什么会出现这个问题。你可以用这样的方环来完成：

#pragma omp parallel for
for(int i=0; i<n; i++) {
    if(i==j) continue;
    phys_vector sum = 0;
    for(int j=0; j<n; j++) {
        //calculate deltaf
        sum += deltaf;
    }
    forces[i] = sum;
}

其使用CCD_ 3迭代并且易于并行化。

但由于force(i,j) = -force(j,i)，我们可以在一半的迭代中做到这一点，n*(n-1)/2，使用三角形循环（这就是您所做的）：

for(int i=0; i<n; i++) {
    phys_vector sum = 0;
    for(int j=i+1; j<n; j++) {
        //calculate deltaf
        sum += deltaf;
        forces[j] -= deltaf;
    }
    forces[i] = sum;
}

问题是，当你进行这种优化时，会使外循环的并行化变得更加困难。有两个问题：写入forces[j]，迭代不再分布良好，即第一个线程比最后一个线程运行更多的迭代。

简单的解决方案是将内环并行化

#pragma omp parallel
for(int i=0; i<n; i++) {
    phys_vector sum = 0;
    #pragma omp for
    for(int j=i+1; j<n; j++) {
        //calculate deltaf
        sum += deltaf;
        forces[j] -= deltaf;
    }
    #pragma omp critical
    forces[i] += sum;
}

这使用了总共n*(n-1)/2次迭代中的n*nthreads次关键操作。因此，关键操作的成本随着n的增大而减小。你可以为每个线程使用一个专用的forces向量，并将它们合并到一个关键部分，但我认为这没有必要，因为关键操作在外循环上，而不是在内循环上。

这里有一个解决方案，它融合了三角形循环，允许每个线程运行相同数量的迭代。

unsigned n = bodies.size();
unsigned r = n*(n-1)/2;
#pragma omp parallel
{
    std::vector<phys_vector> forces_local(bodies.size());
    #pragma omp for schedule(static)
    for(unsigned k=0; k<r; k++) {
        unsigned i  = (1 + sqrt(1.0+8.0*k))/2;
        unsigned j = i - k*(k-1)/2;
        //calculate deltaf
        forces_local[i] += deltaf;
        forces_local[j] -= deltaf;
    }
    #pragma omp critical
    for(unsigned i=0; i<n; i++) forces[i] += forcs_local[i];
}

我对以前融合三角形的方法不满意（因为它需要使用浮点和sqrt函数），所以我基于这个答案想出了一个简单得多的解决方案。

这将三角形映射为矩形，反之亦然。首先，我将其转换为宽度为n但宽度为n*(n-1)/2的矩形（与三角形相同）。然后我计算矩形的（行、列）值，然后使用以下公式映射到三角形（跳过对角线）

//i is the row, j is the column of the rectangle
if(j<=i) {
    i = n - i - 2;
    j = n - j - 1;
}

让我们选择一个例子。考虑以下n=5三角环路对

(0,1), (0,2), (0,3), (0,4)
       (1,2), (1,3), (1,4)
              (2,3), (2,4)
                     (3,4)

将其映射到矩形成为

(3,4), (0,1), (0,2), (0,3), (0,4)
(2,4), (2,3), (1,2), (1,3), (1,4)

具有偶数值的三角形循环以相同的方式工作，尽管可能不那么明显。例如n = 4。

(0,1), (0,2), (0,3)
       (1,2), (1,3)
              (2,3)

这就变成了

(2,3), (0,1), (0,2), (0,3)
(1,2), (1,3)

这并不完全是一个矩形，但映射的工作原理是一样的。我本可以把它映射成

 (0,1), (0,2), (0,3)
 (2,3), (1,2), (1,3)

它是一个矩形，但我需要两个奇数和偶数三角形大小的公式。

以下是使用矩形到三角形映射的新代码。

unsigned n = bodies.size();
#pragma omp parallel
{
    std::vector<phys_vector> forces_local(bodies.size());
    #pragma omp for schedule(static)
    for(unsigned k=0; k<n*(n-1)/2; k++) {
        unsigned i = k/n;
        unsigned j = k%n;
        if(j<=i) {
            i = n - i - 2;
            j = n - j - 1;
        }
        //calculate deltaf
        forces_local[i] += deltaf;
        forces_local[j] -= deltaf;
    }
    #pragma omp critical
    for(unsigned i=0; i<n; i++) forces[i] += forcs_local[i];
}

需要重申的是，您的bodies变量并不是凭空随机出现的；你应该弄清楚它在哪里声明，定义为什么。但是，因为你只访问bodies的元素，而从不更改它们，所以这个变量无论如何都应该是shared，所以这不是你的问题。

您的实际问题来自forces变量。您必须确保不同线程不会更改同一j的变量forces[j]。如果遵循循环的逻辑，就可以确保forces[i]只由不同的线程访问，因此它们之间没有争用。但是相同j的forces[j]可以很容易地通过并行i循环的不同迭代来修改。您需要做的是按照StackOverflow链接中的一个答案来减少您的数组。

相关文章：