OpenMP和C++：私有变量

Question

我对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：与当前的代码，执行结果不同于非并行执行。

Answer 1

需要重申的是，您的bodies变量不会随机出现;你应该确切地知道它在哪里被声明以及它被定义为什么。但是，因为您只访问bodies的元素并且永远不会更改它们，所以该变量无论如何都应该是shared，所以不是您的问题。

您的实际问题来自forces变量。您必须确保不同线程不会更改相同j的变量forces[j]。如果您遵循循环的逻辑，则可以确保forces[i]仅由不同的线程访问，因此它们之间不存在争用。但forces[j]对于相同的j可以非常容易地通过您的并行i循环的不同迭代进行修改。您需要按照StackOverflow链接的答案之一执行reduce on your array。

Answer 2

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; 
} 

它采用n*(n-1)迭代和易于并行。

但由于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]，并且迭代不再被很好地分配，即第一个线程比最后一个线程运行更多的迭代。

的简单的解决方案是parellelize内环

​​

这使用n*nthreads关键操作在总共n*(n-1)/2迭代。因此，随着n变大，关键操作的成本会变得更小。您可以为每个线程使用私有的forces向量，并将它们合并到关键部分，但我不认为这是必需的，因为关键操作位于外部循环而不是内部循环。

---

这里是一个解决方案，它fuses the triangular loop允许每个线程运行在相同数目的迭代。

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函数），所以我想出了基于this answer一个更简单的解决方法。

这将三角形映射到矩形，反之亦然。首先，我将其转换为一个宽度为n但带有n*(n-1)/2（与三角形相同）的矩形。然后I I利用以下公式

//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]; 
}