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