蒙特卡羅pi逼近的並行化

Question

我正在寫一個c腳本來平行pi近似與OpenMp。我認爲我的代碼工作正常，有一個令人信服的輸出。我現在用4個線程運行它。我不確定的是，如果這段代碼容易受到競爭狀況的影響？如果是這樣，我如何協調此代碼中的線程操作？

的代碼如下：

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <math.h> 
#include <omp.h> 

double sample_interval(double a, double b) { 

    double x = ((double) rand())/((double) RAND_MAX); 
    return (b-a)*x + a; 

} 

int main (int argc, char **argv) { 


    int N = atoi(argv[1]); // convert command-line input to N = number of points 
    int i; 
    int NumThreads = 4; 
    const double pi = 3.141592653589793; 
    double x, y, z; 
    double counter = 0; 



    #pragma omp parallel firstprivate(x, y, z, i) reduction(+:counter) num_threads(NumThreads) 
{ 
    srand(time(NULL)); 
    for (int i=0; i < N; ++i) 
    { 
    x = sample_interval(-1.,1.); 
    y = sample_interval(-1.,1.); 
    z = ((x*x)+(y*y)); 

    if (z<= 1) 
    { 
     counter++; 
    } 
    } 

} 
    double approx_pi = 4.0 * counter/ (double)N; 

    printf("%i %1.6e %1.6e\n ", N, 4.0 * counter/ (double)N, fabs(4.0 * counter/ (double)N - pi)/pi); 


    return 0; 

} 

此外，我在想，如果隨機數種子應內部或外部並行申報。我的輸出如下所示：

10 3.600000e+00 1.459156e-01 
100 3.160000e+00 5.859240e-03 
1000 3.108000e+00 1.069287e-02 
10000 3.142400e+00 2.569863e-04 
100000 3.144120e+00 8.044793e-04 
1000000 3.142628e+00 3.295610e-04 
10000000 3.141379e+00 6.794439e-05 
100000000 3.141467e+00 3.994585e-05 
1000000000 3.141686e+00 2.971945e-05 

現在看起來不錯。你對種族條件和種子安置的建議是最受歡迎的。

Answer 1

您可以看到代碼中存在一些問題。主要的是從我的觀點來看，它並不是並行的。或者更確切地說，在編譯OpenMP時，您並未啓用並行機制。這裏是一個可以看到的方式：

代碼並行的方式，主要for循環應全部由所有線程執行（這裏沒有工作共享，不#pragma omp parallel for，只有#pragma omp parallel）。因此，考慮到將線程數設置爲4，全局迭代次數應爲4*N。因此，你的輸出應該慢慢收斂到4 * Pi，而不是Pi。

確實，我在我的筆記本電腦上試過了你的代碼，使用OpenMP支持編譯它，而這正是我得到的。但是，當我不啓用OpenMP時，我得到的輸出與您的類似。因此，最後，您需要：

1. 在編譯時啓用OpenMP以獲取並行版本的代碼。
2. 通過NumThreads把你的結果得到Pi的一個「有效」近似（或以上N有#pragma omp for例如分發您的循環）

但是，這是如果/當你的代碼是正確的其他地方，它還沒有。由於BitTickler已經暗示，rand()不是線程安全的。所以你必須去另一個隨機數發生器，這將允許你私有化它的狀態。例如，這可能是rand_r()。儘管如此，這仍然有不少問題：

1. rand()/rand_r()是隨機性和週期性的任期可怕 RNG。在增加嘗試次數的同時，您將快速超過RNG的週期，並一遍又一遍地重複相同的順序。你需要更強大的東西來遠程執行任何嚴肅的事情。
2. 即使使用「良好」的RNG，並行性方面也可能是一個問題，因爲您希望並行的序列在彼此之間不相關。而只是用每個線程不同的種子值並不保證該給你（雖然寬足以RNG，你有一點餘量爲的）

無論如何，底線是：

• 使用更好的線程安全的RNG（我發現drand48_r()或random_r()是在Linux上的玩具代碼OK）
• 初始化其狀態每個線程例如基於線程ID，同時牢記這不會在某些情況下確保隨機序列的正確解相關（並且您調用函數的次數越多，則越多很可能你最終會有重疊系列）。

這個工作（有一些小的修正一起），你的代碼變得舉例如下：

#include <stdlib.h> 
#include <stdio.h> 
#include <time.h> 
#include <math.h> 
#include <omp.h> 

typedef struct drand48_data RNGstate; 

double sample_interval(double a, double b, RNGstate *state) { 
    double x; 
    drand48_r(state, &x); 
    return (b-a)*x + a; 
} 

int main (int argc, char **argv) { 

    int N = atoi(argv[1]); // convert command-line input to N = number of points 
    int NumThreads = 4; 
    const double pi = 3.141592653589793; 
    double x, y, z; 
    double counter = 0; 
    time_t ctime = time(NULL); 

    #pragma omp parallel private(x, y, z) reduction(+:counter) num_threads(NumThreads) 
    { 
     RNGstate state; 
     srand48_r(ctime+omp_get_thread_num(), &state); 
     for (int i=0; i < N; ++i) { 
      x = sample_interval(-1, 1, &state); 
      y = sample_interval(-1, 1, &state); 
      z = ((x*x)+(y*y)); 

      if (z<= 1) { 
       counter++; 
      } 
     } 

    } 
    double approx_pi = 4.0 * counter/(NumThreads * N); 

    printf("%i %1.6e %1.6e\n ", N, approx_pi, fabs(approx_pi - pi)/pi); 

    return 0; 
} 

其中我這樣進行編譯：

gcc -std=gnu99 -fopenmp -O3 -Wall pi.c -o pi_omp