OpenMP並列での擬似乱数の生成方法

出発点: シリアルコード

Boost.Randomを使うと、(0,1]の実数乱数を生成するコードは以下のように書ける。

serial.cpp

#include <boost/random.hpp>
#include <iostream>

int main(int argc, char **argv) {
  unsigned long seed = 1234;
  boost::minstd_rand seed_gen(seed);
  boost::mt19937 gen(seed_gen);

  int count = 16;
  boost::uniform_real<> dist;
  for (int i = 0; i < count; ++i) {
    std::cout << i << ' ' << dist(gen) << std::endl;
  }
}

boost::mt19937は、初期化の際に直接シード(整数)を渡すこともできるが、ここではseedで初期化したシーケンスseed_genを用いて初期化している。

OpenMPによる並列化

スレッド毎に擬似乱数発生機を準備し、それぞれをシード列の異なる部分を使って初期化を行う。

parallel.cpp

#include <boost/random.hpp>
#include <boost/shared_ptr.hpp>
#include <iostream>
#include <vector>

#ifdef _OPENMP
# include <omp.h>
#else
  int omp_get_max_threads() { return 1; }
  int omp_get_thread_num() { return 0; }
#endif

int main(int argc, char **argv) {
  int num_threads = omp_get_max_threads();
  std::cout << "# num_threads = " << num_threads << std::endl;
  
  unsigned long seed = 1234;
  boost::minstd_rand seed_gen(seed);
  std::vector<boost::shared_ptr<boost::mt19937> > gens(num_threads);
  #pragma omp parallel
  {
    int tid = omp_get_thread_num();
    for (int p = 0; p < num_threads; ++p) {
      if (p == tid)
        gens[p].reset(new boost::mt19937(seed_gen));
      #pragma omp barrier
    }
  }

  int count = 16;
  #pragma omp parallel
  {
    int tid = omp_get_thread_num();
    boost::mt19937& gen = *gens[tid];
    boost::uniform_real<> dist;
    for (int i = 0; i < count; ++i) {
      #pragma omp critical
      std::cout << tid << ' ' << i << ' ' << dist(gen) << std::endl;
    }
  }
}

gensはnum_threads個の乱数発生器へのポインタを保持する配列である。乱数発生機の初期化

  #pragma omp parallel
  {
    int tid = omp_get_thread_num();
    for (int p = 0; p < num_threads; ++p) {
      if (p == tid)
        gens[p].reset(new boost::mt19937(seed_gen));
      #pragma omp barrier
    }
  }

は単に

  for (int p = 0; p < num_threads; ++p) {
    gens[p].reset(new boost::mt19937(seed_gen));
  }

と書いても論理的には同じであり、実際、実行結果も等しくなるが、各乱数発生器を各スレッドの「近く」のメモリ上に配置するために、parallel領域の中で順番にオブジェクトをnewしている。(ファーストタッチ)

乱数生成器を使う際は、

  #pragma omp parallel
  {
    int tid = omp_get_thread_num();
    boost::mt19937& gen = *gens[tid];
    // do some work
  }

のようにリファレンスgenを定義すると、その後のコードがシリアルの場合と全く同じように書けて便利である。

なお、上のコード例parallel.cppでは、各スレッドからの出力が1行の中で混じり合ってしまわないように、criticalディレクティブを挿入しているが、実際の計算では不要であり、性能低下の原因となるので入れるべきではない。