38474


C++/c++11の環境で疑似乱数生成のメルセンヌ・ツイスターを並列化して動作させる方法をご教示下さい。  
boost/random.hppとopenmpの組合せで実装出来るのかとも考えいろいろ試したのですが、  
簡単には並列化できないようです。

[code]

    #include <boost/random.hpp>
    #include <omp.h>
    
    using namespace std;
    using namespace boost;
    mt19937 gen;
    uniform_real<> uniform_dst(0,1);//min,max
    variate_generator< mt19937,uniform_real<> > uniform(gen,uniform_dst);
    
    void test1(){
      uniform.engine().seed(123);
      int NN=1600000;
      int i;
      double x;
    #pragma omp parallel for private(i,x,gen)
      for(i=0;i<NN;++i){
        x=uniform();
        printf("%f in thread:%d\n",x,omp_get_thread_num());
      }
    }
    int main(){
        test1();
        return 0;
    }
    
[/code]

乱数生成器のmt19937をベースにしているのですが、openmpのループの中からこの乱数生成器を呼び出しても、乱数生成器自体は並列化されているわけではないようです。uniform()は個別のスレッドから呼び出されているようですが、乱数生成器の部分は共通して動作しており、並列化はされていないようです。この乱数生成器自体を並列化したいと考えています。

よろしくお願い致します。


CC BY-SA 3.0
2017-10-05T02:46:26.277
38446
2019-03-05T22:33:45.973
2019-03-05T22:33:45.973
4236
25645
question
c++ boost openmp
疑似乱数生成（メルセンヌ・ツイスター）の並列化


乱数呼び出し部分を排他制御するか、  
あるいはスレッドセーフなメルセンヌ・ツイスターの実装を利用してはどうでしょうか。

例えば、mt-liteはスレッドセーフでありマルチスレッドで使えるとあります。  
<http://mt-lite.sourceforge.net/index.html.ja>


CC BY-SA 3.0
2017-10-05T03:42:17.393
38453
2017-10-05T03:42:17.393
3655
38446
answer


OpenMP並列処理でワーカスレッド毎にオブジェクトを複製したい場合、`firstprivate`節にて実現できます。このケースでは、それぞれのスレッドにて同一擬似乱数系列が得られることに注意ください。

[code]

    #pragma omp parallel for private(i,x) firstprivate(uniform)
    for(i=0;i<NN;++i){
      ...
    }
    
[/code]


CC BY-SA 3.0
2017-10-05T05:42:52.533
38457
2017-10-05T05:52:16.670
2017-10-05T05:52:16.670
49
49
38446
answer


メルセンヌ・ツイスターについては、疑似乱数発生法に関する業績で"日本IBM科学賞"を受賞された松本 眞氏の[Mersenne Twister Home
Page](http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/mt.html)が詳しいです。参考になる情報を得られると思います。


CC BY-SA 3.0
2017-10-05T05:51:22.833
38458
2017-10-05T05:51:22.833
217
38446
answer


本質的でない`uniform_real`と`variate_generator`を取り除いて説明します。

[code]

    void test1() {
        printf("%s\n", __func__);
        std::mt19937 gen{ 123 };
        #pragma omp parallel for private(gen)
        for (int i = 0; i < 16; ++i) {
            unsigned int x = gen();
            printf("%u in thread:%d\n", x, omp_get_thread_num());
        }
    }
    
[/code]

メルセンヌ・ツイスターに限らず全ての疑似乱数生成器は

  * ある状態で得られる乱数は決まっている
  * ある状態から遷移する次の状態は決まっている
  * 最初の状態を決定する値をseedと呼ぶ

という特徴があります。  
ですので`private(gen)`を指定したとしても、得られる乱数はスレッド間で同じになりますし、遷移する次の状態も同じになります。結果的にスレッドごとに得られる乱数列は同一のものとなってしまいます。これを解消するにはまずはスレッド間排他を行うことです。

[code]

    void test2() {
        printf("%s\n", __func__);
        std::mt19937 gen{ 123 };
    #pragma omp parallel for
        for (int i = 0; i < 16; ++i) {
            unsigned int x;
    #pragma omp critical
            {
                x = gen();
            }
            printf("%u in thread:%d\n", x, omp_get_thread_num());
        }
    }
    
[/code]

ただし、排他のコストが高すぎるため並列化している意義が薄くなります。そこで排他するよりもスレッド毎に疑似乱数生成器を用意することです。その際、個々の疑似乱数生成器のseedは異なる値に設定する必要があります。

[code]

    void test3() {
        printf("%s\n", __func__);
    #pragma omp parallel
        {
            std::mt19937 gen{ std::random_device{}() };
    #pragma omp for
            for (int i = 0; i < 16; ++i) {
                unsigned int x = gen();
                printf("%u in thread:%d\n", x, omp_get_thread_num());
            }
        }
    }
    
[/code]

* * *

以上を踏まえて質問のコードは次のように書けます。

[code]

    #include <random>
    #include <omp.h>
    
    void test1() {
        int const NN = 1600000;
    #pragma omp parallel
        {
            std::mt19937 gen{ std::random_device{}() };
            std::uniform_real_distribution<> dst{ 0, 1 };
            auto uniform = [&] { return dst(gen); };
    #pragma omp for
            for (int i = 0; i < NN; ++i) {
                auto x = uniform();
                printf("%f in thread:%d\n", x, omp_get_thread_num());
            }
        }
    }
    int main() {
        test1();
        return 0;
    }
    
[/code]


CC BY-SA 3.0
2017-10-06T00:17:31.143
38474
2017-10-06T00:17:31.143
4236
38446
answer