openmpで並列して乱数生成する

はじめに

確率微分方程式を解くときに並列計算ができるときがあります。毎ステップ、乱数の生成を行うわけですが、単純に並列化してしまうと乱数のseedが並列したスレッドで同じになってしまうので正確な計算ができなくなります。
そこで、再現性を持たせてシミュレーションするために乱数のseedの固定法について書きます。

並列化について

ここではOpenmpでの並列化を考えます。
おさらいとして簡単に並列化を確認します。

openmp.c

void func(int i, int a[]){
    int j;
    for(j=0; j<M*i; j++){
        a[i]++;
    }
}

int main(void){
    int i = 0;
    int a[N] = {};
    int j = 0;
    #pragma omp parallel for
    for(i=0; i<N; i++){
        func(i, a);
        printf("thread: %d, a[%d]=%d\n", omp_get_thread_num(), i, a[i]);
    }

    return 0;
}

意味のない計算ですが、これでfunc()の計算は並列化できます。

乱数生成の並列化

乱数を生成する際にははじめにseedを固定しますが、単純にopenmpで並列化すると、今の変数がスレッドごとにコピーされるので全く同じ挙動をしてしまうはずです。
本家のsfmtを入れると簡単にseedの固定ができます。

以下のようなcodeで行うと、スレッドごとに乱数生成のためが固定されて、schedule(static)でiがどれかのスレッドに固定されて、毎回実行されることになります。

random.c

#include<stdio.h>
#include<omp.h>
#include"sfmt.h"

#define N 8

int main(void){
    int i=0, t=0;
    sfmt_t sfmt[N];
    #pragma omp parallel for schedule(static)
    for(i=0; i<omp_get_max_threads(); i++){
        sfmt_init_gen_rand(sfmt+i, i%4);
    }
    for(t=0; t<3; t++){
        printf("step: %d\n", t);
        #pragma omp parallel for schedule(static)
        for(i=0; i<N; i++){
            printf("thread: %d, num: %d, rand: %f\n", omp_get_thread_num(), i, sfmt_genrand_real3(sfmt+omp_get_thread_num()));
        }
    }
    return 0;
}