

お世話になります。  
現在モンテカルロ法を高速化しようとOpen MPを使って並列化処理を行なっているのですが、Open
MPを使ったときと使わなかった時を比べると、実行速度は若干速くなり、逆に結果の精度は悪くなりました。  
この原因は何でしょうか？わかる方教えてください。よろしくお願いいたします。

下記プログラムはモンテカルロ法を使って円周率を求めるものです。

【実行結果】  
・openMPなしの場合  
3.141556216 120948

・openMPありの場合（同じコードで複数回実行）  
3.102521972 103864  
3.104310832 108247  
3.106061276 109388

【実行環境】  
g++ (clangではなくGNUの方です）  
CPUは2コア4スレッドです。

[code]

    #include <iostream>
    #include <random>
    #include <chrono>
    #include <iomanip>
    #include <omp.h> //openMPを使わない時はコメントアウトします
    
    using namespace std;
    const long int N=1'000'000'000;
    
    random_device seed_gen;
    std::mt19937 engine(seed_gen());
    std::uniform_real_distribution<double> dist(0.0,1.0);
    
    
    int main(void)
    {
        double x,y,val=0.0;
        long int count=0;
        chrono::system_clock::time_point start, end;
    
        start = chrono::system_clock::now();
    
        #pragma omp parallel for reduction(+:count) //同様
        for (long int loop=0; loop < N; ++loop)
        {
          x=dist(engine);
          y=dist(engine);
          if(x*x+y*y<=1.) count++;
        }
         end = std::chrono::system_clock::now();
         double elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
        val=4*double(count)/N;
    
      cout<<setprecision(10)<<val<<" "<<elapsed<<endl;
    
      return 0;
    }
    
[/code]


CC BY-SA 4.0
2019-03-05T17:01:04.870
53216
2019-03-05T22:33:27.653
2019-03-05T22:33:27.653
4236
19263
question
c++ openmp
Open MPで並列化したら精度が落ちるモンテカルロ法


[疑似乱数生成（メルセンヌ・ツイスター）の並列化](https://ja.stackoverflow.com/q/38446/4236)とほぼ同じ問題を抱えています。

OpenMPの並列ループ内で乱数`engine`を参照しています。マルチスレッド安全ではない気がします。最善でもスレッド間で同じ乱数値を取得することになりその時点で精度は１スレッド分程度に落ちます。つまり８並列であれば`N/8`回の施行しか意味を持たなくなる。最悪な場合、`std::mt19937`の内部状態がスレッド間で不適切に更新され、乱数として機能しないもしくは乱数の精度が極端に落ちることも考えられます。

ところで、C言語及びC++言語ではBasic風に変数宣言を関数先頭に列挙すべきではありません。特にC++言語では宣言個所でコンストラクターが動作するため不適切です（更にはデストラクターの実行順が不要になった順ではなく、宣言と逆順に固定される点も不適切）。更に今回の場合、変数`x`、`y`はループ内で並列に更新するわけですから、ループ内で宣言すべきです。

[code]

    #pragma omp parallel
    {
        std::mt19937 engine{ std::random_device{}() };
        std::uniform_real_distribution<double> dist{ 0.0, 1.0 };
    #pragma omp for reduction(+:count)
        for (long int loop=0; loop < N; ++loop) {
           auto x=dist(engine);
           auto y=dist(engine);
           if(x*x+y*y<=1.) count++;
        }
    }
    
[/code]


CC BY-SA 4.0
2019-03-05T22:32:38.903
53219
2019-03-05T22:32:38.903
4236
53216
answer