C++11のthread safeな乱数生成器 (random number generator; RNG) のOpenMPでの使い方

00CXXThreadSafeRNG.md

C++11のthread safeな乱数生成器 (random number generator; RNG) のOpenMPでの使い方

C++11の乱数生成器 (random number generator; RNG) である<random>はとてもよくできている。 single threadで動くプログラムなら自作の珍妙なラッパを介さず直接使うべき。 thread safeなのでOpenMPからも使うことができる。 高速にかつ２カ所以上でRNGを使うことを前提にして、 ここgistに自作の珍妙なラッパを公開してみる。

指針

thread数ぶん作った乱数生成エンジンを２カ所 (main.cpp と timing.h) で使う。

乱数生成エンジン

• 乱数生成エンジンを thread local storage (TLS) として確保し、externで大域変数とする。

浮動小数点数一様分布乱数への変換

• std::uniform_real_distribution<double> は乱数生成エンジンを引数にとって、所望の精度（ここではdouble）の一様分布乱数を返す。 これもスレッド毎に確保するが、大域変数とはしない。
• 内部でstd::generate_canonicalを通している。 所望の精度（ここではdouble）に達するまで乱数生成エンジン (32 or 64 bit) が必要な回数だけ繰り返し呼び出される。
• doubleの一様分布乱数が欲しい場合は32 bitより64 bitの乱数生成エンジンを用いるほうが高速といわれている。

浮動小数点数の一様分布乱数の生成速度

• 2012年頃発売の2.5 GHz Intel Core i5で、倍精度なら約50M個／秒／core。
• 2018年頃発売のIntel Xeonで、倍精度なら約100M個／秒／core、単精度なら約300M個／秒／core。
• Hyper Threadingは乱数の発生だけならわりと有効に機能するが、たいていの場合において乱数を使う他の計算で十分には機能しない。
• <random>はもともとBoostの一部として開発されC++11の一部として取り込まれたので、Boostにもほぼ同じ擬似乱数生成器がある。 g++のversion 4か5まではBoostの生成器のほうが高速だった。g++ version 7以降はほぼ変わらない生成速度になっている。

(0.0, 1.0]の一様分布乱数

（書きかけ）

• 1-r
• std::uniform_real_distributionの引数で

(0.0, 1.0]の一様分布でない分布

• Boostには球面

ライセンス

Makefile

# -*-Makefile-*- for randomtest.cpp
##
CXX=g++
CXXFLAGS=-fopenmp -std=c++11 -O2 -g -Wall -pipe
LDFLAGS=

randomtest: randomtest.o rng_engines.o timing.o
	$(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@

%.pdf: %.cpp
	LANG=C a2ps --medium=a4 --header='Takeshi Nishimatsu' --underlay=$<\
	 --center-title=$< --prologue=color --portrait --columns=1\
	 --borders=off -f 11.0 --pretty-print=cpp -o - $< | ps2pdf -sPAPERSIZE=a4 - $@

clean:
	rm -rf *.o core *.s *.pdf randomtest.dSYM randomtest

randomtest.cpp

// randomtest.cpp for learning C++11 thread safe RNG
////
#include <omp.h>
#include <iostream>
#include <random>
#include "rng_engines.h"
#include "timing.h"

int main()
{
  const size_t N = 20;
  double* ary = new double[N];  // allocate them in heap
  int*    thd = new    int[N];

  init_rng_engines(12345678901234567890ULL);

# pragma omp parallel
  {
    std::uniform_real_distribution<double> uni0(0.0,1.0);
#   pragma omp for
    for (size_t i=0; i<N; i++) {
      thd[i] = omp_get_thread_num();
      ary[i] = uni0(*rng_engine_per_thread);
    }
  }

  for (size_t i=0; i<N; i++) {
    printf("%5zu %2d   %10.8f\n", i, thd[i], ary[i]);
  }

  timing();

  destruct_rng_engines();
  delete [] thd;
  delete [] ary;
  return 0;
}

rng_engines.cpp

// rng_engines.cpp
////
#include <omp.h>
#include <random>
#include "rng_engines.h"

std::mt19937_64 *rng_engine_per_thread;

void init_rng_engines(uint64_t seed)
{
  #pragma omp parallel
  {
    rng_engine_per_thread = new std::mt19937_64(seed+omp_get_thread_num());
  }
}

void destruct_rng_engines()
{
  #pragma omp parallel
  {
    delete rng_engine_per_thread;
  }
}

rng_engines.h

// rng_engine.h
////
#include <random>
#ifndef RNG_ENGINES_H
#define RNG_ENGINES_H
extern std::mt19937_64 *rng_engine_per_thread;
#pragma omp threadprivate(rng_engine_per_thread)
void init_rng_engines(uint64_t seed);
void destruct_rng_engines();
#endif

timing.cpp

// timing.cpp
////
#include <omp.h>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <random>
#include "rng_engines.h"

template<typename _Rep, typename _Period>
void report_timing(const std::chrono::duration<_Rep, _Period>& dur, const size_t nn)
{
  auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(dur).count();
  std::cout << std::setw( 5)
            << std::chrono::duration_cast<std::chrono::milliseconds>(dur).count()
            << " ms, " << std::setprecision(4)
            << (double)nn / ms / omp_get_max_threads() / 1000 << " M/s/thread" << std::endl;
}

void timing()
{
  const size_t NN = 2000000000ull;
  auto start = std::chrono::system_clock::now();
  double sum=0.0;
# pragma omp parallel
  {
    std::uniform_real_distribution<double> uni0(0.0,1.0);
#   pragma omp for reduction(+:sum) 
    for (size_t i=0; i<NN; i++) {
      sum += uni0(*rng_engine_per_thread);
    }
  }
  auto end = std::chrono::system_clock::now();
  report_timing(end-start, NN);
}

timing.h

// timing.h
////
#include <chrono>

#ifndef TIMING_H
#define TIMING_H
template<typename _Rep, typename _Period>
  void report_timing(const std::chrono::duration<_Rep, _Period>& dur, const size_t nn);
void timing();
#endif