炉物理の研究

第 55 号 (2003 年 3 月)

＜特集：炉物理研究への PC クラスタの利用・並列計算プログラミング超入門＞

並列計算プログラミング超入門
佐々木誠(株)日本総合研究所

sasaki.makoto@jri.co.jp

さて、ここまでの記事であなたの手元には PC クラスターが構築されているでしょう。ただ、そのままでは単なる PC をネットワークでつないだシステムにすぎません。これからこの上で「並列計算」を行なうソフトウェアを自ら構築するか、他所から導入するかするわ
けですが、そのためには並列処理を手助けするソフトウェアを導入する必要があります。
ここではそのようなソフトウェアとして MPI(Message Passing Interface)を導入し、あわせて簡単な例で MPI を用いて並列計算を行なうためにはどのようなプログラミングを行
なうかを示します。
１. MPI の導入
MPI(Message Passing Interface)は並列処理を行なうプログラミングを行なうためのプ
ログラミングインターフェイスの規格であり、米欧日を中心とした企業や機関を会員とす
る非営利機関である MPI Forum で議論され決定されています。現在、MPI-2.0 までの規格が決定されています(http://www.mpi-forum.org/)。
MPI Forum はプログラミングのためのインターフェイス、すなわち C 言語関数およびFORTRAN のサブルーチンや関数を決めているだけで、それらの関数ライブラリの構築や、どのようにして並列計算を行なわせるかについては別途それらの実装者を想定しています。
代 表 的 な 実 装 と し て

MPICH(http://www-unix.mcs.anl.gov/mpi/mpich/) やLAM(http://www.lam-mpi.org/)があります。MPICH は各種メーカー製並列計算機用ライ
ブラリーのベースにも採用されている等の実績があり利用者ベースが大きいと思われます
ので、ここでも MPICH を導入する事にして解説します。
●MPICH のコンパイルとインストール
クラスターの OS を Linux としたとき、RPM などのパッケージとして MPICH が提供されている場合もあると思いますが、ここではソースからコンパイルしてインストールする方法を解説します。まず MPICH の上記の WWW ページからソースをダウンロードします("Download MPICH"をクリックしてダウンロードのページに行きます)。ファイル名はmpich.tar.gz でこの記事の執筆中の最新バージョンは mpich-1.2.4 です。
さて mpich.tar.gz を入手したら適当な場所でこれを展開します；

tar zxvf mpich.tar.gz

39炉物理の研究第 55 号 (2003 年 3 月)

mpich-1.2.4 という名前のディレクトリが作られているはずですのでそこに移動してconfigure コマンドを実行します。

cd mpich-1.2.4./configure -prefix=/home/mpich
ここで"-prefix=/home/mpich"という「オプション」をつけたのは MPI がインストールされる先を現在の場所ではなくて違う場所にしたいときに有効です。特にそのインストール
先が(今の場合は/home/mpich)クラスターを構成するマシン間で NFS などによるファイル
共有されている場所であれば、クラスターの各マシンでいちいち MPICH のインストールを行なわなくてもどれかのマシンで一回だけインストールを行なえばよいということになります。
configure コマンドの実行で Makefile が作成されますのでコンパイルを行ないます。
makeコンパイルが終了したらインストールを行ないます(インストール先が root ユーザーでな
ければ書き込みが出来ないところなら su コマンドで root ユーザーになっておきます)。
make install
インストールが終了したら MPICH 関係のコマンドが実行出来るように環境変数
PATH(bash をログインシェルにしている場合)あるいはシェル変数 path(tcsh をログインシ
ェルにしている場合)に MPICH をインストールしたディレクトリ以下の bin ディレクト
リが含まれるようにする必要があります。ログインシェルに応じて以下の内容をログイン
シェルの初期化コマンドファイルに追加します。
・bash の場合： ファイル".bashrc"に以下を追加します。
PATH=${PATH}:/home/mpich/bin
export PATH・tcsh の場合： ファイル".cshrc"もしくは".tcshrc"に以下を追加します。
path=($path /home/mpich/bin)
上記で/home/mpich/bin は適当に貴方のインストール先に読み変えて下さい。
●MPICH の設定変更

40

炉物理の研究第 55 号 (2003 年 3 月)

MPICH をインストールした後で設定を変えておいた方がよいことがあります。インスト
ー ル し た 先 の デ ィ レ ク ト リ に share と い う デ ィ レ ク ト リ が あ り 、 そ の 中 に
machines.LINUX というファイルがあるはずです。そしてたとえばインストールを行なった PC のマシン名が"hogehoge1"であった場合その内容に以下のような行があるはずです。
hogehoge1
hogehoge1
hogehoge1
hogehoge1
hogehoge1
このままでは並列計算コマンドで計算を実行するとき、計算を行なうマシンが hogehoge1
ばかりになり、複数マシン並列計算になりません。もし PC クラスターが hogehoge1 から
hogehoge8 までの 8 台の PC で構成されている場合には以下のように書き換えておきます。
hogehoge1
hogehoge2
hogehoge3
hogehoge4
hogehoge5
hogehoge6
hogehoge7
hogehoge8●rsh、rlogin が実行できる環境にする。
MPICH が動作するためには rsh(remote shell)と rlogin(remote login)がクラスタを構成する PC 間でパスワード無しで実行出来る環境でなければなりません。そのような設定を可能にするには /etc/hosts.equiv にマシン名のリストを記述します(root ユーザ権限が必要です)。
●RedHat7.2 での rsh、rlogin の問題
各 PC にインストールされた Linux が RedHat7.2 以降の場合、デフォルトでは rsh、rloginが使えない設定になっています。これらのコマンドの実行を可能にするには/etc/xinetd.d/rsh と /etc/xinetd.d/rlogin に含まれる以下のような行を；
disable= yes以下のように書き換えておきます；
disable= no

このあと xinetd デーモンの再起動を行ないます；

41

炉物理の研究第 55 号 (2003 年 3 月)

/etc/rc.d/init.d/xinetd restart
なおこれらの処理をおこなうには root ユーザー権限が必要です。
同様のケースが RedHat Linux 以外にもあるかも知れませんので注意してください。
２.MPI による並列計算プログラミング(1)MPI では並列処理を行なうための様々な関数(サブルーティン)が用意されていて、並列計
算プログラミングを行なう人はそれらを利用したプログラムをコーディングします。まず簡単な例を FORTRAN で書いたものを以下に示します。これは MPICH のサンプルとしてついてくるものと同様のものですが、数値積分によって円周率πの計算を行ないます。
C============================================================
implicit double precision (a-h,o-z)
include 'mpif.h'Ccall MPI_Init(ierr)call MPI_Comm_Size(MPI_COMM_WORLD,nproc,ierr)call MPI_Comm_Rank(MPI_COMM_WORLD,myrank,ierr)C
if( myrank.eq.0 ) thenread(*,*) nwrite(*,*) 'Number of divisions ',nend ifC
call MPI_Bcast(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)C
t1 = MPI_Wtime()Ch = 1.0d0/n
pisum = 0
do i=myrank,n-1,nproc
x = (i+0.5)*hpisum = pisum + sqrt(1-x*x)
end dopisum = 4.0*h*pisumCcall MPI_Reduce(pisum,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,&
MPI_COMM_WORLD,ierr)C
t2 = MPI_Wtime()Cif( myrank.eq.0 ) thenwrite(*,'(a,1p,d18.10)') 'Calculated Pi is ',pi
write(*,*) 'Calculattion time is ',t2-t1
end ifC
call MPI_Finalize(ierr)C
end
C================================================================

42

炉物理の研究第 55 号 (2003 年 3 月)

このサンプルコードで使われている"MPI_"で始まるサブルーチンが MPI の FORTRAN
サブルーティンです。また include 'mpif.h' があることに注意して下さい。FORTRAN のMPI サブルーチンの引数に MPI_INTEGER、MPI_COMM_WORLD といった定数が必要になりますが mpif.h はこれらを定義しているファイルで、一般にはこれをインクルードします(C 言語の場合にも同様に'mpi.h'をインクルードします)。
このプログラムをコンパイルするには MPICH でインストールされるコマンド mpif77 を用います。たとえば上記のプログラムが ex1.f というソースファイルにあるのなら；

mpif77 ex1.f
一般的には。
mpif77 [options] sources
また C 言語のプログラム用にはコマンド mpicc、Fortran90 用には mpif90、C++用にはmpiCC があります。
コンパイルで得られたバイナリファイルを a.out とすると、これを４つのプロセスで並列計算させるには mpirun コマンドを用います；

mpirun -np 4 a.out
一般には；

mpirun -np

number_of_process program [arguments]

コンパイルで得られたバイナリファイルの存在するディレクトリとそれを実行するディ
レクトリが異なる場合には program はフルパス名で指定する必要があります。たとえば
a.out の フ ル パ ス 名 が /home/mpitest/a.out で 、 mpirun を 実 行 す る デ ィ レ ク ト リ が/home/mpitest でない場合には以下のようにしなければなりません；

mpirun -np 4 /home/mpitest/a.out
さてプログラムソースに戻って、なぜこのプログラムが並列処理になっているのかを各

43

炉物理の研究第 55 号 (2003 年 3 月)

MPI サブルーチンの説明を通して解説します。
◎ MPI_Init(ierr)
これが呼ばれることによって MPI でのメッセージ交換が出来るようになります。実際には mpirun によるプログラムの実行の段階で同じプログラムがすでに複数動いているのですが、MPI_Init()を呼び出すことで、動いている各プログラム(これらを「プロセス」とよぶことにします)が MPI のプロセスとしてメッセージの交換などができるようになります。
MPI での並列計算では必ず call します。
ierr はエラーコードで MPI の FORTRAN プログラムではほとんどのものがこれを最後の引数とします(C 言語の場合には関数の戻り値になります)。
◎ MPI_Comm_Size(communicator,nproc,ierr)
並列計算がいくつのプロセスで行なわれているかを知るために呼び出します。並列プロセ
ス数をこれで与えるわけではありません。並列プロセス数は mpirun コマンド呼び出し時に決まっています。
communicator という引数がありますが、これは並列プロセスのまとまりに対して与えられた識別番号のようなもので、各種の MPI 呼び出し関数で指定されます。大抵の場合ここ
で は デ フ ォ ル ト で 決 ま っ て い て mpirun で 生 成 す る 全 て の プ ロ セ ス を 含 ん で い る
MPI_COMM_WORLD を 使 用 し ま す 。 自 分 で プ ロ セ ス の 別 の 束 ね 方 を 指 定 し て
communicator を作るということもできますが、これはかなり凝った処理をする場合に必要でしょう。
◎ MPI_Comm_Rank(communicator,rank,ierr)
自分が communicator 中の何番目のプロセスであるかを得ます。番号は 0 番から順に付けられます。rank が 0 のプロセスは mpirun を実行したマシンのプロセスになるので、計算の入出力などを担当させることが多くなります。また計算の分担を決めるとき等にもrank を使用します。
◎ MPI_Bcast(buffer,number_of_data,data_type,source,communicator,ierr)
配列 buffer に含まれる number_of_data 個のデータを rank 値が source のプロセスから
communicator の全プロセスに同報します。すなわち送り手のプロセス source の配列 bufferの値が、そのプロセス以外のすべてのプロセスの配列 buffer にコピーされます。例では積44

炉物理の研究第 55 号 (2003 年 3 月)

分区間の分割数 n だけを同報していますので number_of_data は 1、また data_type は整数 型 な の で MPI_INTEGER と な り ま す 。 デ ー タ 型 を 指 定 す る 定 数 と し て は 、 他 に
MPI_REAL, MPI_DOUBLE_PRECISION な ど も あ り ま す (C 言 語 で は MPI_INT 、
MPI_FLOAT、MPI_DOUBLE などとなります)。
◎ MPI_Reduce(buffer1,buffer2,number_of_data,data_type,operation,dest,
communicator,ierr)
配列 buffer1 に含まれる各データに対して operation で指定される「縮約」操作をおこな
い、それをランクが dest の配列 buffer2 に返します(buffer1 と buffer2 にメモリ上で重なりがあってはいけません)。例では operation が MPI_SUM となっていて総和計算を行なっています。operation には他に MPI_MAX, MPI_MIN(最大値、最小値)、MPI_PROD(総乗)などがあります。
ここで例にあげたプログラムでは pisum という変数の全プロセス総和を rank が 0 のプロ
セスの pi という変数に入れています。この pisum という変数の計算のされ方を見てみましょう。通常なら；
pisum = 0
do i=0,n-1....
end doと、n 分割された積分区間の全てを計算しないと答えになりませんが、例では；
pisum = 0
do i=myrank,n-1,nproc....
end doとなっていて図 1 のようにプロセス毎にひとつずつずれた開始インデックスから nproc おきに計算していて、各プロセスの計算量は前の場合の 1/nproc になっています。そしてこれらの総和をとってはじめて積分が完了します。このため並列計算での各タスクの計算量が減って、全体としての計算時間が小さくなります。
もちろん、以下のような分割でも同じようなことになります；
nn = (n+nproc-1)/nproc
n1 = myrank*nn
n2 = n1+nn-1
if(myrank.eq.nproc-1) n2 = n-1
pisum = 0
do i=n1,n2....
end do45

第 55 号 (2003 年 3 月)

炉物理の研究

Y

O

0

図１

1 2 0

1

2 0

1 2

0 1
1

X

積分区間のプロセスへの分解（３プロセスの場合）

各プロセスでの計算値の総和を行なうことで答えが得られるものには、数値積分以外に、放射線輸送モンテカルロ計算などがあります。
◎ MPI_Wtime()
経過時間(CPU ではなくていわゆる wall clock time)を返します。MPI_Wtime()は MPIのライブラリーでは例外的に FORTRAN でも関数呼び出しになっています。
◎ MPI_Finalize(ierr)
これが呼ばれることによって MPI でのメッセージ交換が終了します。これ以降は
MPI_Reduce などの MPI 処理を行なうことはできません。
MPI のライブラリー関数(サブルーチン)はまだまだありますが、詳しくは MPI フォーラ
ムのホームページから入手できるドキュメント、あるいは Web ぺージを参照して下さい。
３.MPI による並列計算プログラミング(2)
前項で挙げたものよりもう少し複雑な例を挙げてみます。この例では以下のようなポアッ
ソン方程式を２次元の差分法で離散化し SOR 法で解くものです。

46

第 55 号 (2003 年 3 月)

炉物理の研究∇ 2φ (r ) = S (r )境界条件は図 2 のようにしています。２次元のメッシュは図 3 のように Y 方向に分割され
ます。
y

∂φ
=0
∂n
S=0

∂φ
=0
∂n

∂φ
=0
∂n

S=0.1

φ=0
図２

x

例題の境界条件とソース分布y

プロセス 3プロセス 2プロセス 1プロセス 0

図３x

並列プロセスへの分解（４プロセスの場合）

まずはメインルーチンと計算結果の出力ルーチンです。
implicit double precision (a-h,o-z)
include 'mpif.h'Cparameter (memsize = 500000)
common /marea/ a(memsize)C47

炉物理の研究第 55 号 (2003 年 3 月)

call MPI_Init(ierr)call MPI_Comm_Size(MPI_COMM_WORLD,nproc,ierr)call MPI_Comm_Rank(MPI_COMM_WORLD,myrank,ierr)C
if( myrank.eq.0 ) thenread(*,*) nx,ny
write(*,'(a,2i10)') 'Number of divisions ',nx,ny
end ifC
call MPI_Bcast(nx,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)call MPI_Bcast(ny,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)C
nn = (ny+nproc-1)/nproc
ny1 = myrank*nn + 1
ny2 = ny1 + nn - 1
if ( myrank.eq.nproc-1 ) ny2 = nyClphi = 1
lphi2 = lphi + nx*(ny2-ny1+1+2)
lmem = lphi2 + nx*(ny2-ny1+1)
if ( lmem > memsize+1 ) thenwrite(*,*) '== process ',myrank,' memory over ',lmem-1
stop 1
end ifC
call MPI_Barrier(MPI_COMM_WORLD,ierr)C
call poisson( nx,ny,ny1,ny2,myrank,nproc,a(lphi),a(lphi2))Cnymax = nnClphi3 = lphi2
lmem = lphi3 + nx*nymax
if ( lmem > memsize+1 ) thenwrite(*,*) '== process ',myrank,' memory over ',lmem-1
stop 1
end ifC
call output(nx,ny,ny1,ny2,myrank,nproc,&a(lphi),a(lphi3),nymax)
Ccall MPI_Finalize(ierr)
C
stop
end
C--------------------------------------------------------------subroutine output(nx,ny,ny1,ny2,myrank,nproc,phi,phi3,nymax)
include 'mpif.h'
double precision phi(nx,ny1-1:ny2+1)
double precision phi3(nx,nymax)
integer mpistat(MPI_STATUS_SIZE)
nn = (ny+nproc-1)/nproc
do ir = 0,nproc-1
iny1 = ir*nn + 1
iny2 = iny1 + nn - 1
if ( ir.eq.nproc-1 ) iny2 = ny
kk = iny2-iny1+1
if( ir.gt.0.and.ir.eq.myrank ) thencall MPI_Send(phi(1,ny1),nx*kk,MPI_DOUBLE_PRECISION,0,&123+ir,MPI_COMM_WORLD,ierr)
end ifif ( myrank.eq.0 ) then48

炉物理の研究第 55 号 (2003 年 3 月)

if ( ir.gt.0 ) thencall MPI_Recv(phi3(1,1),nx*kk,MPI_DOUBLE_PRECISION,ir,&
123+ir,MPI_COMM_WORLD,mpistat,ierr)else
do j=1,kk
do i=1,nx
phi3(i,j) = phi(i,ny1+j-1)
end doend doendif
do j = 1,kk
write(10,*) (phi3(i,j),i=1,nx)
end doendif
end doreturnend
ここで新しい関数 MPI_Barrier が登場しました。
◎ MPI_Barrier(communicator,ierr)communicator に所属するすべてのプロセスがこの関数を実行するまで待ち合わせを行います。すべてのプロセスが同時に何かの処理を始めなければいけないような場合に利用します。ここでは後続の poisson ルーチン内で実行時間測定を行うために、poisson ルーチンの開始がすべてのプロセスで同期するようにするために使用しています。
そ してソ ース を与え る関 数とポ アッ ソン方 程式を SOR 法 で解く 本体 サブル ーチン
poisson()です。
C--------------------------------------------------------------double precision function S(i,j,nx,ny)
C ... source
S = 0
if ( i.ge.nx/4.and.i.le.nx/2&.and. j.ge.ny/4.and.j.le.3*ny/4) thenS = 0.1
end ifreturnend
C--------------------------------------------------------------subroutine poisson( nx,ny,ny1,ny2,myrank,nproc,phi,phi2)
implicit double precision (a-h,o-z)
include 'mpif.h'
double precision phi(nx,ny1-1:ny2+1)
double precision phi2(nx,ny1:ny2)
integer mpistat(MPI_STATUS_SIZE)C... SOR acceleration factor ....
w = 1.9d0
w1 = 1.0d0 - w
C
t1 = MPI_WTime()lcount = 049炉物理の研究第 55 号 (2003 年 3 月)

100 continue
lcount = lcount + 1
C ..... send/recieve process boundary values
if ( ny1.ne.1 ) thencall MPI_Send(phi(1,ny1),nx,MPI_DOUBLE_PRECISION,myrank-1,&
1,MPI_COMM_WORLD,ierr)
end ifif ( ny2.ne.ny ) thencall MPI_Send(phi(1,ny2),nx,MPI_DOUBLE_PRECISION,myrank+1,&
2,MPI_COMM_WORLD,ierr)
end ifif ( ny1.ne.1 ) thencall MPI_Recv(phi(1,ny1-1),nx,MPI_DOUBLE_PRECISION,myrank-1,&
2,MPI_COMM_WORLD,mpistat,ierr)
end ifif ( ny2.ne.ny ) thencall MPI_Recv(phi(1,ny2+1),nx,MPI_DOUBLE_PRECISION,myrank+1,&
1,MPI_COMM_WORLD,mpistat,ierr)
end ifC
do j=ny1,ny2
do i=1,nx
phi2(i,j) = phi(i,j)
end doend doCC

if ( ny1.eq.1 ) then.... Dirichlet condition.
do i=1,nx
phi(i,1) = 0.0
end doelse
phi(1,ny1) = w1*phi(1,ny1) +&w*(phi(2,ny1) + phi(1,ny1-1) + phi(1,ny1+1)+&S(1,ny1,nx,ny))/3d0
do i=2,nx-1
phi(i,ny1) = w1*phi(i,ny1) +&w*(phi(i-1,ny1) + phi(i+1,ny1)&+phi(i,ny1-1) + phi(i,ny1+1)+&S(i,ny1,nx,ny))/4d0
end dophi(nx,ny1) = w1*phi(nx,ny1) +&w*(phi(nx-1,ny1) + phi(nx,ny1-1) + phi(nx,ny1+1) +&S(nx,ny1,nx,ny))/3d0
end ifC
do j = ny1+1,ny2-1
do i=1,nx
C
&
&
&
C
C
C

if ( i.ne.1.and.i.ne.nx ) thenphi(i,j) = w1*phi(i,j) +w *(phi(i-1,j) + phi(i+1,j)+phi(i,j-1) + phi(i,j+1) +S(i,j,nx,ny))/4d0

x = 1 : neumann conditionx = nx : neumann condition
else if(i.eq.1) thenphi(1,j) = w1*phi(1,j) +&w*(phi(2,j)+phi(1,j-1)+phi(1,j+1)+&
S(1,j,nx,ny))/3d050炉物理の研究第 55 号 (2003 年 3 月)&
else if(i.eq.nx) thenphi(nx,j) = w1*phi(nx,j) +&w*(phi(nx-1,j)+phi(nx,j-1)+phi(nx,j+1)+&
S(nx,j,nx,ny))/3d0
end ifend doend doCif ( ny2.eq.ny ) thenphi(1,ny) = w1*phi(1,ny) +&w*(phi(2,ny)+phi(1,ny-1)+S(1,ny,nx,ny))/2d0
do i=2,nx-1
phi(i,ny) = w1*phi(i,ny) +&w*(phi(i-1,ny)+phi(i+1,ny)+phi(i,ny-1)+&S(i,ny,nx,ny))/3d0
end dophi(nx,ny) = w1*phi(nx,ny) +&w*(phi(nx-1,ny)+phi(nx,ny-1)+S(nx,ny,nx,ny))/2d0
else
phi(1,ny2) = w1*phi(1,ny2) +&w*(phi(2,ny2) +phi(1,ny2-1) + phi(1,ny2+1)+&
S(1,ny2,nx,ny))/3d0
do i=2,nx-1
phi(i,ny2) = w1*phi(i,ny2) +&w*(phi(i-1,ny2) + phi(i+1,ny2)&+phi(i,ny2-1) + phi(i,ny2+1)+&
S(i,ny2,nx,ny))/4d0
end dophi(nx,ny2) = w1*phi(nx,ny2) +&w*(phi(nx-1,ny2)+phi(nx,ny2-1)+phi(nx,ny2+1)+&
S(nx,ny2,nx,ny))/3d0
end ifC
devmax = -1e30
do j=ny1,ny2
do i=1,nx
d = abs(phi2(i,j) - phi(i,j))
if ( d .gt. devmax ) devmax = d
end doend doC
call MPI_Allreduce(devmax,dd,1,MPI_DOUBLE_PRECISION,MPI_MAX,&MPI_COMM_WORLD, ierr)C
if ( dd.gt.1.0e-7.and.lcount.lt.100000 ) go to 100
C
t2 = MPI_WTime()
if ( myrank.eq.0 ) thenwrite(*,*) '== Converged : deviation max. ',devmax,&' iteration ',lcountwrite(*,*) '== Calculation time ',t2-t1
end ifreturnend

この例ではプロセス間の一対一通信を行なうルーティンが使用されています。

51

炉物理の研究第 55 号 (2003 年 3 月)

◎ MPI_Send(buffer,number_of_data,data_type,target,tag,communicator,ierr)
配列 buffer に含まれる number_of_data 個のデータを rank が target のプロセスに送ります。tag はメッセージを識別するのに使用するタグで、同じプロセスに複数のメッセージ
通信を続けて行う場合に各メッセージを受け側で対応するメッセージを正しく受け取れる
ように異なった値の tag を与えます。
◎ MPI_Recv(buffer,number_of_data,data_type,source,tag,communicator,
mpistatus,ierr)
配列 buffer に含まれる number_of_data 個のデータを rank が source のプロセスから受けとります。tag はメッセージを識別するのに使用するタグで、対応する MPI_Send と同じ値を指定します。mpistatus は FORTRAN では長さ MPI_STATUS_SIZE の配列で、正常にデータを受け取れたかどうかなどを示します。
◎ MPI_Allreduce(buffer1,buffer2,number_of_data,data_type,operation,communicator,ierr)
配列 buffer1 に含まれる各データに対して operation で指定される「縮約」操作をおこな
い配列 buffer2 に格納します。communicator の全てのプロセスが縮約計算の答えを受けとるのが MPI_Reduce と異なる点です。この例では oepration を MPI_MAX としてプロセス間の最大値を求めています。
MPI_Send と MPI_Recv はプロセス毎に分解されたメッシュ群の境界で値を交換するのに 使用されています。各プロセスは y が ny1 から ny2 までの計算を担当します。
phi(nx,ny1-1:ny2+1)のように y 方向にひとつづつ担当分メッシュより多いメッシュをもっていて、隣のプロセスの境界値を受け取れるようにしています。
ここで対応する MPI_Send と MPI_Recv において、MPI_Recv の方を先に call したら何が起こるでしょうか？MPI_Recv はデータを受け取り終えるまで待たなければならないので、２つの隣接するプロセス間で相手のデータがやってくるのをずっと待ち続ける状態になります。このような状態を「デッド・ロック」(dead lock)とよびます。並列処理ルーチ
ンを使用する際にはデッド・ロックのような状態が発生するのを注意深く避ける必要があります。この例の場合ではこのような問題に気を使わないですむ MPI_Sendrecv というデータ交換専門のルーチンもありますので、それを使うのも手でしょう。せっかくですので
この例での計算時間の例を表１に示します。この計算は Pentium4/2GHz の PC を52

炉物理の研究第 55 号 (2003 年 3 月)

100Base-TX のネットワークボードとスイッチングハブで接続したクラスターによるものです。メッシュ数は 300×300 です。
表１

プロセス
数

例題の Poisson ソルバーの計算時間

時間（秒）

収束までの

繰り返し１万回繰り返し数あたりの計算時間

1

368.96

54577

67.6

2

219.01

56200

38.97

3

172.70

57810

29.87

4

147.64

59416

24.84

MPI_Recv がデータを受け取り終えるまで待たなければならないことについて書きまし
たが、データ受け取り終了を「待たない」MPI_Irecv というサブルーチンもあります。これはデータ受け取りの終了を待つ MPI_Wait というサブルーチンと対になって使用されます。
◎ MPI_Irecv(buffer,number_of_data,data_type,source,tag,communicator,
mpistatus,request,ierr)
MPI_Recv と同様にデータを受け取るルーチンですが、データ受け取り終了を待たないですぐにサブルーティンを抜け出す点が異なります。request は request handle と呼ばれるもので MPI_Irecv によって与えられます。
◎ MPI_Wait(request,mpistatus,ierr)
MPI_Irecv(または MPI_Isend)の終了を待ちます。request は対応する MPI_Irecv などで設定されたハンドルです。
これらを使用した poisson ルーチンを以下に示します。
C--------------------------------------------------------------subroutine poisson( nx,ny,ny1,ny2,myrank,nproc,phi,phi2)
implicit double precision (a-h,o-z)
include 'mpif.h'
double precision phi(nx,ny1-1:ny2+1)
double precision phi2(nx,ny1:ny2)
integer mpistat(MPI_STATUS_SIZE)C... SOR acceleration factor ....
w = 1.9d0
w1 = 1.0d0 - w
C53炉物理の研究第 55 号 (2003 年 3 月)

t1 = MPI_WTime()lcount = 0
100 continue
lcount = lcount + 1
C ..... send/recieve process boundary values
if ( ny1.ne.1 ) thencall MPI_Irecv(phi(1,ny1-1),nx,MPI_DOUBLE_PRECISION,&
myrank-1,2,MPI_COMM_WORLD,mpistat,ireq1,ierr)
end ifif ( ny2.ne.ny ) thencall MPI_Irecv(phi(1,ny2+1),nx,MPI_DOUBLE_PRECISION,&
myrank+1,1,MPI_COMM_WORLD,mpistat,ireq2,ierr)
end ifif ( ny1.ne.1 ) thencall MPI_Send(phi(1,ny1),nx,MPI_DOUBLE_PRECISION,myrank-1,&
1,MPI_COMM_WORLD,ierr)
end ifif ( ny2.ne.ny ) thencall MPI_Send(phi(1,ny2),nx,MPI_DOUBLE_PRECISION,myrank+1,&
2,MPI_COMM_WORLD,ierr)
end ifC
do j=ny1,ny2
do i=1,nx
phi2(i,j) = phi(i,j)
end doend doCCC

if ( ny1.eq.1 ) then.... Dirichlet condition.
do i=1,nx
phi(i,1) = 0.0
end doend ifC
do j = ny1+1,ny2-1
do i=1,nx
C
&
&
&
C
C
C

if ( i.ne.1.and.i.ne.nx ) thenphi(i,j) = w1*phi(i,j) +w *(phi(i-1,j) + phi(i+1,j)+phi(i,j-1) + phi(i,j+1) +S(i,j,nx,ny))/4d0

x = 1 : neumann conditionx = nx : neumann condition
else if(i.eq.1) thenphi(1,j) = w1*phi(1,j) +&w*(phi(2,j)+phi(1,j-1)+phi(1,j+1)+&
S(1,j,nx,ny))/3d0&
else if(i.eq.nx) thenphi(nx,j) = w1*phi(nx,j) +&w*(phi(nx-1,j)+phi(nx,j-1)+phi(nx,j+1)+&
S(nx,j,nx,ny))/3d0
end ifend doend doCif ( ny2.eq.ny ) then54

炉物理の研究第 55 号 (2003 年 3 月)

phi(1,ny) = w1*phi(1,ny) +w*(phi(2,ny)+phi(1,ny-1)+S(1,ny,nx,ny))/2d0
do i=2,nx-1
phi(i,ny) = w1*phi(i,ny) +&w*(phi(i-1,ny)+phi(i+1,ny)+phi(i,ny-1)+&S(i,ny,nx,ny))/3d0
end dophi(nx,ny) = w1*phi(nx,ny) +&w*(phi(nx-1,ny)+phi(nx,ny-1)+S(nx,ny,nx,ny))/2d0
end if&Cif( ny1.ne.1 ) call MPI_Wait(ireq1,mpistat,ierr)
if( ny2.ne.ny ) call MPI_Wait(ireq2,mpistat,ierr)C
if( ny1.ne.1 ) thenphi(1,ny1) = w1*phi(1,ny1) +&w*(phi(2,ny1) + phi(1,ny1-1) + phi(1,ny1+1)+&S(1,ny1,nx,ny))/3d0
do i=2,nx-1
phi(i,ny1) = w1*phi(i,ny1) +&w*(phi(i-1,ny1) + phi(i+1,ny1)&+phi(i,ny1-1) + phi(i,ny1+1)+&S(i,ny1,nx,ny))/4d0
end dophi(nx,ny1) = w1*phi(nx,ny1) +&w*(phi(nx-1,ny1) + phi(nx,ny1-1) + phi(nx,ny1+1) +&S(nx,ny1,nx,ny))/3d0
end ifC
if ( ny2.ne.ny ) thenphi(1,ny2) = w1*phi(1,ny2) +&w*(phi(2,ny2) +phi(1,ny2-1) + phi(1,ny2+1)+&
S(1,ny2,nx,ny))/3d0
do i=2,nx-1
phi(i,ny2) = w1*phi(i,ny2) +&w*(phi(i-1,ny2) + phi(i+1,ny2)&+phi(i,ny2-1) + phi(i,ny2+1)+&
S(i,ny2,nx,ny))/4d0
end dophi(nx,ny2) = w1*phi(nx,ny2) +&w*(phi(nx-1,ny2)+phi(nx,ny2-1)+phi(nx,ny2+1)+&
S(nx,ny2,nx,ny))/3d0
end ifC
devmax = -1e30
do j=ny1,ny2
do i=1,nx
d = abs(phi2(i,j) - phi(i,j))
if ( d .gt. devmax ) devmax = d
end doend doC
call MPI_Allreduce(devmax,dd,1,MPI_DOUBLE_PRECISION,MPI_MAX,&MPI_COMM_WORLD, ierr)C
if ( dd.gt.1.0e-7.and.lcount.lt.100000 ) go to 100
C
t2 = MPI_WTime()
if ( myrank.eq.0 ) thenwrite(*,*) '== Converged : deviation max. ',devmax,&' iteration ',lcountwrite(*,*) '== Calculation time ',t2-t1

55

炉物理の研究第 55 号 (2003 年 3 月)

end ifreturnend

では何のために MPI_Irecv と MPI_Wait を用いたのでしょうか？この例では MPI_Irecvと MPI_Wait の間に、プロセス間境界メッシュ以外でのメッシュでの処理がおかれていることが分かります。つまり MPI_Irecv で指定されたデータ通信と、その通信にかかわるデータを参照しない計算を「同時に」行なうことで MPI_Recv を使用した場合に比べて計算時間を短縮できるということになります。計算時間を表２に示します。表１の計算と同じ環境によるものです。このような通信の方法を「非同期通信」または「ノンブロッキング
通信」とよびます。100Mb Ether ボードとスイッチングハブといった安価な通信手段を用いた PC クラスターの場合、通信にかかる時間の負荷は並列処理専用に設計された計算機や
Myrinet などの高速の通信手段を用いた場合に比べて大きなものになり、あまり並列計算による高速化の効果が得られない場合があります。しかし非同期通信のような技法で計算時間短縮が行なえる場合もありますので、可能ならば使用してみるベきでしょう。
今回挙げた SOR 法による計算の例では繰り返し計算のメッシュスイープの順序が変わってしまってプロセス数ごとに収束までの繰り返し数が変化してしまい、あまり時間短縮効
果はありませんでしたが、単位繰り返しあたりの秒数を見ると効果があることはわかるでしょう。
この例のように並列プロセス数が変わると収束までの繰り返し数が変わってしまうような繰り返し方法は並列化を行う際に問題となります。解法を SOR 法でなくヤコビ法にすれば繰り返し回数はプロセス数と無関係になり、もっと並列化の効果が上がったかもしれま
せ ん 。 実 際 に ヤ コ ビ 法 を 用 い た ベ ン チ マ ー ク 計 算 も あ り ま す（http://w3cic.riken.go.jp/HPC/HimenoBMT/）。
表２

例題の Poisson ソルバーの計算時間（非同期通信を用いた場合）

プロセス数

計算時間（秒）

収束までの繰り返し１万回繰り返し数あたりの計算時間

1

367.57

54577

67.35

2

213.76

59498

35.93

3

161.47

64311

25.11

4

133.79

69068

19.37

４. おわりに56

炉物理の研究第 55 号 (2003 年 3 月)

並列計算の入門ということで MPI を用いたプログラミングを紹介しました。並列プログ
ラミングに慣れていない方にはなにやら奇妙で面倒臭いプログラムの書き方に思えたかもしれません。並列処理専用に設計された計算機の場合にはいろいろ並列プログラミングの手段があるでしょうが、PC クラスターでは現在のところ MPI あるいは同様のメッセージ
通信ルーチンによる PVM などを用いた並列計算が当分は主流であると思われます。この文章が PC クラスター利用の一助となれば幸いです。
参考文献、ウエッブサイトなど(1) 湯浅太一、他編 「はじめての並列プログラミング」(共立出版)(2) MPI Forum： http://www.mpi-forum.org(3) MPICH A Portable MPI Implementation：
http://www-unix.mcs.anl.gov/mpi/mpich/(4) 青山幸也 「並列プログラミング虎の巻 MPI 版」(日本 IBM、非売品)(5) MPI Primer/Developing with LAM：
http://www.lam-mpi.org/download/files/lam61.nol.doc.pdf

57