1

並列プログラミング超入門講習会
MPIコース
九州大学情報基盤研究開発センター

2

並列計算機の構成
計算ノード

ネットワーク

CPU

…

コア
メモリ

アクセラレータ（GPU等）

例： スーパーコンピュータシステム ITO サブシステム B
ノード数

128

CPU数 / ノード

2

コア数 / CPU

GPU数 / ノード

18

4

3

MPI (Message Passing Interface)による並列計算• 主に複数の計算ノードによるクラスタ型計算機向け• 複数のメモリ領域の間で「通信」をしながら計算• ノード内をさらに複数のメモリ領域に分けることも可能• ノード内で通信が必要となるが, その方が高速な場合がある

各ノードに一つずつ
メモリ領域配置

各ノードに２つずつ
メモリ領域配置• OpenMPと組み合わせて利用することが多い• OpenMP:一つのメモリ領域内で並列計算

各ノードに４つずつ
メモリ領域配置

4

MPIコースの内容• 第一部：• MPIプログラムの基本構成• コンパイルと実行• MPIの初期化と終了• 実習１• 第二部：• 計算とデータの分割• 実習２• 第三部：• 通信の記述• 実習３• 第四部：• より高度な MPIプログラミン
グに向けて
• 実習４

5

プログラム中の通信• インターネットプロトコル（TCP, UDP）の場合• 接続：• socket, bind, listen, connect, accept, ...

• ホストの識別：• ホスト名, IPアドレス, ポート番号, ...

• 転送

TCPによるサーバプログラム例• バイト単位• 基本的に一対一通信のみ

TCPによるクライアントプログラム例

6

MPI (Message Passing Interface)• 並列計算向けの通信関数群• C, C++, Fortranのプログラムから呼び出し• ほぼ全ての並列計算機で利用可能

直感的に並列プログラムを記述できるように, 通信を抽象化• 接続：• MPI_Init• ホストの識別：• MPI_Comm_rank()（識別番号の取得）• 転送• データ型単位• 新たにデータ型を定義可能• 一対一, 一対多, 多対多
MPIのプログラム例

7
MPIのプログラム例（C/C++）#include <stdio.h>#include "mpi.h"
int main(int argc, char *argv[]){
int myid, procs, i;double myval, val;
MPI_Status status;
FILE *fp;

MPIの準備MPI_Init(&argc, &argv);
自分のプロセス番号（＝ランク）を取得
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
実行に参加しているプロセス数を取得
MPI_Comm_size(MPI_COMM_WORLD, &procs);
ランク0か否か
if (myid == 0) {fp = fopen("test.dat", "r");最初のデータは, myvalに格納
fscanf(fp, "%lf", &myval);i = １～procs－１
for (i = 1; i < procs; i++){
次のデータを読み込み, valに格納
fscanf(fp, "%lf", &val);
MPI_Send(&val, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);}MPI_Sendにより, valの値をランク iに送信
fclose(fp);} elseMPI_Recv(&myval, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
ランク0以外のプロセスは, MPI_Recvでランク0から受信し, myvalに格納
printf("PROCS: %d, MYID: %d, MYVAL: %e¥n", procs, myid, myval);
MPI_Finalize();
myvalの値を表示
return 0;

MPIの終了処理}
7

8program ex1
implicit none
include "mpif.h"

MPIのプログラム例（Fortran）

integer :: myid, procs, ireal(8) :: myval, valinteger :: ierrinteger, dimension(MPI_STATUS_SIZE) :: statuscall MPI_Init(ierr)call MPI_Comm_rank(MPI_COMM_WORLD, myid, ierr)call MPI_Comm_size(MPI_COMM_WORLD, procs, ierr)

MPIの準備自分のプロセス番号（＝ランク）を取得

実行に参加しているプロセス数を取得
ランク0か否か

if (myid == 0) then
最初のデータは, myvalに格納
open(10, file="test.dat")read(10, *) myval
i = １～procs－１
do i = 1, procs-1
次のデータを読み込み, valに格納
read(10, *) val
call MPI_Send(val, 1, MPI_DOUBLE_PRECISION, i, 0, MPI_COMM_WORLD, ierr)end doMPI_Sendにより, valの値をランク iに送信
close(10)elsecall MPI_Recv(myval, 1, MPI_DOUBLE_PRECISION, 0, 0, MPI_COMM_WORLD, status, ierr)
end if
ランク0以外のプロセスは, MPI_Recvでランク0から受信し, myvalに格納
print *, "PROCS: ", procs, " MYID: ", myid, " MYVAL: ", myval
call MPI_Finalize(ierr)end program

myvalの値を表示

MPIの終了処理
8

9

プログラム例の実行の流れ• 複数の"プロセス"が，自分の番号（ランク）に応じて実行
rank 0
データreadmyval

rank 1

rank 2

ランク０から受信データreadval

ランク１にvalを送信
データread

ランク２にvalを送信myval表示

ランク０から受信データの到着待ち

myval
val

myval表示
データの到着待ち

myval
myval表示10

実行例• 各プロセスがそれぞれ勝手に表示するので, 表示の順番は毎回変わる可能性がある

PROCS: 4 MYID: 1 MYVAL: 20.0000000000000000
PROCS: 4 MYID: 2 MYVAL: 30.0000000000000000
PROCS: 4 MYID: 0 MYVAL: 10.0000000000000000
PROCS: 4 MYID: 3 MYVAL: 40.0000000000000000

rank 1
rank 2
rank 0rank 3

11

MPIインタフェースの特徴• C/C++, Fortranプログラムから呼び出す関数（サブルーチン）• 基本的に, 各プロセスが同じプログラムを実行する• ランク（＝プロセス番号）を使って，プロセス毎に違う仕事を実行• 他のプロセスの変数を直接読み書きすることはできない
Rank 0

Read file
Read file

myval
val

Rank 2
Receive

Send
Read file

Rank 1

val

Receivemyval

Print myval

Send
Print myval

myval
Print myval

11

12

MPIプログラムの基本構造
必須行#include <stdio.h>#include "mpi.h"

ヘッダファイル"mpi.h"int main(int argc, char *argv[]){...
MPI_Init(&argc, &argv);

MPIの準備...
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Comm_size(MPI_COMM_WORLD, &procs);...
MPI_Send(&val, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);...
MPI_Recv(&myval, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);

MPI関数を使用可能な範囲...
MPI_Finalize();
return 0;}

MPIの終了13

MPIライブラリ• MPI関数の実体は，MPIライブラリに格納されている• MPI用コンパイルコマンド（mpicc等）でMPIライブラリをプログラムに結合

mpicc
main(){
MPI_Init(...);...
MPI_Comm_rank(...);...
MPI_Send(...);...
}

source program

翻訳結合

MPI_Init
MPI_Comm_rank...

MPIライブラリ

実行ファイル

14

MPIプログラムのコンパイル• MPIライブラリが提供するコンパイルコマンドを利用
MPIライブラリ

コンパイルコマンド
C/C++

FortranIntel MPI

mpiiccmpiifort

Open MPI

mpicc

mpifortMVAPICH2

mpicc

mpifort富士通

mpifccmpifrt• 例）
mpiicc test.c -o test

15

MPIプログラムの実行• 通常, mpiexec コマンドで実行• -np オプションでプロセス数を指定• 例）
mpiexec -np 8 ./test• 計算機や MPIライブラリによって,コマンド名, オプション, 環境変数が違う
• マニュアル参照

16

ITOでの MPIプログラム実行• 「バッチ型」による利用

ログインノード
ITOポータル

フロントエンド• 基本的な流れ：

バックエンド１．ジョブスクリプト作成
２．ジョブ投入（システムの空き状況に応じて,順にジョブを実行）
３．完了を待って, 出力ファイル確認

計算結果

17

バッチシステムの仕組み• 処理してほしい内容を記述したファイルを投入• ジョブとして受付• 資源の空き状況に応じて順に処理される• 要求内容や空き状況によっては先を越されることも#!/bin/sh#PJM -L "vnode=4"#PJM -L "vnode-core=36"#PJM -L "rscunit=ito-b"#PJM -L "rscgrp=ito-g-16"#PJM -L "elapse=10:00"
mpiexec -np 16 ./a.out

ジョブスクリプトの例

リソースグループ（待ち行列）
ジョブの使用資源量（CPU,メモリ）で選択

ジョブ

小規模ジョブ用

中規模ジョブ用

大規模ジョブ用計算機

18

バッチ処理に用いるコマンド• バッチジョブの投入

pjsub• バッチジョブの状況

pjstat• バッチジョブのキャンセル

pjdel

19

pjsub バッチジョブの投入• コマンド$ pjsub オプション ジョブスクリプトファイル名• オプション：使用する資源等に関する指定• いつも同じ指定をするのであれば, ジョブスクリプトファイルの中に記述• pjsubコマンドでのオプション指定が優先• ジョブスクリプトファイル：依頼する処理内容• シェルスクリプトとして記述• 例）ジョブスクリプトファイル test.sh を投入$ pjsub test.sh[INFO] PJM 0000 pjsub Job 28246 submitted.

ジョブ ID

20

今回の実習で使うジョブスクリプト#!/bin/bash#PJM -L "rscunit=ito-a"#PJM -L "rscgrp=ito-a-lecture"#PJM -L "vnode=2" 使用ノード数 2#PJM -L "vnode-core=36" ノード当たり使用可能最大コア数 36#PJM -L "elapse=00:05:00" 実行時間 5分以内
#PJM -jmodule load intel/2017

Intel MPIのバージョン

export I_MPI_FABRICS=shm:ofa
export I_MPI_HYDRA_BOOTSTRAP=rsh
export I_MPI_HYDRA_BOOTSTRAP_EXEC=/bin/pjrsh
export I_MPI_HYDRA_HOST_FILE=${PJM_O_NODEINF}
export I_MPI_PERHOST=2Intel MPIの設定

ノードあたりのプロセス数

mpiexec.hydra -np 4 ./ex1

実行コマンド（プロセス数 <= vnode * I_MPI_PERHOSTとなるように, 設定）

21

MPIプログラムの作成• 以下の逐次プログラムを例に, 並列化の手順例を紹介
a = (double *)malloc(N*sizeof(double));newa = (double *)malloc(N*sizeof(double));allocate(a(0:n-1))
allocate(newa(0:n-1))

for (i = 1; i < N-1; i++)a[i] = 0.0;a[0] = 100.0;
a[N-1] = 10.0;a(1:n-2) = 0.0
a(0) = 100.0a(n-1) = 10.0

for (j = 0; j < REPEAT; j++){for (i = 1; i <= N-2; i++)newa[i] = (a[i-1]+a[i]+a[i+1])/3.0;
for (i = 1; i <= N-2; i++)a[i] = newa[i];printf("Step %2d: ", j);for (i = 0; i < N; i++)printf(" %6.2f", a[i]);printf("¥n");}[0]

do j = 1, repeat
do i = 1, n-2
newa(i) = (a(i-1) + a(i) + a(i+1)) / 3.0
end doa(1:n-2) = newa(1:n-2)
write(*,'(a5,i2,a1)',advance='no') "Step ",j,":"do i = 0, n-1
write(*,'(f7.2)',advance='no') a(i)end dowrite(*,*)
end do[19]

22

並列化の目標• 並列化前[0][19]• 並列化後rank 0[0][5]rank 1 [0][6]rank 2 [0][6]通信
rank 3 [0]

[5]

23

並列化の手順（例）• ステップ１：MPIの必須関数追加• ステップ２：計算の並列化とデータの分割• ループをプロセスで分担• （必要に応じて）ループの割当に合わせて配列を分割• ステップ３：通信関数追加• 隣接プロセスの計算結果取得• 全プロセスの計算結果結合24

ステップ１：MPIの必須関数追加• MPI_Init• MPIの初期化処理• MPI_Finalize• MPIの終了処理• ヘッダファイルの include• C/C++: mpi.h• Fortran: mpif.h25

MPI_Init
C/C++:
int MPI_Init(int *argc, char **argv);

Fortran:call MPI_Init(ierr)integer :: ierr

• MPI利用開始• プロセスの起動やプロセス間通信路の確立等• 他のMPI関数を呼ぶ前に, 必ずこの関数を呼ぶ• 引数• C/C++• argc, argv : main関数の2つの引数へのポインタ• 各プロセス起動時に実行ファイル名やオプションを共有するために参照• Fortran• ierr :エラー番号を返す整数変数• 全ての Fortran用 MPIルーチンに共通

26

MPI_Finalize
C/C++:
int MPI_Finalize();

Fortran:call MPI_Finalize(ierr)integer :: ierr

• MPI利用終了• このルーチン実行後はMPIルーチンを呼び出せない• プログラム終了前に全プロセスで必ずこのルーチンを実行

27

MPI_Init, _Finalizeの追加/* Initialize MPI */MPI_Init(&argc, &argv);! Initialize MPI
call MPI_Init(ierr);/* Finalize MPI */MPI_Finalize();! Finalize MPI
call MPI_Finalize(ierr);

実習時に追加するコード（C/C++）

実習時に追加するコード（Fortran）• 全プロセスが同じプログラムを実行
rank 0[0][19]

rank 1[0][19]

rank 2[0][19]

rank 3

[0]

[19]

28

実習１ 準備• ITOにログイン• /home/tmp/mpi/mpiex-2018 を自分のホームにコ

ピーし, mpiex-2018ディレクトリに移動
cp -r /home/tmp/mpi/mpiex-2018 .
cd mpiex-2018

29

実習1 MPI_Init, _Finalizeを追加• ex1.c もしくは ex1.f90 を編集し, 全プロセスで同じプ

ログラムを並列実行するプログラムに変更• ヘッダファイル, 変数宣言も追加/* Add header file */#include "mpi.h"! Add header fileinclude 'mpif.h'! Add necessary variable
integer :: ierr

• コンパイル後, ジョブ投入
module load intel/2017
mpiicc ex1.c -o ex1
cat ex1.sh
pjsub ex1.sh• 結果の確認ls
cat ex1.sh.oジョブ番号

module load intel/2017
mpiifort ex1.f90 -o ex1

30

ステップ２：
計算の並列化とデータの分割
• 分割前
[0]

[19]

• 分割後
rank 0[0][5]rank 1 [0][6]rank 2 [0]rank 3 [0]

[6]

[5]

31

計算の並列化• ループをプロセスに割り当て• 今回は「ブロック分割」を適用• ループを部分ループに分割してプロセスに割り当て

rank 0[0][19]

rank 1[0][19]

rank 2[0][19]

rank 3

[0]

[19]

32

他の割当方法• サイクリック
rank 0[0][19]

rank 1[0][19]

rank 2[0][19]

rank 3[0][19]• ブロックサイクリック
rank 0[0][19]

rank 1[0][19]

rank 2[0][19]

rank 3

[0]

[19]

33

（必要に応じて）データの分割• １台の計算ノードではメモリが不足する場合,

データ分割• 通常, 各プロセスが, 自分の担当範囲のデータを配置

rank 0[0][19]

rank 1[0][19]

rank 2[0][19]

rank 3

[0]

[19]

34

（必要に応じて）データ分割を調整• 今回のプログラムでは, 隣接の要素も参照するので,少し余裕を持って配列を配置しておく• あとで通信に利用する

rank 0[0][5]rank 1 [0][6]rank 2 [0]rank 3 [0]

[6]

[5]

35

Nがプロセス数で割り切れる場合のループとデータ分割• ランク myid, プロセス数 procs を利用
ランク

配列の大きさループ範囲

myid == 0

N/procs + 1

1 ... N/procs - 1

1 <= myid <= procs - 2

N/procs + 2

1 ... N/procs

myid == procs - 1

N/procs + 1

1 ... N/procs - 1

36

MPI_Comm_rank
C/C++:
int MPI_Comm_rank(MPI_Comm comm,int *rank);

Fortran:call MPI_Comm_rank(comm, rank, ierr)integer :: comm, rank, ierr

• そのプロセスのランク取得• 引数• comm : “コミュニケータ”• rank : ランクを格納する場所• コミュニケータ• プロセスのグループを表す識別子• 通常は，MPI_COMM_WORLD を指定• MPI_COMM_WORLD: 実行に参加する全プロセスによるグループ• プロセスを複数の小グループに分けて,それぞれ別の仕事をさせることも可能

37

MPI_Comm_size
C/C++:
int MPI_Comm_size(MPI_Comm comm,int *size);• プロセス数取得• 引数• comm : “コミュニケータ”• size : プロセス数を格納する場所

Fortran:call MPI_Comm_size(comm, size, ierr)integer :: comm, size, ierr

38

Nがプロセス数で割り切れる場合の配列サイズとループ終了値の計算例/* Get myid and procs */MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Comm_size(MPI_COMM_WORLD, &procs);! Get myid and procs
call MPI_Comm_rank(MPI_COMM_WORLD, myid, ierr)call MPI_Comm_size(MPI_COMM_WORLD, procs, ierr)/* Divide N into procs */divN = N / procs;
modN = N % procs;
if (modN != 0)
printf("ERROR: modN is not 0 ¥n");! Divide N into procs
divN = floor(real(N/procs))
modN = mod(N, procs)
if (modN /= 0) thenwrite(*, *) "Error modN is not 0"
end ifif ((myid == 0) || (myid == procs-1)){
sizea = divN + 1;
loopend = divN - 1;
} else {sizea = divN + 2;
loopend = divN;}

if ((myid == 0) .or. (myid == procs-1)) thensizea = divN + 1
loopend = divN - 1
else
sizea = divN + 2
loopend = divN
endif

39

Nがプロセス数で割り切れない場合• 並列化前（N = 18）
[0]

[17]

• 並列化後（プロセス数 ４）
rank 0[5][0]rank 1 [0][6]rank 2 [0]rank 3 [0]

[5]

[4]

40

Nがプロセス数で割り切れない場合のループとデータ分割の例• 余った分（= N % procs）を前半のプロセスに分配• divN = N / procs （切り捨て）
ランク

配列の大きさループ範囲

myid == 0

divN + 2

1 ... divN

1 <= myid <= (N % procs) - 1

divN + 3

1 ... divN +1(N % procs) <= myid <= procs - 2

divN + 2

1 ... divN

myid == procs - 1

divN + 1

1 ... divN - 1

41

Nがプロセス数で割り切れない場合の配列サイズとループ終了値の計算例/* Divide N into procs */divN = N / procs;
modN = N % procs;! Divide N into procs
divN = floor(real(N/procs))
modN = mod(N, procs)

if ((myid == 0) || (myid == procs-1) {sizea = divN + 1;
loopend = divN - 1;
} else {sizea = divN + 2;
loopend = divN;}

if ((myid == 0) .or. (myid == procs-1))
thensizea = divN + 1
loopend = divN - 1
else
sizea = divN + 2
loopend = divN
endif

if (myid < modN) {sizea++;
loopend++;}if (myid < modN) thensizea = sizea + 1
loopend = divN + 1
end if42

配列の初期化例• 配列サイズ変更• 初期化ループの範囲変更/* Initialize local a */a = (double *)malloc(sizea*sizeof(double));newa = (double *)malloc(sizea*sizeof(double));for (i = 1; i <= loopend; i++)a[i] = 0.0;
if (myid == 0)a[0] = 100.0;
if (myid == procs - 1)
a[sizea - 1] = 10.0;! Initialize local a
allocate(a(0:sizea-1))
allocate(newa(0:sizea-1))a(1:loopend) = 0.0
if (myid == 0) thena(0) = 100.0
end ifif (myid == procs - 1) thena(sizea-1) = 10.0
end if

43

ループの並列化例• ループの範囲修正• 表示部では, とりあえずランク0の配列のみ表示することにする/* Calculate local area */for (i = 1; i <= loopend; i++)newa[i] = (a[i-1]+a[i]+a[i+1])/3.0;! Calculate local area
do i = 1, loopend
newa(i) = (a(i-1) + a(i) + a(i+1)) / 3.0
end do/* Update data */for (i = 1; i <= loopend; i++)a[i] = newa[i];! Update dataa(1:loopend) = newa(1:loopend)/* Print local area */if (myid == 0) {printf("Step %2d: ", j);for (i = 0; i < sizea; i++)printf(" %6.2f", a[i]);printf("¥n");}! Print local area
if (myid == 0) thenwrite(*,'(a5,i2,a1)',advance='no') "Step ",j,":"do i = 0, sizea-1
write(*,'(f7.2)',advance='no') a(i)end dowrite(*,*)
end if

44

実習2 計算の並列化とデータ分割• ex2.c もしくは ex2.f90 を編集し, 計算を並列化してデータを分割するプログラムに変更• 必要な変数宣言を追加/* Add necessary variables */double *a, *newa;
int i,• j;
int myid, procs, modN, divN, sizea, loopend;! Add necessary variable
integer :: i, j, ierrinteger :: myid, procs, modN, divN, sizea, loopend• ex2-answer.c, ex2-answer.f90に回答例があります• コンパイル後, ジョブ投入
mpiicc ex2.c -o ex2
pjsub ex2.sh• 結果の確認
pjstatls
cat ex2.sh.oジョブ番号

mpiifort ex2.f90 -o ex2
pjsub ex2.sh45

ステップ３：
通信の追加• 他のプロセスの計算結果を参照• 今回のプログラムで使用する通信：• 一対一通信• 隣のプロセスと,境界部分の値を交換• 集団通信（グループ通信）• 全プロセスの計算結果の連結に利用• 配列表示用46

隣のプロセスとの値交換• 左のプロセス（myid-1）との交換（ランク 0以外）：• a[1] を送信• a[0] に受信• 右のプロセス（myid+1）との交換（ランクprocs-1 以外）：
• a[loopend] を送信• a[loopend+1] に受信

47

一対一通信関数• MPI_Send, MPI_Recv• ブロッキング送信, 受信• MPI_Isend, MPI_Irecv• 非ブロッキング送信, 受信• MPI_Wait, MPI_Waitall• 非ブロッキング送信, 受信の完了待ち

48

ブロッキング通信
vs 非ブロッキング通信• ブロッキング通信：通信の完了を待って, 次の命令へ• 非ブロッキング通信：通信の完了を待たずに, 次の命令へ• 通信完了は, 別途, MPI_Wait関数等で待つ
Blocking

Non-Blocking
MPI_Recv

Wait for thearrival of data

MPI_Irecv

data
next instructions

Proceed to the next
instruction without
waiting for the data

next instructions
MPI_Wait

data

49

非ブロッキング通信の利点• デッドロックの回避• お互いに相手のデータの送信待ち, という状態を回避• 通信時間の隠蔽• データが転送されている間に他の計算を進める50

非ブロッキング通信によるデッドロックの回避• デッドロック：

何らかの理由でプログラムを進行できなくなった状態• MPIプログラムにおけるデッドロックの例：• お互い, 相手からデータが送信されないと自分の送信を始めない• 非ブロッキング通信で回避：• 受信開始後, 完了を待たずに送信開始

if (myid == 0){
MPI_Recv from rank 1
MPI_Send to rank 1}
if (myid == 1){
MPI_Recv from rank 0
MPI_Send to rank 1}
if (myid == 0){
MPI_Irecv from rank 1
MPI_Send to rank 1
MPI_Wait}if (myid == 1){
MPI_Irecv from rank 0
MPI_Send to rank 0
MPI_Wait}

51

非ブロッキング通信による通信時間の隠蔽• 非ブロッキング通信開始後, その通信と無関係な計算実行⇒ 見かけ上, 通信時間が ０に近づく• 今回のプログラム• 通信中に, 隣のプロセスの値を必要としない範囲の計算が可能
rank 0[0][5]rank 1 [0][6]rank 2 [0]rank 3 [0]

[6]

[5]

52

非ブロッキング通信の注意点• 通信開始後, 完了待ち（Wait）までは,通信対象領域を読み書きしない• MPI_Isend• 完了前に送信対象領域を書き換えると,書き換え前と後のどちらのデータが送信されるか不明

MPI_Isend(A)A = ...
MPI_Wait()• MPI_Irecv• 完了前に受信対象領域を参照すると,データ到着前と後のどちらのデータを参照するか不明

MPI_Irecv(A)... = A
MPI_Wait()

53

通信を隠蔽するプログラムの流れ• 計算を, 内部と境界部に分ける/* Exchange values */if (myid > 0) { /* left */MPI_Isend(左にa[1] を送信開始);
MPI_Irecv(左からa[0] に受信開始);}
if (myid < procs-1) {/* right */MPI_Isend(右にa[loopend]を送信開始);
MPI_Irecv(右からa[loopend+1]に受信開始);}! Exchange values
if (myid > 0) thencall MPI_Isend(左にa(1)を送信開始)
call MPI_Irecv(左からa(0)に受信開始)
end ifif (myid < procs - 1) thencall MPI_Isend(右にa(loopend)を送信開始)
call MPI_Irecv(右からa(loopend+1を受信開始)
end if/* Modify loop to calculate internal area only */for (i = 2; i <= loopend-1; i++)newa[i] = (a[i-1]+a[i]+a[i+1])/3.0;! Modify loop to calculate internal area only
do i = 2, loopend-1
newa(i) = (a(i-1) + a(i) + a(i+1)) / 3.0
end do/* Wait for all non-blocking communications */MPI_Waitall(全ての非ブロッキング通信);/* Calculate edges */newa[1] =(a[0]+a[1]+a[2])/3.0;
newa[loopend] = (a[loopend-1]+a[loopend]+a[loopend+1])/3.0;! Wait for all non-blocking communications
call MPI_Waitall(全ての非ブロッキング通信)! Calculate edges
newa(1) = (a(0) + a(1) + a(2)) / 3.0
newa(loopend) = (a(loopend-1) + a(loopend) + &a(loopend+1)) / 3.0

54

MPI_Send
C/C++:
int MPI_Send(void *b, int c, MPI_Datatype d,
int dest, int t, MPI_Comm comm);• ブロッキング送信• 引数• b : 送信データの先頭• c : 要素数
• d : データ型
• dest : 送信先ランク• t : タグ• メッセージにつける番号• 不規則な通信の記述に利用• 通常は 0 で可• comm : コミュニケータ
• r : リクエスト情報の格納場所

Fortran:call MPI_Send(b, c, d, dest, t, comm, ierr)<type> :: binteger :: c, d, dest, t, comm, ierr

55

MPIの主なデータ型
データ型

MPIのデータ型（C/C++） MPIのデータ型（Fortran）

整数

MPI_INTMPI_INTEGER

単精度実数

MPI_FLOATMPI_REAL

倍精度実数

MPI_DOUBLEMPI_DOUBLE_PRECISION

単精度複素数MPI_COMPLEX

倍精度複素数MPI_DOUBLE_COMPLEX文字

MPI_CHARMPI_CHARACTER

56

MPI_Recv
C/C++:
int MPI_Recv(void *b, int c, MPI_Datatype d,
int dest, int t, MPI_Comm comm,
MPI_Status *s);

Fortran:call MPI_Isend(b, c, d, dest, t, comm, s, ierr)<type> :: binteger :: c, d, dest, t, comm, ierrinteger,dimension(MPI_STATUS_SIZE) :: s

• ブロッキング受信• 引数• b : 受信データの格納場所の先頭• c : 要素数
• d : データ型
• dest : 受信元ランク• t : タグ• comm : コミュニケータ
• s : 受信したデータの情報の格納場所• 送信元ランク, タグの値, 等• 不規則な通信に利用• 通常は MPI_STATUS_IGNORE を指定

57

MPI_Isend
C/C++:
int MPI_Isend(void *b, int c, MPI_Datatype d,
int dest, int t, MPI_Comm comm,
MPI_Request *r);

• 非ブロッキング送信• 引数• b : 送信データの先頭• c : 要素数
• d : データ型
• dest : 送信先ランク• t : タグ• comm : コミュニケータ
• r : リクエスト情報の格納場所

Fortran:call MPI_Isend(b, c, d, dest, t, comm, r, ierr)<type> :: binteger :: c, d, dest, t, comm, r, ierr

58

MPI_Irecv
C/C++:
int MPI_Irecv(void *b, int c, MPI_Datatype d,
int dest, int t, MPI_Comm comm,
MPI_Request *r);

Fortran:call MPI_Irecv(b, c, d, dest, t, comm, r, ierr)<type> :: binteger :: c, d, dest, t, comm, r, ierr

• 非ブロッキング受信• 引数• b : 受信データの格納場所の先頭• c : 要素数
• d : データ型
• dest : 受信元ランク• t : タグ• comm : コミュニケータ
• r : リクエスト情報の格納場所

59

MPI_Wait
C/C++:
int MPI_Wait (MPI_Request *r, MPI_Status *s);

Fortran:call MPI_Wait(r, s)
integer :: rinteger,dimension(MPI_STATUS_SIZE) :: s

• 非ブロッキング通信の完了待ち• 引数• r : リクエスト情報の格納場所• s : 受信したデータの情報の格納場所• MPI_STATUS_IGNORE指定可

60

リクエスト情報• 非ブロッキング通信の完了待ちに必要な情報• 各非ブロッキング通信発行毎に, 別の場所に保存• 完了待ち後は, 再利用可能
for () {MPI_Isend(..., &r1);
MPI_Irecv(..., &r2);...
MPI_Wait(&r1, MPI_IGNORE_STATUS);
MPI_Wait(&r2, MPI_IGNORE_STATUS);}
非ブロッキング通信の完了待ち（C/C++）docall MPI_Isend(..., r1, ...);
call MPI_Irecv(..., r2, ...);...
call MPI_Wait(r1, MPI_IGNORE_STATUS);
call MPI_Wait(r2, MPI_IGNORE_STATUS);
end do非ブロッキング通信の完了待ち（Fortran）

61

MPI_Waitall
C/C++:
int MPI_Waitall (int c, MPI_Request *r,
MPI_Status *s);

Fortran:call MPI_Irecv(c, r, s)
integer :: cinteger,dimension(:) :: rinteger,dimension(MPI_STATUS_SIZE,:) :: s

• 複数の非ブロッキング通信の完了待ち• 引数• c : 待つリクエストの数• r : リクエスト情報の配列• s : 受信したデータの情報の格納場所の配列• MPI_STATUS_IGNORE指定可

62

非ブロッキング通信によるデータ交換と並列計算の例（C/C++）/* Allocate an array of requests */reqs = (MPI_Request *)malloc(4*sizeof(MPI_REQUEST));/* Exchange values */nreqs = 0
if (myid > 0) { /* left */MPI_Isend(&(a[1]), 1, MPI_DOUBLE, myid-1, 0, MPI_COMM_WORLD, &(reqs[nreqs]));nreqs++;
MPI_Irecv(&(a[0]), 1, MPI_DOUBLE, myid-1, 0, MPI_COMM_WORLD, &(reqs[nreqs]));nreqs++;}
if (myid < procs-1) {/* right */MPI_Isend(&(a[loopend]), 1, MPI_DOUBLE, myid+1, 0, MPI_COMM_WORLD, &(reqs[nreqs]));nreqs++;
MPI_Irecv(&(a[loopend+1]), 1, MPI_DOUBLE, myid+1, 0, MPI_COMM_WORLD, &(reqs[nreqs]));nreqs++;}/* Modify loop to calculate internal area only */for (i = 2; i <= loopend - 1; i++)newa[i] = (a[i-1] + a[i] + a[i+1])/3.0;/* Wait for all non-blocking communications */MPI_Waitall(nreqs, reqs, MPI_STATUS_IGNORE);/* Calculate edges */newa[1] = (a[0] + a[1] + a[2])/3.0;
newa[loopend] = (a[loopend-1] + a[loopend] + a[loopend+1])/3.0;

63

非ブロッキング通信によるデータ交換と並列計算の例（Fortran）! Allocate an array of requests
allocate(reqs(0:3))! Exchange values
nreqs = 0
if (myid > 0) thencall MPI_Isend(a(1), 1, MPI_DOUBLE_PRECISION, myid - 1, 0, MPI_COMM_WORLD, reqs(nreqs), ierr)nreqs = nreqs + 1
call MPI_Irecv(a(0), 1, MPI_DOUBLE_PRECISION, myid - 1, 0, MPI_COMM_WORLD, reqs(nreqs), ierr)nreqs = nreqs + 1
end ifif (myid < procs - 1) thencall MPI_Isend(a(loopend), 1, MPI_DOUBLE_PRECISION, myid + 1, 0, MPI_COMM_WORLD, reqs(nreqs), ierr)nreqs = nreqs + 1
call MPI_Irecv(a(loopend+1), 1, MPI_DOUBLE_PRECISION, myid + 1, 0, MPI_COMM_WORLD, reqs(nreqs), ierr)nreqs = nreqs + 1
end if! Modify loop to calculate internal area only
do i = 2, loopend-1
newa(i) = (a(i-1) + a(i) + a(i+1)) / 3.0
end do! Wait for all non-blocking communications
call MPI_Waitall(nreqs, reqs, MPI_STATUS_IGNORE, ierr)! Calculate edges
newa(1) = (a(0) + a(1) + a(2)) / 3.0
newa(loopend) = (a(loopend-1) + a(loopend) + a(loopend+1)) / 3.0

64

全プロセスの計算結果の収集• 今回のプログラムでは,計算結果を一つのプロセスにまとめて, 表示したい• 各プロセスで表示するとバラバラになるため• 集団通信関数の一つ, MPI_Gather関数を利用
[0]

[19]

rank 0[5][0]rank 1 [0][6]rank 2 [0][6]rank 3 [0]

[5]

65

MPI_Gather
C/C++:
int MPI_Gather(void *b1, int c1, MPI_Datatype d1,
void *b2, int c2, MPI_Datatype d2,
int root, MPI_Comm comm);

Fortran:call MPI_Gather(b1, c1, d1, b2, c2, d2, root,
comm, ierr)<type>,dimension(:) :: b1, b2
integer :: c1, d1, c2, d2, root, comm, ierr

• 全プロセスのデータを一つの配列にプロセス順に連結して格納• 引数• b1 : 送信データの先頭• c1 : 送信データの要素数• d1 : 送信データのデータ型• b2 : 受信データの先頭• c2 : 受信データの要素数• d2 : 受信データのデータ型• root : 収集したデータを格納するランク• comm : コミュニケータ
• 注意）b2, c2, d2 は, rootのみで意味を持つ

66

データを集めて表示する例（C/C++）• rank0のみMPI_Gatherの先頭が違う/* Allocate an array for MPI_Gather */if (myid == 0)
worka = (double *)malloc(N*sizeof(double));/* Set start address for MPI_Gather */if (myid == 0)starta = 0;
elsestarta = 1;/* Gather data to rank 0 and print */MPI_Gather(&(a[starta]), divN, MPI_DOUBLE, worka, divN, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if (myid == 0) {printf("Step %2d: ", j);for (i = 0; i < N; i++)printf(" %6.2f", worka[i]);printf("¥n");}/* Free the array for MPI_Gather */if (myid == 0)free(worka);

67

データを集めて表示する例（Fortran）• rank0のみMPI_Gatherの先頭が違う! Allocate an array for MPI_Gather
if (myid == 0) thenallocate(worka(0:n-1))
end if! Set start address for MPI_Gather
if (myid == 0) thenstarta = 0
else
starta = 1
end if!Gather data into rank 0 and print
call MPI_Gather(a(starta), divN, MPI_DOUBLE_PRECISION, worka, divN, MPI_DOUBLE_PRECISION, 0,&
MPI_COMM_WORLD, ierr)
if (myid == 0) thenwrite(*,'(a5,i2,a1)',advance='no') "Step ",j,":"do i = 0, n-1
write(*,'(f7.2)',advance='no') worka(i)end dowrite(*,*)
end if! Free the array for MPI_Gather
if (myid == 0) thendeallocate(worka)
end if68

終了時の注意• 動的配列の解放は全プロセスの通信が完了後• 未完了の通信による領域違反を防ぐ• 同期関数を利用
free(a);free(newa);

deallocate(a)deallocate(newa)

MPI_Finalize();call MPI_Finalize(ierr)

動的配列の解放部分（C/C++）

動的配列の解放部分（Fortran）

69

MPI_Barrier
C/C++:
int MPI_Barrier(MPI_Comm comm);

Fortran:call MPI_Barrier(comm, ierr)integer :: comm, ierr

• 全プロセスが到着するまで待つ（同期）• 引数• comm : コミュニケータ

/* Wait for other process before free */MPI_Barrier(MPI_COMM_WORLD);! Wait for other process before freecall MPI_Barrier(MPI_COMM_WORLD, ierr)70

実習3 通信• ex3.c もしくは ex3.f90 を編集し, 隣のプロセスとのデータ

交換, および計算結果の収集のための通信を行うプログラムに変更• 変数定義の追加/* Add necessary variables */double *a, *newa, *worka;
int i, j, nreqs, starta;
int myid, procs, modN, divN, sizea, loopend;
MPI_Request *reqs;! Add necessary variables
real(8), dimension(:), allocatable :: a, newa, worka
integer :: i, j, ierr, nreqs, starta
integer :: myid, procs, modN, divN, sizea, loopend
integer, dimension(:), allocatable :: reqs

• ex3-answer.c, ex3-answer.f90 に回答例があります• コンパイル後, ジョブ投入
mpiicc ex3.c -o ex3
pjsub ex3.sh• 結果の確認
pjstatls
cat ex3.sh.oジョブ番号

mpiifort ex3.f90 -o ex3
pjsub ex3.sh

71

Nがプロセス数で割り切れない場合の全プロセスの計算結果の収集• 送信データの要素数が, プロセスごとに異なる
[17]

[0]

rank 0[5][0]rank 1 [0][6]rank 2 [0][5]rank 3 [0]

[4]

72

MPI_Gatherv
C/C++:
int MPI_Gatherv(
void *b1, int c1, MPI_Datatype d1,
void *b2, const int x[], const int y[],MPI_Datatype d2,
int root, MPI_Comm comm);

Fortran:call MPI_Gatherv(b1, c1, d1, b2, x, y, d2, root,
comm, ierr)<type>,dimension(:) :: b1, b2
integer :: c1, d1, d2, root, comm, ierrinteger,dimension(:) :: x, y

• プロセスごとにサイズが違うデータの連結• 引数• b1 : 送信データの先頭• c1 : 送信データの要素数• d1 : 送信データのデータ型• b2 : 受信データの先頭• x : 各ランクからの送信データ要素数の配列• y : 各ランクからの送信データを格納する位置の配列• d2 : 受信データのデータ型• root : 収集したデータを格納するランク• comm : コミュニケータ

73

Nがプロセス数で割り切れない場合にデータを集めて表示する例（C/C++）/* Allocate an array for MPI_Gather */if (myid == 0)
worka = (double *)malloc(N * sizeof(double));/* Prepare tables and gatherlen for MPI_Gatherv */elems = (int *)malloc(procs*sizeof(int));
displs = (int *)malloc(procs*sizeof(int));
for (i = 0; i < modN; i++)elems[i] = divN + 1;
for (i = modN; i < procs; i++)elems[i] = divN;
displs[0] = 0;
for (i = 1; i < procs; i++)displs[i] = displs[i-1] + elems[i-1];
if (myid < modN)
gatherlen = divN + 1;
else
gatherlen = divN;/* Gather data into rank 0 and print */if (modN == 0)
MPI_Gather(&(a[starta]), divN, MPI_DOUBLE, worka, divN, MPI_DOUBLE, 0, MPI_COMM_WORLD);
else
MPI_Gatherv(&(a[starta]), gatherlen, MPI_DOUBLE, worka, elems, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if (myid == 0) {printf("Step %2d: ", j);for (i = 0; i < N; i++)printf(" %6.2f", worka[i]);printf("¥n");}

74

Nがプロセス数で割り切れない場合にデータを集めて表示する例（Fortran）! Allocate an array for MPI_Gather
if (myid == 0) thenallocate(worka(0:n-1))
end if! Prepare tables and gatherlen for MPI_Gatherv
allocate(elems(0:procs-1))
allocate(displs(0:procs-1))
do i = 0, modN-1
elems(i) = divN + 1
end dodo i = modN, procs - 1
elems(i) = divN
end dodispls(0) = 0
do i = 1, procs - 1
displs(i) = displs(i-1) + elems(i-1)
end doif (myid < modN) thengatherlen = divN + 1
else
gatherlen = divN
end if!Gather data into rank 0 and print
if (modN == 0) thencall MPI_Gather(a(starta), divN, MPI_DOUBLE_PRECISION, worka, divN, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)else
call MPI_Gatherv(a(starta), divN, MPI_DOUBLE_PRECISION, worka, elems, displs, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
end ifif (myid == 0) thenwrite(*, '(a5,i2,a1)', advance='no') "Step ", j, ":"do i = 0, n-1
write(*,'(f7.2)',advance='no') worka(i)end dowrite(*,*)
end if

75

その他の集団通信の例• MPI_Bcast• 全プロセスにコピー

Rank 0• MPI_Scatter• 全プロセスに分散

Rank 0Rank 13 1 8 23 1 8 2

Rank 23 1 8 2

Rank 1

7

5

Rank 297 5 9• MPI_Reduce• 全プロセスの値を集約

Rank 0Rank 1

1 2 3
12

15

4 5 6

18

• MPI_Allgather, MPI_Allreduce• MPI_Gather, MPI_Reduceの結果を全プロセスにコピー

Rank 27 8 9

76

集団通信の利用に当たって• 同じ関数を全プロセスが実行するよう, 記述する
• 例えば MPI_Bcastは，root rankが送信するデータを他のランクが受信する• このように集団通信は全プロセスが参加する• 送信データと受信データの場所を別々に指定するタイプの集団通信では, 送信データの範囲と受信データの範囲が重ならないように指定する• MPI_Gather, MPI_Allgather, MPI_Gatherv, MPI_Allgatherv,

MPI_Recude, MPI_Allreduce, MPI_Alltoall, MPI_Alltoallv, etc.

77

並列処理に対する期待と現実• プログラマ:「ＣＰＵを 4台使うんだから，並列化で 4倍速くなって欲しい」• 計算機製作者:「ＣＰＵ 4台で 3倍くらい速くなれば十分だろう」• Why?
• アムダールの法則• 負荷のバランス• 通信のコスト

78

アムダールの法則• プログラム中の高速化した部分しか高速化されない• 並列化にあてはめて考えると：並列化による性能向上率の理論的な限界＝1/((1-P)+P/N)• P: プログラム中の並列化対象部分が全処理時間に占める割合• N: プロセス数• Example) N=4 で 3.5倍以上高速化するためには 95%以上の部分の並列化が必要

79

負荷のバランス• 並列プログラムの処理時間は「最も遅いプロセスの処理時間」
である
Execution timeof this program

Execution timeof this program

Rank 0Rank 0Rank 1

Rank 1

Rank 2

Rank 2

Rank 3

Rank 380

通信時間• 並列化前は不要だった時間＝

並列化によるオーバーヘッド

Rank 0Rank 1
Rank 2
Rank 381

並列化が難しいプログラム• 実行の順序によって結果が変わる• Example:
前に計算した値を参照して計算

for (i = 1; i < 12; i++)a[i] = a[i] + a[i-1];

Rank 0 Rank 1 Rank 2 Rank 3
a

1

2

3

4

5 6

a

1

3

6

?

?

7

8

? ? ?

9 10 11 12
?

?

?

?

81

82

MPIプログラムの時間計測• MPI_Wtime• 現在時間（秒）を実数で返す関数 Returns the current time inseconds.
• Example)

Measuretime here...
double t1, t2;...
t1 = MPI_Wtime();処理
t2 = MPI_Wtime();printf("Elapsed time: %e sec.¥n", t2 – t1);

83

並列プログラムにおける時間計測の問題• プロセス毎に違う時間を測定：

どの時間が本当の所要時間か?

Rank 0t1 = MPI_Wtime();Rank 1

Read

Measuretime here

t1 = MPI_Wtime();
Read
Send
Read

Rank 2
t1 = MPI_Wtime();Receive
Receive
t1 = MPI_Wtime();Send
t1 = MPI_Wtime();t1 = MPI_Wtime();84

集団通信 MPI_Barrierを使った解決策• 時間計測前にMPI_Barrierで同期

Rank 0

MPI_Barrier

t1 = MPI_Wtime();

Measuretime here

Rank 1

MPI_Barrier
Receive

Rank 2

MPI_Barrier
Receive

Read
Read
Send
Read

MPI_Barrier

Send

MPI_Barrier
t1 = MPI_Wtime();MPI_Barrier

85

実習4 性能計測• 配列サイズを大きくして時間計測関数を追加してある

ex4.c もしくは ex4.f90をコンパイル後, ジョブ投入
mpiicc ex4.c -o ex4
pjsub ex4.sh• 結果の確認
pjstatls
cat ex4.sh.oジョブ番号

mpiifort ex4.f90 -o ex4
pjsub ex4.sh

86

MPI + OpenMP• MPIの各プロセスを, OpenMPのスレッドで並列化可能• 並列リージョン内での MPI関数利用は要注意/* Exchange values */if (myid > 0) { /* left */MPI_Isend(左にlocala[1] を送信開始);
MPI_Irecv(左からlocala[0] に受信開始);}
if (myid < procs-1) {/* right */MPI_Isend(右にlocala[loopend]を送信開始);
MPI_Irecv(右からlocala[loopend+1]に受信開始);}/* Modify loop to calculate internal area only */#pragma omp parallel forfor (i = 2; i <= loopend-1; i++)localnewa[i] =(locala[i-1]+locala[i]+locala[i+1])/3.0;! Exchange values
if (myid > 0) thencall MPI_Isend(左にa(1)を送信開始)
call MPI_Irecv(左からa(0)に受信開始)
end ifif (myid < procs - 1) thencall MPI_Isend(右にa(loopend)を送信開始)
call MPI_Irecv(右からa(loopend+1を受信開始)
end if/* Wait for all non-blockings */MPI_Waitall(全ての非ブロッキング通信);! Wait for all non-blocking communications
call MPI_Waitall(全ての非ブロッキング通信)/* Calculate edges */localnewa[1] =(locala[0]+locala[1]+locala[2])/3.0;
localnewa[loopend] = (locala[loopend-1]+locala[loopend]+locala[loopend+1])/3.0;! Modify loop to calculate internal area only!$omp parallel dodo i = 2, loopend-1
newa(i) = (a(i-1) + a(i) + a(i+1)) / 3.0
end do! Calculate edges
newa(1) = (a(0) + a(1) + a(2)) / 3.0
newa(loopend) = (a(loopend-1) + a(loopend) + &a(loopend+1)) / 3.0

87

MPI + OpenMP の実行• プロセス数とスレッド数を調整

2プロセス x 4スレッド

4プロセス x 2スレッド• 実行例）４プロセス x ２スレッド

8プロセス（OpenMP無し）

export OMP_NUM_THREADS=2
mpiexec -np 4 ./test• 並列リージョンで MPI関数を呼ぶプログラム• 使用しているMPIライブラリが対応しているか, 管理者に確認• MPI_Initを MPI_Init_threadに変更
MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &p);
if (p != MPI_THREAD_MULTIPLE)
printf("Error¥n");88

MPIの参考資料• 片桐 孝洋「スパコンプログラミング入門: 並列処理とMPIの学習」, 東京大学出版会, 2013

• P. Pacheco 「MPI並列プログラミング」, 培風館, 2001
• M. J. Quinn "Parallel Programming in C with MPI and

OpenMP", McGraw-Hill, 2003