概要
数値計算をやってる人だったら「モンテカルロ法」の名前くらいは聞いたことがあると思います。分子のふるまいとかそういう「真面目にプログラム組むと結構めんどくさい」シミュレーションに威力を発揮するそうです。
モンテカルロ法自体は並列計算向き(ほかのCPUとデータのやりとりをあまりせず黙々と計算できる)ですが、問題は計算に使う乱数の発生をどうやって並列化するかです。
インターネットで検索するとこの分野については論文がたくさん見つかります。色々工夫するところがあるということだと思います。
本稿では、乱数発生の並列化にあたっての問題点と初歩的なMPI並列化の方法を検討します。
なお、乱数の発生のしかたそのものについては今回は深入りせず、シンプルな乗算合同式法で説明します。
乱数発生の並列化にあたっての問題点
Cの標準ライブラリなどで使用されている乱数は「線形合同法」といい、次の式で乱数を求めます。
X(N+1) = (A * X(N) + B) mod M
乗算合同式法はこの式のBが0の場合にあたります。
X(N+1) = (A * X(N)) mod M
Park & Millerの式ではA=48271、M=2**31-1(32ビットの正の最大の整数)となります。
これらの式からわかるように、乱数の計算式は1ステップ前で求めた乱数の値をもとに次の乱数を求める漸化式になっており、並列計算機が最も苦手とするパターンです(データ依存性がある)。
MPI並列でプログラムを組む場合であれば、1ステップ前の乱数の値を持っているプロセスとその都度通信して値をもらってくる必要がありますが、計算に比べて通信は圧倒的にコストが高いので実用的でないことは明らかです。
どうやって並列化するか
MPI並列において乱数発生を並列化する場合に、大きく分けて2通りのアプローチがありそうです。
- 乱数を発生するプロセスを決めておき、計算が一巡するごとに残りのプロセスに乱数をブロードキャストする
- 乱数の種をあらかじめランク数分のサイズの配列として用意しておき、各プロセスでは自分のランクのところの種を使って計算する
最初のアプローチでは乱数担当のランクをあらかじめ決めておき、そのランクで残りのプロセスが必要とする乱数(計算1回分)を求め、MPI_BCAST()で全プロセスに対してその乱数を放送します。
INCLUDE 'mpif.h'
DOUBLE PRECISION TSTART,TEND,TD
DIMENSION XRAND(1024)
DIMENSION ISTAT(MPI_STATUS_SIZE)
N = 20000000
M = 0
IA = 48828125
IR = 584287
MTTL = 0
C INITIALIZE MPI
CALL MPI_INIT(IERR)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR)
IF(IRANK.NE.0) GO TO 5
WRITE(6,100) ISIZE,N
100 FORMAT(1H ,13HPROCESSORS : ,16X,I4/1H ,13HN : ,
18X,I12)
C MAIN LOOP
5 TSTART = MPI_WTIME()
DO 10 I=IRANK+1,N,ISIZE
IF(IRANK.NE.0) GO TO 20
CALL GENRND(ISIZE,IA,IR,XRAND,2)
20 CALL MPI_BCAST(XRAND,2*ISIZE,MPI_REAL,0,MPI_COMM_WORLD,
1IERR)
X = XRAND(IRANK*2+1)
Y = XRAND(IRANK*2+2)
Z = X*X+Y*Y
IF(Z.LE.1.0) M = M + 1
10 CONTINUE
TEND = MPI_WTIME()
CALL MPI_REDUCE(M,MTTL,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,
1IERR)
IF(IRANK.NE.0) GO TO 30
PI = 4.0*REAL(MTTL)/REAL(N)
TD = TEND - TSTART
FLOPS = 7.0 * REAL(N) / TD
WRITE(6,110) PI,TD,FLOPS/1.0E6
110 FORMAT(1H ,13HPI : ,1PE20.12,
1/1H ,13HTIME(SEC) : ,10X,0PF10.3,
2/1H ,13HMFLOPS : ,10X,F10.3)
30 CALL MPI_FINALIZE(IERR)
STOP
END
SUBROUTINE GENRND(ISIZE,IA,IR,DEST,N)
DIMENSION DEST(ISIZE*N)
ML = (2**30-1)+2**30
BNORM = 1.0/(REAL(ML)+1.0)
DO 10 I=1,ISIZE*N
IR = IA*IR
IF(IR.LT.0) IR=IR+ML+1
DEST(I) = BNORM * REAL(IR)
10 CONTINUE
RETURN
END
MPI_BCASTは全プロセスでこのサブルーチンを呼び出したのを待って全プロセスにルートプロセス(この場合はランク0)の持っているデータを放送しますので、全プロセスで呼び出す必要があります。(呼び出さないプロセスがあるとそこで計算が止まります)
このプログラムを実行すると、ループ1回ごとに集団通信(MPI_BCAST)を行うため、実行時間が非常にかかります。しかも、完走せずに途中で止まることもしばしばあります。
MPI_BCASTはコミュニケータに属する全プロセスにデータを放送する集団通信関数で、本来ならば全プロセスがここで同期して次のループに進むという実行形態になるので「完走せず止まる」ということはないはずなのですが....OpenMPIのバグでしょうか。
この方法ではうまくない(たとえ完走したとしても実行時間が非常に長くなり、しかも並列数が増すほど遅くなる可能性が大)ので、第2の方法を試してみます。
C MONTE_PI_MPI.F
C CALCULATE PI BY MONTE CARLO METHOD
INCLUDE 'mpif.h'
DIMENSION IA(1024)
DIMENSION XRAND(2)
DOUBLE PRECISION TSTART,TEND,TD
N = 2000000000
M = 0
IR = 584287
MTTL = 0
C INITILAIZE MPI
CALL MPI_INIT(IERR)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,ISIZE,IERR)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,IRANK,IERR)
IF(IRANK.NE.0) GO TO 5
WRITE(6,100) ISIZE,N
100 FORMAT(1H ,13HPROCESSORS : ,16X,I4/1H ,13HN : ,
18X,I12)
C INITIALIZE RANDOM NUMBER
5 CALL INIRND(ISIZE,IA)
TSTART = MPI_WTIME()
DO 10 I=IRANK+1,N,ISIZE
CALL VRAND(ISIZE,IRANK,IA,IR,XRAND,2)
X = XRAND(1)
Y = XRAND(2)
Z = X*X+Y*Y
IF(Z.LE.1.0) M = M+1
10 CONTINUE
TEND = MPI_WTIME()
CALL MPI_REDUCE(M,MTTL,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,
1IERR)
IF(IRANK.NE.0) GO TO 40
PI = 4.0*REAL(MTTL)/REAL(N)
TD = TEND - TSTART
FLOPS = 7.0 * REAL(N) / TD
WRITE(6,110) PI,TD,FLOPS/1.0E6
110 FORMAT(1H ,13HPI : ,1PE20.12
1/1H ,13HTIME(SEC) : ,10X,0PF10.3,
2/1H ,13HMFLOPS : ,10X,F10.3)
40 CALL MPI_FINALIZE(IERR)
STOP
END
SUBROUTINE INIRND(ISIZE,IA)
DIMENSION IA(ISIZE)
ML = (2**30-1)+2**30
IA(1) = 48828125
DO 10 I=2,ISIZE
IA(I) = 48828125*IA(I-1)
IF(IA(I).LT.0) IA(I)=(IA(I)+ML)+1
10 CONTINUE
RETURN
END
SUBROUTINE VRAND(ISIZE,IRANK,IA,IR,DEST,N)
DIMENSION IA(ISIZE)
DIMENSION DEST(N)
ML = (2**30-1)+2**30
BNORM = 1.0/(FLOAT(ML)+1.0)
DO 10 I=1,N
IR = IA(IRANK+1)*IR
IF(IR.LT.0) IR=IR+ML+1
DEST(I) = BNORM * REAL(IR)
10 CONTINUE
RETURN
END
このプログラムでは、シミュレーションのループに入る前にINIRND()サブルーチンを呼び出し、並列実行のランク数だけ乱数の乗数(冒頭の式におけるAの値)を発生し、これを配列IAにストアします。
MPIではすべてのプロセスが同一のプログラムを実行するので、別にランク0で値を生成してMPI_BCAST()で放送しなくても、すべてのプロセスで同一の値を発生させることができます。実際に各プロセスで使用するのは自分のランクのところの値だけで他の部分は使わないので、少々もったいないような気もしますが、ここの部分はさほどメモリを喰わないので気にしないこととします。
気になる場合はランク0でAの値の配列を生成し、MPI_SEND()で各ランクに配信するようにすればよいでしょう。
シミュレーションのループでは、配列IAの自分のランクのところから乗数を取り出して乱数を求めます(VRAND()サブルーチン)。
実行結果
第2のアプローチで乱数を発生させ、モンテカルロ法によりπの値を計算したときの実行結果です。20億回ループを回し、コア数を変えて実行時間がどう変わるかを見ました。
なお、計算に使用したスパコンは以前製作したNanoPi-NEO4x3台並列のSUGIBAC S-820です(18コア)。
pi@S820-1:~$ mpiexec -np 1 --hostfile ./nodelist ./monte_pi_mpi
PROCESSORS : 1
N : 2000000000
PI : 3.141581773758E+00
TIME(SEC) : 41.279
MFLOPS : 339.156
pi@S820-1:~$ mpiexec -np 2 --hostfile ./nodelist ./monte_pi_mpi
PROCESSORS : 2
N : 2000000000
PI : 3.141594648361E+00
TIME(SEC) : 20.650
MFLOPS : 677.958
pi@S820-1:~$ mpiexec -np 6 --hostfile ./nodelist ./monte_pi_mpi
PROCESSORS : 6
N : 2000000000
PI : 3.141594648361E+00
TIME(SEC) : 2.109
MFLOPS : 6637.961
pi@S820-1:~$ mpiexec -np 12 --hostfile ./nodelist ./monte_pi_mpi
PROCESSORS : 12
N : 2000000000
PI : 3.141606092453E+00
TIME(SEC) : 1.139
MFLOPS : 12293.441
pi@S820-1:~$ mpiexec -np 18 --hostfile ./nodelist ./monte_pi_mpi
PROCESSORS : 18
N : 2000000000
PI : 3.141597509384E+00
TIME(SEC) : 0.743
MFLOPS : 18844.020
1コアで実行したときは約40秒かかっていますが、2コアで約20秒、6コアで約2秒....とコア数なりに速くなっています。18コアで実行したときはわずか0.75秒という驚きのスピードです。
所要時間がきれいにコア数に反比例していないのはNanoPi-NEO4のクロックが負荷に応じて自動的に最適な値に調節される(負荷が低いときはスピードを落として電力消費を抑える)ためです。18コアでの実行時はクロックが最大になる前に計算が終わってしまうこともあります。
実際に18コアで並列実行したとき、最大では上の例の倍のスピードが出ることがあります(約36GFLOPS)。理論性能約50GFLOPSのマシンですから、良好な実行能率といえそうです。
まとめ
並列計算で乱数(乗算合同式法)を用いるときは、ループ1周ごとにランク0で乱数を発生して全プロセスにMPI_BCAST()で放送するより、あらかじめランク数だけ乗数Aの値を配列で用意しておき、あとは各ランクでその値を使って独立に計算させたほうが通信量が少なく高速に計算できます。
なお、今回は乱数の質は無視して最も素朴な乗算合同式法で試しましたが、ほかの乱数でも同様の手法が使えると思います。