mpi4pyを使ってみよう

MPIについて

MPI標準のリリース履歴

この標準に準拠する実装は多數あります。

mpi4py について

インストールと設定

mpi4py は次のようにインストールすることができます。

code: bash

$ pip install mpi4py

pipでインストールする場合、コンパイル環境が必要になる場合があります。

conda でインストールする方が楽です。

code: bash

$ conda install mpi4py

次の形式のホストファイルを用意しておきます。

code: hostfile

localhost slots=4

このホストファイルは４CPUのローカルマシンで実行する場合です。

mpi4py の概要

MPIプログラムの基本は次のものがあります。

ランク０はマスタープロセスとなる

C/C++ と mpi4py でのコードの記述の違いを見てましょう。

code: C/C++

int main(int argc, char **argv)

{

int n, myid, numproc, i;

MPI_Init(&argc, &argv);

MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

MPI_COMM_rank(MPI_COMM_WORLD, &myid);

printf("Hello World %d\n", myid);

MPI_Finalize();

}

code: mpi4py_hello.py

import sys

from mpi4py import MPI

size = MPI.COMM_WORLD.Get_size()　　# SPMD環境の規模を取得

rank = MPI.COMM_WORLD.Get_rank() # このプロセスのランク(番号)

name = MPI.Get_processor_name() # このプロセスのプロセス名

sys.stdout.write(

f"Hello, World! I am process {rank} of {size} on {name}\n"

)

code: bash

$ mpirun --host=hostfile -n 4 python3 hello.py

これは、"python3 hello.py" のコマンドを４つのプロセスで実行するということになります。

Pythonオブジェクトを送受信

バッファオブジェクトの送受信

補足説明：C型配列とF型配列

科学技術計算で良く使用されるプログラム言語に、FOTRANがあります。

このFOTRANの配列が実際のメモリ上にデータを配置する方法が、C言語とは異なります。FOTRAN形式（F型)配列、C形式(C型)配列と言われます。

https://gyazo.com/64068d1d10c612f26b30b1c508700241

出典: 東京大学理学部地球惑星物理学科 Fortran演習 (地球惑星物理学演習)

MPIプログラムでの通信パターン

MPIプログラムでの通信パターンにはつぎのものがあります。

send/recv

broadcast

scatter/gather/reduce

Send/Recv通信

Send/Recv通信は、ポイント・ツー・ポイント（Point-to-Point: P2P）通信とも言われ、２端点間で通信を行うことを表しています。

code: mpi4py_send_rev.py

from mpi4py import MPI

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

size = comm.Get_size()

# マスタープロセス

if rank == 0:

data = {'x': 1, 'y': 2.0}

# マスタープロセスは、すべてのワーカープロセスのランクを通して、

# ワーカープロセスにデータを送信する

for i in range(1, size):

comm.send(data, dest=i, tag=i)

print('Process {} sent data:'.format(rank), data)

# ワーカープロセス

else:

# 各ワーカープロセスはマスタープロセスからデータを受け取る

data = comm.recv(source=0, tag=rank)

print(f'Process {rank} received data: {data}')

ノンブロッキング通信

code:mpi4py_send_recv_nonblocking.py

from mpi4py import MPI

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

size = comm.Get_size()

# マスタープロセス

if rank == 0:

data = {'x': 1, 'y': 2.0}

# マスタープロセスは、すべてのワーカープロセスのランクを通して、

# ワーカープロセスにデータを送信する

for i in range(1, size):

req = comm.isend(data, dest=i, tag=i)

req.wait()

print('Process {} sent data:'.format(rank), data)

# ワーカープロセス

else:

# 各ワーカープロセスはマスタープロセスからデータを受け取る

req = comm.recv(source=0, tag=rank)

data = req.wait()

print(f'Process {rank} received data: {data}')

集合通信

bcast(): ブロードキャスト

https://gyazo.com/27877a2d4ec418bf345e3573cc00c3f9

code: mpi4py_broadcast.py

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

size = comm.Get_size()

if rank == 0:

data = np.arange(4.0)

else:

data = None

data = comm.bcast(data, root=0)

if rank == 0:

print(f'Process {rank} broadcast data: {data}')

else:

print(f'Process {rank} received data: {data}')

code: bash

% mpirun --hostfile hostfile -n 4 python mpi4py_broadcasst.py

NumPy array のブロードキャスト

code: mpi4py_numpy_bcast.py

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

if rank == 0:

data = np.arange(100, dtype='i')

else:

data = np.empty(100, dtype='i')

comm.Bcast(data, root=0)

for i in range(100):

scatter(): 分散

Scatterは配列を取得し、その連続するセクションをコミュニケーターのランク全体に分散します。

つまり、タスクを分割してサブタスクをプロセスに分散するすることができます。 ただし、プロセッサの数より多くの要素を配布することはできないことに注意してください。

https://gyazo.com/5c5f6805b625bb6f610c9c04704393e3

code: mpi4py_scatter.py

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

nprocs = comm.Get_size()

if rank == 0:

data = np.arange(15.0)

# 各サブタスクのサイズを決定

ave, res = divmod(data.size, nprocs)

# 各サブタスクの開始インデックスと終了インデックスを決定

# データを配列のリストに変換

else:

data = None

data = comm.scatter(data, root=0)

print(f'Process {rank} has data: {data}')

code: bash

$ mpirun --hostfile hostfile -n 4 python mpi4py_scatter.py

合計15個の要素があるため、最初の3つのプロセスのそれぞれが4つの要素を受け取り、最後のプロセスが3つの要素を受け取ったことがわかります。

NumPy array を 分散

code: mpi4py_numpy_scatter.py

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

size = comm.Get_size()

rank = comm.Get_rank()

sendbuf = None

if rank == 0:

recvbuf = np.empty(100, dtype='i')

comm.Scatter(sendbuf, recvbuf, root=0)

assert np.allclose(recvbuf, rank)

gather(): 収集

https://gyazo.com/d782e16cd108f1ad7e14a754a3ab069f

code: mpi4py_gather.py

from mpi4py import MPI

import time

import math

t0 = time.time()

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

nprocs = comm.Get_size()

nsteps = 10000000

dx = 1.0 / nsteps

if rank == 0:

# 各サブタスクのサイズを決定

ave, res = divmod(nsteps, nprocs)

# 各サブタスクの開始インデックスと終了インデックスを決定

# 開始インデックスと終了インデックスをデータに保存

else:

data = None

data = comm.scatter(data, root=0)

# 各プロセスに分配された円周率の計算を実行

partial_pi = 0.0

x = (i + 0.5) * dx

partial_pi += 4.0 / (1.0 + x * x)

partial_pi *= dx

partial_pi = comm.gather(partial_pi, root=0)

if rank == 0:

pi = sum(partial_pi)

walltime = time.time() - t0

print(f'pi computed in {walltime:.3f} sec')

err_pi = abs(pi - math.pi)

print(f'{pi}: error is {err_pi}')

code: bash

$ mpirun --hostfile hostfile -n 4 python mpi4py_gather.py

pi computed in 1.437 sec

3.1415926535896697: error is 1.234568003383174e-13

reduce(): 縮約

MPI reduce 操作は、各プロセスの配列から値を取り込み、マスタープロセスでそれらを単一の結果に削減します。 これは基本的に、各プロセスからマスタープロセスへのやや複雑な送信コマンドを実行し、マスタープロセスにreduce 操作を実行させるようなものです。 MPI reduceは、これらすべてを1つの簡潔なコマンドで実行することができます。

https://gyazo.com/1e915529fd4f2fadaae4e087545f4876

code: mpi4py_reduce.py

from mpi4py import MPI

import time

import math

t0 = time.time()

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

nprocs = comm.Get_size()

nsteps = 10000000

dx = 1.0 / nsteps

if rank == 0:

# 各サブタスクのサイズを決定

ave, res = divmod(nsteps, nprocs)

# 各サブタスクの開始インデックスと終了インデックスを決定

# 開始インデックスと終了インデックスをデータに保存

else:

data = None

data = comm.scatter(data, root=0)

# 各プロセスに分配された円周率の計算を実行

partial_pi = 0.0

x = (i + 0.5) * dx

partial_pi += 4.0 / (1.0 + x * x)

partial_pi *= dx

pi = comm.reduce(partial_pi, op=MPI.SUM, root=0)

if rank == 0:

walltime = time.time() - t0

print(f'pi computed in {walltime:.3f} sec')

err_pi = abs(pi - math.pi)

print(f'{pi}: error is {err_pi}')

NumPy array を収集

code: mpi4py_numpy_gather.py

rom mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

size = comm.Get_size()

rank = comm.Get_rank()

sendbuf = np.zeros(100, dtype='i') + rank

recvbuf = None

if rank == 0:

comm.Gather(sendbuf, recvbuf, root=0)

if rank == 0:

for i in range(size):

並列での行列 x 配列

code: mpi4py_numpy_matrix.py

from mpi4py import MPI

import numpy

def matvec(comm, A, x):

p = comm.Get_size()

xg = numpy.zeros(m*p, dtype='d')

y = numpy.dot(A, xg)

return y

MPI-IO

MPIプログラムで複数プロセスからファイルに書き出す場合は、これまでの単純な方法では、I/Oが集中してしまうことでI/O待ちが発生して、せっかく苦労して並列化を行っているのに性能を引き出すことができなくなります。特にHPC計算で性能を出すためには、バックエンドにある並列ファイルシステムの特性について理解しておかなければなりません。

MPI 標準に含まれるファイル I/O のための APIで、これを使用すると、それぞれのプロセスがファイルのどの位置からどの位置まで書き込むかを与えて、分散してI/Oを行うことができます。

HDF5や Parallel netCDF などのアプリケーション向けのI/Oライブラリから、MPI-IOで並列I/Oを行うこともできます。

MPI-IOの基本コンセプト

すべてのプロセスで個別にすることも、プロセス間で共有することもできます

ファイルハンドルを介してアクセス

ファイルへの効率的な非連続アクセスを可能にします

集団的I/O(Collective I/O): MPIはプロセスのリード/ライトを調整します

独立I/O(Independent I/O): MPIはファイルI/Oに何も調整はしない

基本の操作

ファイルのオープン/クローズ/削除

クローズは

fh = MPI.File.Open(comm, filename, amode=MODE_RDONLY, info=INFO_NULL)

fh.Close()

fh.Delete()

ファイルポインタ

ファイルポインタは、それぞれのMPIプロセスは個別に保持します。

共有ファイルポインタは、それぞれのMPIプロセスで共有されます。

fh.Seek(offset, whence=SEEK_SET)

fh.Seek_share(offset, whence=SEEK_SET)

リード

実行後はファイルポインタは移動します。スレッドセーフではありません。

fh.Read(buf, status)

他にも多数のリード関連メソッドがあります。

Read_all()、Read_all_begin()、Read_all_end()、Read_at()、Read_at_all()、Read_at_all_begin()、

Read_at_all_end()、Read_ordered()、Read_ordered_begin()、Read_ordered_end()、Read_shared()

ライト

実行後はファイルポインタは移動します。スレッドセーフではありません。

fh.Write(buf, status)

他にも多数のWrite関連メソッドがあります。

Write_all()、Write_all_begin()、Write_all_end()、Write_at()、Write_at_all()、Write_at_all_begin()、

Write_at_all_end()、Write_ordered()、Write_ordered_begin()、Write_ordered_end()、Write_shared()

モード

ファイルビュー

デフォルトでは、ファイルはバイトで構成されているものとして扱われ、プロセスはファイル内の任意のバイトにアクセス（リード/ライト）できます。ファイルビューは、ファイルのどの部分をI/O対象にするかを定義します。

fh.Set_view(disp=0, etype=None, filetype=None, datarep=None, info=INFO_NULL)

デフォルトはNoneで、メモリと同じ形式で保存

サンプル

NumPy array を並列で書き出す

code: mpi4py_numpy_io.py

from mpi4py import MPI

import numpy as np

amode = MPI.MODE_WRONLY|MPI.MODE_CREATE

comm = MPI.COMM_WORLD

fh = MPI.File.Open(comm, "./datafile.contig", amode)

buffer = np.empty(10, dtype=np.int)

offset = comm.Get_rank()*buffer.nbytes

fh.Write_at_all(offset, buffer)

fh.Close()

NumPy array とデータを非連続に並列で書き出す

code: mpi4py_numpay_noncontiguous_io.py

from mpi4py import MPI

import numpy as np

comm = MPI.COMM_WORLD

rank = comm.Get_rank()

size = comm.Get_size()

amode = MPI.MODE_WRONLY|MPI.MODE_CREATE

fh = MPI.File.Open(comm, "./datafile.noncontig", amode)

item_count = 10

buffer = np.empty(item_count, dtype='i')

filetype = MPI.INT.Create_vector(item_count, 1, size)

filetype.Commit()

displacement = MPI.INT.Get_size()*rank

fh.Set_view(displacement, filetype=filetype)

fh.Write_all(buffer)

filetype.Free()

fh.Close()

SWIGとの連携

SWIG は C/C++、FORTRAN、Ruby、Perl など他の言語で記述したコードを相互に取り込むための仕組みを提供しています。よく使われるケースは、数値計算などの求解するときに速度の問題でアルゴリズムをC/C++、FORTRANで実装したり、過去の資産としてC/C++、FORTRANで記述されたライブラリがあり、それを呼び出したいといったものがあります。

今、C言語で記述されたコードがあるとします。

code: helloworld.c

/* file: helloworld.c */

void sayhello(MPI_Comm comm)

{

int size, rank;

MPI_Comm_size(comm, &size);

MPI_Comm_rank(comm, &rank);

printf("Hello, World! "

"I am process %d of %d.\n",

rank, size);

}

SWIGのインタフェースファイルを次のように記述します。

code: helloworld.i

// file: helloworld.i

%module helloworld

%{

}%

%include mpi4py/mpi4py.i

%mpi4py_typemap(Comm, MPI_Comm);

void sayhello(MPI_Comm comm);

次のようにして、Pythonから利用することができます。

code: python

>> from mpi4py import MPI

>> import helloworld

>> helloworld.sayhello(MPI.COMM_WORLD)

Hello, World! I am process 0 of 1.

mpi4py.futures

mpi4py.futures は、プロセス間通信にMPIを使用して、ワーカープロセスのプールで呼び出し可能オブジェクトを非同期的に実行するための高レベルのインターフェイスを提供します。

MPIPoolExecutor

注意

MPICommExecutor

従来のMPI-1または部分的なMPI-2実装における動的プロセス管理機能の欠如を軽減することに加えて、

MPICommExecutorコンテキストマネージャーは、mpi4py.futuresパッケージで利用可能なシンプルなタスクベースのマスター/ワーカーアプローチを利用することをいとわない古典的なMPIベースのPythonアプリケーションで役立つ場合があります。

mpiexecコマンド

MPI実装での動的プロセス管理機能のサポートの欠如に関連する問題を想定して、mpi4py.futuresは、Pythonコードがmpi4pyのコマンドラインから実行される代替の使用パターンをサポートします。

もっと簡単に言うと、mpi4py.futures は次の方法で呼び出すことができます。

code: bash

mpi4py.futures の例

code: mpi4py_julia.py

from mpi4py.futures import MPIPoolExecutor

x0, x1, w = -2.0, +2.0, 640*2

y0, y1, h = -1.5, +1.5, 480*2

dx = (x1 - x0) / w

dy = (y1 - y0) / h

c = complex(0, 0.65)

def julia(x, y):

z = complex(x, y)

n = 255

while abs(z) < 3 and n > 1:

z = z**2 + c

n -= 1

return n

def julia_line(k):

line = bytearray(w)

y = y1 - k * dy

for j in range(w):

x = x0 + j * dx

return line

if __name__ == '__main__':

with MPIPoolExecutor() as executor:

image = executor.map(julia_line, range(h))

with open('julia.pgm', 'wb') as f:

f.write(b'P5 %d %d %d\n' % (w, h, 255))

for line in image:

f.write(line)

mpiexec コマンドを使って、このスクリプトを実行します。

code: bash

$ mpiexec -n 1 -usize 17 python mpi4py_julia.py

mpiexecコマンドは、Pythonインタープリターを実行してメインスクリプトを実行する単一のMPIプロセス（マスター）を起動します。 必要に応じて、mpi4py.futuresは16個の追加のMPIプロセス（子プロセス）を生成して、ワーカーのプールを動的に割り当てます。 マスターは子プロセスにタスクを提出し、結果を待ちます。 子プロセスはタスクを受け取り、それらを実行して、結果をマスターに送り返します。

code: bash

$ mpiexec -n 17 python -m mpi4py.futures mpi4py_julia.py

前に説明したように、17のプロセスは1つのマスターと16のワーカーに分割されます。 マスタープロセスはメインスクリプトを実行し、ワーカーはマスターから送信されたタスクを実行します。

mpi4py.run

バージョン3.0.0の新機能。

残念ながら、自動MPIファイナライズと未処理の例外の相互作用により、デッドロックが発生する可能性があります。バッチ処理で実行すると、これらのデッドロックによってラップトップのバッテリーが消耗したり、スーパーコンピューティング施設で貴重な割り当て時間が消費されたりします。

次のPythonコードのスニペットについて考えてみます。このコードが標準のPythonスクリプトファイルに保存され、2つ以上のプロセスでmpiexecを使用して実行されると想定します。

code: python

from mpi4py import MPI

assert MPI.COMM_WORLD.Get_size() > 1

rank = MPI.COMM_WORLD.Get_rank()

if rank == 0:

1/0

MPI.COMM_WORLD.send(None, dest=1, tag=42)

elif rank == 1:

MPI.COMM_WORLD.recv(source=0, tag=42)

インターフェイスオプション

code: bash

まとめ

mpi4py を利用するとMPIを利用した並列計算をPython で簡単に記述、実行できるようになります。

ただし、mpi4py は MPI をラップしているだけなので、実際のプログラミングではMPIの知識が必要になってしまいます。また、実際の動作ではバックエンドのMPI実装の知識も必要になります。

それでも、マルチノードでの並列処理を簡単にできることは特筆するべきだと考えています。

参考