内容に関する質問はkatagiri@cc.u-tokyo.ac.jp
まで

第４回 Hybrid並列化技法(MPIとOpenMPの応用)
東京大学情報基盤センター

1

２０１３年度

計算科学技術特論A

片桐孝洋

講義日程と内容について
 ２０１３年度

計算科学技術特論Ａ（１学期：木曜3限 ）第１回：プログラム高速化の基礎、2013年4月11日
 イントロダクション、ループアンローリング、キャッシュブロック化、
数値計算ライブラリの利用、その他
 第２回：MPIの基礎、2013年4月18日
 並列処理の基礎、MPIインターフェース、MPI通信の種類、その他
 第３回：OpenMPの基礎、2013年4月25日
 OpenMPの基礎、利用方法、その他
 第４回：Hybrid並列化技法(MPIとOpenMPの応用)、2013年5月9日
 背景、Hybrid並列化の適用事例、利用上の注意、その他
 第５回：プログラム高速化の応用、2013年5月16日
 プログラムの性能ボトルネック に関する考えかた（I/O、単体性能(演算機ネック、メモリネック)、並列性能(バランス))、性能プロファイル、
その他


2

２０１３年度

計算科学技術特論A

実際の並列計算機構成例

3２０１３年度

計算科学技術特論A

東京大学情報基盤センタースパコン
T2Kオープンスパコン（東大版）（HA8000クラスタシステム）
Total Peak performance: 140 TFLOPS
Total number of nodes: 952
Total memory: 32000 GB
Peak performance per node : 147.2 GFLOPSMain memory per node
: 32 GB, 128 GB
Disk capacity: 1 PB
AMD Quad Core Opteron （2.3GHz）
製品名：HITACHI HA8000-tc/RS425

4２０１３年度

計算科学技術特論A
4

各ＣＰＵの内部構成

T2K東大 ノード構成（タイプＡ群）

L3（4コア共有、2MB）

ソケット、ノードとは
8GB
ノード内
合計
メモリ量：
8GB×４
＝３２GB

Memory #0

8GB

512MB

L2

L2

L2

L2

64KB

L1

L1

L1

L1

Core
#0

Core
#1

Core
#2

Core
#3

Memory #1

コア
AMD
Quad Core
Opteron#1

AMD QuadCore Opteron
#0

ソケット

South
BridgeMyrinetMyrinetソケット
総合BW:42.6 GB/秒

AMD QuadCore Opteron
#2

ソケット
8GB
5

Memory #2

AMD QuadCore Opteron
#3

ソケット
8GB

Memory #3

２０１３年度

計算科学技術特論ASouthBridgeMyrinetMyrinetGbE
RAID

T2K（東大）での全体メモリ構成図
コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

…

コア
１２

コア
１３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

…

コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

…

Ｌ３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３
メインメモリ

コア
１２

コア
１３

コア
１４

コア
１５

コア
０

コア
１

コア
２

コア
３

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３

… …

…

Ｌ３

メインメモリ

コア
１２

コア
１３

コア
１４

コア
１５

Ｌ１

Ｌ１

Ｌ１

Ｌ１

Ｌ２

Ｌ２

Ｌ２

Ｌ２

Ｌ３
メインメモリ

メモリが多段に階層化（L1、L2、L3、分散メモリ）
6

コア
１３

Ｌ３

メインメモリ

コア
０

…

コア
１２

２０１３年度

計算科学技術特論A

高速ネットワーク（５Ｇバイト／秒
×双方向）（タイプＡ群）

ノード構成（T2K東大、タイプＡ群）

各ＣＰＵの内部構成
L3

ｃｃNUMA構成とはMemory #0

Memory #1

AMD QuadCore Opteron#0

AMD
Quad Core
Opteron
#1

アクセス
速い

共有メモリ
でアクセス
時間が
不均一
（ccNUMA）

L2

L2

L2

L1

L1

L1

L1

Core
#0

Core
#1

Core
#2

Core
#3

South
BridgeMyrinetMyrinet

アクセス遅い
AMD QuadCore Opteron#2

Memory #2
7

L2

AMD QuadCore Opteron#3

Memory #3
２０１３年度

計算科学技術特論ASouthBridgeMyrinetMyrinetGbE
RAID

東京大学情報基盤センタースパコン
FX10スーパーコンピュータシステム
Total Peak performance: 1.13 PFLOPS
Total number of nodes: 4,800
Total memory: 150TB
Peak performance per node : 236.5 GFLOPSMain memory per node
: 32 GBDisk capacity: 2.1 PBSPARC64 IXfx （1.848GHz）
製品名：Fujitsu PRIMEHPC FX10

2012年4月運用開始
8

２０１３年度

計算科学技術特論A8

FX10計算ノードの構成

TOFU
Network

１ソケットのみ
各ＣＰＵの内部構成

20GB/秒

Core
#0

Core
#1

Core
#2

Core
#3

…

Core
#12

Core
#13

Core
#14

Core
#15

L1

L1

L1

L1

: L1データキャッシュ32KB

L1

L1

L1

L1

ICC

L2 (16コアで共有、12MB)

85GB/秒=(8Byte×1333MHz
×8 channel)

MemoryMemoryMemoryMemory4GB ×2枚

4GB ×2枚

4GB ×2枚

4GB ×2枚

ノード内合計メモリ量：8GB×4＝32GB
9

２０１３年度

計算科学技術特論A

DDR3 DIMM

FX10の通信網（１TOFU単位）
１TOFU単位ノード

計算ノード内
ノード

ノード

ノード
ノード
ノード

ノード

ノード

1TOFU単位
間の結合用

ノード

ノード

ノード

10

ノード

ノード

２０１３年度

計算科学技術特論A

6本それぞれ5GB/秒（双方向）FX10の通信網（１TOFU単位間の結合）
３次元接続
１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

11

１
TOFU
単位

１１
TOFU
TOFU
単位
単位

１１
TOFU
TOFU
単位
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１１
TOFU
TOFU
単位
単位

１１
TOFU
TOFU
単位
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

１
TOFU
単位

２０１３年度 ユーザから見ると、X軸、Y軸、Z軸について、
奥の1TOFUと、手前の１TOFUは、繋がってみえます（３次元トーラス接続）
 ただし物理結線では
 X軸はトーラス Y軸はメッシュ
 Z軸はメッシュまたは、
トーラスになっています
計算科学技術特論A

バッチ処理とMPIジョブの投入

12

２０１３年度

計算科学技術特論A

FX10スーパコンピュータシステムでのジョブ実行形態の例
 以下の２通りがあります インタラクティブジョブ実行

ＰＣでの実行のように、コマンドを入力して実行する方法
 スパコン環境では、あまり一般的でない デバック用、大規模実行はできない
 ＦＸ１０では、以下に限定（東大基盤センターの運用方針）


１ノード（16コア）（２時間まで）
 ８ノード（128コア）（10分まで）
 バッチジョブ実行

バッチジョブシステムに処理を依頼して実行する方法
 スパコン環境で一般的
 大規模実行用 ＦＸ１０では、最大1440ノード（23,040コア）(６時間)

13

２０１３年度

計算科学技術特論A

バッチ処理とは スパコン環境では、通常は、インタラクティブ実行（コマンドラ

インで実行すること）はできません。
 ジョブはバッチ処理で実行します。
バッチキュー
バッチ処理
システムがジョブを取り出す

ジョブの依頼

実行
ユーザ
14

スパコン
２０１３年度

計算科学技術特論A

コンパイラの種類とインタラクティブ実行
およびバッチ実行の例（FX10）


インタラクティブ実行、およびバッチ実行で、利用するコンパイラ（Ｃ言語、Ｃ＋＋言語、Ｆｏｒｔｒａｎ９０言語）の種類が違います

インタラクティブ実行では
 オウンコンパイラ（そのノードで実行する実行ファイルを生成するコンパイラ）を使います

バッチ実行では クロスコンパイラ（そのノードでは実行できないが、バッチ実行する時のノードで実行できる実行ファイルを生成するコンパイラ）を使います

それぞれの形式
 オウンコンパイラ： ＜コンパイラの種類名＞
 クロスコンパイラ： ＜コンパイラの種類名＞px
 例）富士通Fortran90コンパイラ
オウンコンパイラ： frt
 クロスコンパイラ： frtpx


15

２０１３年度

計算科学技術特論A

バッチキューの設定のしかた（FX10の例）
 バッチ処理は、富士通社のバッチシステムで管理されている。
 以下、主要コマンドを説明します。
 ジョブの投入：

pjsub <ジョブスクリプトファイル名> -g <プロジェクトコード> 自分が投入したジョブの状況確認： pjstat
 投入ジョブの削除： pjdel <ジョブID>
 バッチキューの状態を見る： pjstat --rsc
 バッチキューの詳細構成を見る： pjstat –-rsc -x
 投げられているジョブ数を見る： pjstat --rsc -b
 過去の投入履歴を見る： pjstat --history
 同時に投入できる数／実行できる数を見る： pjstat --limit
16

２０１３年度

計算科学技術特論A

インタラクティブ実行のやり方の例（FX10スーパーコンピュータシステム）
 コマンドラインで以下を入力
 1ノード実行用$ pjsub --interact
 ８ノード実行用$ pjsub --interact -L “node=8”※インタラクティブ用のノード総数は５０ノードです。
もしユーザにより５０ノードすべて使われている場合、資源が空くまで、ログインできません。
17

２０１３年度

計算科学技術特論A

pjstat --rsc の実行画面例$ pjstat --rsc
RSCGRP
lecturelecture8

STATUS[ENABLE,START][DISABLE,STOP]

NODE:COORD
72:2x3x12
72:2x3x12

使える
キュー名
(リソース
グループ)

現在
使えるか

ノードの
物理構成情報18

２０１３年度

計算科学技術特論A

pjstat --rsc -x の実行画面例$ pjstat --rsc -x
RSCGRP
lecturelecture8

使える
キュー名
(リソース
グループ)

19

STATUSMIN_NODE MAX_NODE ELAPSE MEM(GB) PROJECT[ENABLE,START]1
1200:15:00
28
gt58[DISABLE,STOP]1
12
00:15:00
28
gt58

現在
使えるか

ノードの
実行情報

２０１３年度

計算科学技術特論A

課金情報
（財布）
実習では
１つのみ

pjstat --rsc -b の実行画面例$ pjstat --rsc –b
RSCGRP
lecturelecture8

STATUSTOTAL RUNNING QUEUED HOLD OTHER NODE:COORD[ENABLE,START] 0000072:2x3x12[DISABLE,STOP]
0
0
0
0
0
72:2x3x12

使える
現在
キュー名 使える(リソース かグループ)

20

ジョブ 実行
の
して
総数 いる
ジョブ
の
数
２０１３年度

待た
されて
いる
ジョブ
の
数
計算科学技術特論A

ノードの
物理
構成
情報

JOBスクリプトサンプルの説明（ピュアＭＰＩ）（hello-pure.bash, C言語、Fortran言語共通）
リソースグループ名：lecture#!/bin/bash#PJM -L "rscgrp=lecture"#PJM -L "node=12"#PJM --mpi "proc=192"#PJM -L "elapse=1:00"
mpirun ./hello

利用ノード数
利用コア数（ＭＰＩプロセス数）

MPIジョブを16＊12 = 192 プロセスで実行する。

21

２０１３年度

計算科学技術特論A

実行時間制限：１分

TOFU
Network

ピュアＭＰＩの実行状況（ノード内）
１ソケットのみ

ＭＰＩプロセス

各ＣＰＵの内部構成

20GB/秒

Core
#0

Core
#1

Core
#2

Core
#3

…

Core
#12

Core
#13

Core
#14

Core
#15

L1

L1

L1

L1

: L1データキャッシュ32KB

L1

L1

L1

L1

ICC

L2 (16コアで共有、12MB)

85GB/秒=(8Byte×1333MHz
×8 channel)

MemoryMemoryMemoryMemory4GB ×2枚

4GB ×2枚

4GB ×2枚

4GB ×2枚

DDR3 DIMM

ノード内合計メモリ量：8GB×4＝32GB

FX10計算ノードの構成
22

２０１３年度

計算科学技術特論A

並列版Helloプログラムを実行しよう（ピュアMPI）
1.
2.
3.

4.
5.

23

Helloフォルダ中で以下を実行する$ pjsub hello-pure.bash
自分の導入されたジョブを確認する$ pjstat
実行が終了すると、以下のファイルが生成される
hello-pure.bash.eＸＸＸＸＸＸ
hello-pure.bash.oＸＸＸＸＸＸ （XXXXXXは数字）
上記の標準出力ファイルの中身を見てみる$ cat hello-pure.bash.oＸＸＸＸＸＸ“Hello parallel world!”が、16プロセス*１２ノード=１９２表示されていたら成功。
２０１３年度

計算科学技術特論A

バッチジョブ実行による標準出力、
標準エラー出力
 バッチジョブの実行が終了すると、標準出力ファイルと標

準エラー出力ファイルが、ジョブ投入時のディレクトリに作成されます。
 標準出力ファイルにはジョブ実行中の標準出力、標準エ
ラー出力ファイルにはジョブ実行中のエラーメッセージが出力されます。

ジョブ名.oXXXXX --- 標準出力ファイル
ジョブ名.eXXXXX --- 標準エラー出力ファイル(XXXXX はジョブ投入時に表示されるジョブのジョブID)

24

２０１３年度

計算科学技術特論A

並列版Helloプログラムを実行しよう（ハイブリッドMPI）
1.
2.
3.

4.
5.
25

Helloフォルダ中で以下を実行する$ pjsub hello-hy16.bash
自分の導入されたジョブを確認する$ pjstat
実行が終了すると、以下のファイルが生成される
hello-hy16.bash.eＸＸＸＸＸＸ
hello-hy16.bash.oＸＸＸＸＸＸ （XXXXXXは数字）
上記標準出力ファイルの中身を見てみる$ cat hello-hy16.bash.oＸＸＸＸＸＸ“Hello parallel world!”が、1プロセス*12ノード=12 個表示されていたら成功。
２０１３年度

計算科学技術特論A

JOBスクリプトサンプルの説明（ハイブリッドＭＰＩ）（hello-hy16.bash, C言語、Fortran言語共通）#!/bin/bash#PJM -L "rscgrp=lecture"#PJM -L "node=12"#PJM --mpi "proc=12"#PJM -L "elapse=1:00"
export OMP_NUM_THREADS=16
mpirun ./hello
MPIジョブを1＊12 = 12 プロセスで実行する。

26

２０１３年度

計算科学技術特論A

リソースグループ名：lecture
利用ノード数
利用コア数（ＭＰＩプロセス数）
実行時間制限：１分
１ＭＰＩプロセス当たり
16スレッド生成

TOFU
Network

ピュアＭＰＩの実行状況（ノード内）
１ソケットのみ
各ＣＰＵの内部構成

20GB/秒

Core
#0

Core
#1

Core
#2

Core
#3

…

Core
#12

Core
#13

Core
#14

Core
#15

L1

L1

L1

L1

: L1データキャッシュ32KB

L1

L1

L1

L1

ICC

L2 (16コアで共有、12MB)

85GB/秒=(8Byte×1333MHz
×8 channel)

MemoryMemoryMemoryMemory4GB ×2枚

4GB ×2枚

4GB ×2枚

4GB ×2枚

ＭＰＩプロセス

ノード内合計メモリ量：8GB×4＝32GB

スレッド
27

DDR3 DIMM

FX10計算ノードの構成
２０１３年度

計算科学技術特論A

その他の注意事項（その１） MPI用のコンパイラを使うこと
 MPI用のコンパイラを使わないと、MPI関数が未定義というエラーが出て、コンパイルできな
くなる 例えば、以下のコマンド
 Fortran90言語： mpif90

 C言語： mpicc
 C++言語: mpixx, mpic++ コンパイラオプションは、逐次コンパイラと同じ
28

２０１３年度

計算科学技術特論A

その他の注意事項（その２） ハイブリッドMPIの実行形態

MPIプロセス数＋OpenMPスレッド数 <= 利用コア総数 HT（Intel）やSMT(IBM)などの、物理コア数の定数倍のスレッドが実行できるハードの場合
 スレッド数（論理スレッド数）が上記の利用コア総数

 必ずしも、１ノード内に１MPIプロセス実行が高速とはならない 一般に、OpenMPによる台数効果が8スレッドを超えると悪くなるため。
 効率の良いハイブリッドMPI実行には、効率の良いOpenMP実装が必須
29

２０１３年度

計算科学技術特論A

MPI実行時のリダイレクトについて
 一般に、スーパーコンピュータでは、

MPI実行時の入出力のリダイレクトができません ×例）mpirun ./a.out < in.txt > out.txt
 専用のリダイレクト命令が用意されています。
 FX10でリダイレクトを行う場合、以下のオプションを指定します。
 ○例) mpirun --stdin ./in.txt --ofout out.txt ./a.out

30

２０１３年度

計算科学技術特論A

並列処理の評価指標：
弱スケーリングと強スケーリング

31

２０１３年度

計算科学技術特論A

弱スケーリング (Weak Scaling) ノードあたりの問題サイズを固定し、並列処理時の全体の問題サイズを増加することで、性能評価をする方法
 問題サイズN ときの計算量がO(N )である場合、並列処理のノード数が増加しても、理想的な実行時間は変わらないと期待できる
一般的にノード数が増加すると通信時間が増大するため、そうはならない 該当する処理は


陽解法のシミュレーション全般
 陰解法で、かつ連立一次方程式の解法に反復解法を用いているシミュレーション

１ノードあたりの
問題サイズ
８ノード実行での問題サイズ
32

２０１３年度

６４ノード実行での問題サイズ
計算科学技術特論A

強スケーリング (Strong Scaling) 全体の問題サイズを固定し、ノード数を増加することで性能評価をする方法
 理想的な実行時間は、ノード数に反比例して減少する。
一般的にノード数が増加すると１ノードあたりの問題サイズが減少し、通信時間の占める割合が増大するため、理想的な実行時間の減少は限定
 該当する処理は

計算量が膨大なアプリケーション 例えば、連立一次方程式の解法。データ量O (N 2 )に対して、計算量はO ( N 3 )


固定した
問題サイズ
33

８ノード実行での問題サイズ
２０１３年度

６４ノード実行での問題サイズ
計算科学技術特論A

弱スケーリングと強スケーリング
適用アプリの特徴
 弱スケーリングが適用できるアプリケーションは、原理的に通信が少ないアプリケーション
領域分割法などにより、並列化できるアプリケーション
 主な通信は、隣接するプロセスしかない ノード数を増すことで、実行時間の面で容易に問題サイズを大規模化
 通信時間の占める割合は超並列実行でも少ないアプリケーション
 強スケーリングを適用しないといけないアプリケーションは、計算量が膨大になるアプリケーション
全体の問題サイズは、実行時間の制約から大規模化できない そのため、１ノードあたりの問題サイズは、ノード数が多い状況で小さくなる
 その結果、通信処理の占める時間がほとんどになる 超並列実行時で通信処理の最適化が重要になるアプリケーション


34

２０１３年度

計算科学技術特論A

強スケールアプリケーションの問題
 TOP500で採用されているLINPACK

密行列に対する連立一次方程式の解法のアプリケーション
 ２０１２年１１月のTOP500の問題サイズ
 （２位）Sequoia、
N=12,681,215、#cores=1,572,864、N/#cores=8.06
 （３位）K-computer、
N=11,870,208、#cores=705,024、N/#cores=16.8
 （６位）SuperMUC、
N=5,201,920、#cores=147,456、N/#cores=35.2
 上位のマシンほど、コア当たりの問題サイズが小さい←通信時間の占める割合が大きくなりやすい 今後コア数が増加すると、通信時間の削減が問題になる
35

２０１３年度

計算科学技術特論A

ピュアMPIプログラム開発の基礎

36

２０１３年度

計算科学技術特論A

MPI並列化の大前提（再確認） ＳＰＭＤ 対象のメインプログラムは、
 すべてのコア上で、かつ、
 同時に起動された状態から処理が始まる。
 分散メモリ型並列計算機
 各プロセスは、完全に独立したメモリを持っている。（共有メモリではない）
37

２０１３年度

計算科学技術特論A

並列化の考え方（Ｃ言語） SIMDアルゴリズムの考え方（４プロセスの場合）
各ＰＥで
行列Ａ
重複して
for ( j=0; j<n/4; j++) { 内積( j, i ) }

所有する

プロセス0
for ( j=0; j<n; j++){ 内積( j, i ) }for ( j=n/4; j<(n/4)*2; j++) { 内積( j, i ) }
プロセス１
for ( j=(n/4)*2; j<(n/4)*3; j++) { 内積( j, i ) }プロセス2
for ( j=(n/4)*3; j<n; j++) { 内積( j, i ) }
プロセス3

38

２０１３年度

計算科学技術特論A

ベクトルｘ

並列化の考え方（Fortran言語）
 SIMDアルゴリズムの考え方（４プロセスの場合）
行列Ａ
プロセス0
do j=1, n
内積( j, i )enddo

do j=1, n/4
内積( j, i )enddo

各プロセスで重複して所有する

do j=n/4+1, (n/4)*2
内積( j, i )プロセス１ enddo
do j=(n/4)*2+1, (n/4)*3
内積( j, i )プロセス2 enddo
do j=(n/4)*3+1, n
内積( j, i )プロセス3 enddo

39

２０１３年度

計算科学技術特論A

ベクトルｘ

初心者が注意すること
 各プロセスでは、独立した配列が個別に確保されます。
ＰＥ０
Ａ［N］［N］

ＰＥ１

ＰＥ２

Ａ［N］［N］

Ａ［N］［N］

ＰＥ３
Ａ［N］［N］

 myid変数は、MPI_Init()関数が呼ばれた段階で、各プロセス固有の値になっています。
ＰＥ０

myid = 0

40

ＰＥ１

ＰＥ２

myid = １

２０１３年度myid = 2

計算科学技術特論A

ＰＥ３
myid = 3

並列プログラム開発の指針
1.
2.
3.
4.
5.
6.

正しく動作する逐次プログラムを作成する
１．のプログラムで、適切なテスト問題を作成する
２．のテスト問題の実行について、適切な処理の単位ごとに、正常動作する計算結果を確認する
１．の逐次プログラムを並列化し、並列プログラミングを行う
２．のテスト問題を実行して動作検証する
このとき３．の演算結果と比較し、正常動作をすることを確認する。もし異常であれば、４．に戻りデバックを行う。

41

２０１３年度

計算科学技術特論A

数値計算プログラムの特徴を利用し並列化がなされる
 数値計算プログラムの処理単位は、プログラム上の

基本ブロック（ループ単位など）ではなく、数値計算上の処理単位（数式レベルで記述できる単位）となる


離散化（行列作成）部分、LU分解法部分（LU分解部分、
前進代入部分、後退代入部分）、など

 演算結果は、なんらかの数値解析上の意味において検証

理論解とどれだけ離れているか、考えられる丸め誤差の範囲内にあるか、など
 計算された物理量（例えば流速など）が物理的に妥当な範囲内にあるか、など
 両者が不明な場合でも、数値的に妥当であると思われる逐次の結果と比べ、並列化した結果の誤差が十分に小さいか、など42

２０１３年度

計算科学技術特論A

並列化の方針の例（Ｃ言語）
全プロセスで行列AをN×Nの大きさ、ベクトルｘ、ｙをNの大きさ、確保してよいとする。
各プロセスは、担当の範囲のみ計算するように、
ループの開始値と終了値を変更する。

1.
2.


ブロック分散方式では、以下になる。
（n が numprocs で割り切れる場合）

ib = n / numprocs;for ( j=myid*ib; j<(myid+１)*ib; j++){…}（２の並列化が完全に終了したら）各プロセスで担当のデータ部分しか行列を確保しないように変更する。

3.


43

上記のループは、以下のようになる。
for ( j=0; j<ib; j++) { … }２０１３年度

計算科学技術特論A

並列化の方針の例（Fortran言語）
全プロセスで行列AをN×Nの大きさ、ベクトルｘ、ｙをNの大きさ、確保してよいとする。
各プロセスは、担当の範囲のみ計算するように、
ループの開始値と終了値を変更する。

1.
2.


ブロック分散方式では、以下になる。
（n が numprocs で割り切れる場合）

ib = n / numprocs
do j=myid*ib+１, (myid+１)*ib … enddo（２の並列化が完全に終了したら）各プロセスで担当のデータ部分しか行列を確保しないように変更する。

3.


44

上記のループは、以下のようになる。
do j=１, ib … enddo
２０１３年度

計算科学技術特論A

データ分散方式に関する注意 負荷分散を考慮し、多彩なデータ分散方式を採用可能
 しかし、数学的に単純なデータ分散方式が良い ◎：ブロック分散、サイクリック分散（ブロック幅＝１）
 △～〇：ブロック・サイクリック分散（ブロック幅＝任意）
 理由：
 複雑な（一般的な）データ分散は、各MPIプロセスが所有するデータ分散情報（インデックスリスト）を必要とするため、メモリ量が余分に必要なる
 例：１万並列では、少なくとも１万次元の整数配列が必要
 数学的に単純なデータ分散の場合は、インデックスリストは不要
45

２０１３年度

計算科学技術特論A

並列化の方針（行列-ベクトル積)（Ｃ言語） 全PEでN×N行列を持つ場合
PE0

for ( j=(n/4)*2; j<(n/4)*3; j++) { 内積( j, i ) }for ( j=0; j<(n/4); j++) { 内積( j, i ) }
PE2

PE１

for ( j=(n/4)*3; j<n; j++) { 内積( j, i ) }for ( j=(n/4); j<(n/4)*2; j++) { 内積( j, i ) }PE3※各PEで使われない領域が出るが、担当範囲指定がしやすいので実装がしやすい。
46

２０１３年度

計算科学技術特論A

並列化の方針（行列-ベクトル積）（Fortran 言語）
 全PEでN×N行列を持つ場合
do j=(n/4)*2+1, (n/4)*3
内積( j, i )
enddo

PE0
do j=1, n/4
内積( j, i )
enddo

PE2

do j=(n/4)*3+1, n
内積( j, i )
enddo

PE１
do j=n/4+1, (n/4)*2
内積( j, i )
enddo

PE3

※各PEで使われない領域が出るが、担当範囲指定がしやすいので実装がしやすい。
47

２０１３年度

計算科学技術特論A

並列化の方針（行列-ベクトル積） この方針では、ｙ＝Ａｘのベクトルｙは、以下のように一部分しか計算されないことに注意！
PE0
＝

＝

＝

PE2

＝

PE１

PE3
48

２０１３年度

計算科学技術特論A

並列化の方針のまとめ
 行列全体（A[N][N]）を各プロセスで確保することで、SIMDの考え方を、逐次プログラムに容易に適用できる
ループの開始値、終了値のみ変更すれば、並列化が完成する
 この考え方は、MPI、OpenMPに依存せず、適用できる。
 欠点


 最大実行可能な問題サイズが、利用ノード数によらず、１ノードあたりのメモリ量で制限される（メモリに関するスケーラビリティが無い） ステップ２のデバックの困難性を低減できる

完全な並列化（ステップ２）の際、ステップ１での正しい計算結果を参照できる 数値計算上の処理単位ごとに、ステップごとに完全並列化ができる（モジュールごとに、デバックできる）


49

２０１３年度

計算科学技術特論A

行列‐ベクトル積のピュアMPI並列化の例（C言語）
ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);ierr = MPI_Comm_size(MPI_COMM_WORLD, &numprocs);…
ib = n/numprocs;
ブロック分散を仮定した
jstart = myid * ib;
担当ループ範囲の定義
jend = (myid+1) * ib;
if ( myid == numprocs-1) jend=n;for( j=jstart; j<jend; j++) {y[ j ] = 0.0;
for(i=0; i<n; i++) {y[ j ] += A[ j ][ i ] * x[ i ];}}50

MPIプロセスの担当ごとに縮小したループの構成

２０１３年度

計算科学技術特論A

行列‐ベクトル積のピュアMPI並列化の例（Fortran言語）
call MPI_INIT(ierr)call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)…
ib = n/numprocs
jstart = 1 + myid * ib
ブロック分散を仮定した
jend = (myid+1) * ib
担当ループ範囲の定義
if ( myid .eq. numprocs-1) jend = n
do j = jstart, jendy( j ) = 0.0d0
do i=1, ny( j ) = y( j ) + A( j, i ) * x( i )
enddo
enddo
51

２０１３年度MPIプロセスの担当ごとに縮小したループの構成

計算科学技術特論A

nがMPIプロセス数で割切れない時
 nがプロセス数のnumprocsで割り切れない場合
 配列確保：A（N/numprocs + mod(N, numprocs), N）
 ループ終了値：numprocs-1のみ終了値がnとなるように実装
ib = n / numprocs;
if ( myid == (numprocs - １) ) {i_end = n;
} else {i_end = (myid+１)*ib;}for ( i=myid*ib; i<i_end; i++) { … }52

２０１３年度

計算科学技術特論A

余りが多い場合
 mod(N, numprocs)が大きいと、負荷バランスが悪化
 例：N=10、numprocs=6
 int(10/6)=1なので、プロセス0～5は1個のデータ、プロセス6は4個のデータを持つ
 各プロセスごとの開始値、終了値のリストを持てば改善可能
 プロセス0： i_start(0)=1, i_end(0)=2, 2個
 プロセス1: i_start(1)=3, i_end(1)=4, 2個
 プロセス2: i_start(2)=5, i_end(2)=6, 2個
 プロセス3: i_start(3)=7, i_end(3)=8, 2個
 プロセス4: i_start(4)=9, i_end(4)=9, 1個
 プロセス5: i_start(5)=10, i_end(5)=10, 1個
 欠点：プロセス数が多いと、上記リストのメモリ量が増える
53

２０１３年度

計算科学技術特論A

ハイブリットMPIプログラム開発の基礎

54

２０１３年度

計算科学技術特論A

用語の説明 ピュアMPI実行

並列プログラムでMPIのみ利用
 MPIプロセスのみ ハイブリッドMPI実行

並列プログラムでMPIと何か（X（エックス））を利用
 MPIプロセスと何か（X）の混合
 何か（X）は、OpenMPによるスレッド実行が主流
 「MPI+X」の実行形態

上記のハイブリッドMPI実行と同義として使われる
 Xは、OpenMPや自動並列化によるスレッド実行、CUDAなどのGPU向き実装、OpenACCなどのGPUやメニーコア向き実行、などの組合せがある。今後主流となる形態で変わる。
55

２０１３年度

計算科学技術特論A

ハイブリッドMPI実行の目的
 同一の資源量（総コア数）の利用に対し ピュアMPI実行でのMPIプロセス数に対し、

ハイブリッドMPI実行でのMPIプロセス数を減らすこと
で、通信時間を削減することが主な目的 例）東京大学のFX10
 全系は4,800ノード、76,800コア
 ピュアMPI実行：76,800プロセス実行
 ハイブリッドMPI実行（1ノード16スレッド実行）：4,800プロセス
 プロセス数の比は16倍！
56

２０１３年度

計算科学技術特論A

ハイブリッドMPI並列プログラム
開発の指針
正しく動作するピュアMPIプログラムを開発する
OpenMPを用いて対象カーネルをスレッド並列化する
2．の性能評価をする
3．の評価結果から性能が不十分な場合、対象カーネルについてOpenMPを用いた性能チューニングを行う。
3．へ戻る。
全体性能を検証し、通信時間に問題がある場合、通信処理のチューニングを行う。

1.
2.
3.
4.

5.

57

２０１３年度

計算科学技術特論A

行列‐ベクトル積のハイブリッドMPI並列化の例（C言語）
ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);ierr = MPI_Comm_size(MPI_COMM_WORLD, &numprocs);…
ib = n/numprocs;
ブロック分散を仮定した
jstart = myid * ib;
担当ループ範囲の定義
jend = (myid+1) * ib;
if ( myid == numprocs-1) jend=n;#pragma omp parallel for private(i)この一文を追加するだけ！
for( j=jstart; j<jend; j++) {y[ j ] = 0.0;
MPIプロセスの担当ごとに
for(i=0; i<n; i++) {縮小したループの構成
y[ j ] += A[ j ][ i ] * x[ i ];}}58

２０１３年度

計算科学技術特論A

行列‐ベクトル積のハイブリッドMPI並列化の例（Fortran言語）
call MPI_INIT(ierr)call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)…
ib = n/numprocs
jstart = 1 + myid * ib
ブロック分散を仮定した
jend = (myid+1) * ib
担当ループ範囲の定義
if ( myid .eq. numprocs-1) jend = n!$omp parallel do private(i)この文を追加するだけ！
do j = jstart, jendy( j ) = 0.0d0
do i=1, n
MPIプロセスの担当ごとにy( j ) = y( j ) + A( j, i ) * x( i )縮小したループの構成
enddo
enddo
!$omp end parallel do
59

２０１３年度

計算科学技術特論A

ハイブリッドMPI実行の注意点（その1） ハイブリッドMPI実行では、MPIプロセス数に加えて、スレッド数がチューニングパラメタとなり、複雑化する。
例）１ノード16コア実行

２MPIプロセス、８スレッド実行

４MPIプロセス、４スレッド実行

コ コ コ コ コ コ
コ コ コ コ コ コ コ コ コ コ
ア ア ア ア ア ア
ア ア ア ア ア ア ア ア ア ア
1 1 1 1 1 1
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5

コ コ コ コ コ コ
コ コ コ コ コ コ コ コ コ コ
ア ア ア ア ア ア
ア ア ア ア ア ア ア ア ア ア
1 1 1 1 1 1
0 1 2 3 4 5 6 7 8 9
0 1 2 3 4 5

１つのMPI
プロセス
の割り当て
対象

 ccNUMAの計算機では、ソケット数ごとに１MPIプロセス実行が高速となる可能性がある（ハードウェア的に）
例）T2K（AMD Quad Core Opteron）、4ソケット、16コア
４MPIプロセス、４スレッド実行
コ コ コ コア ア ア ア0 1 2 3ソケット060

コ コ コ コア ア ア ア4 5 6 7

ソケット1

コ コ
コ コ
ア ア
ア ア
1 1
8 9
0 1

ソケット2
２０１３年度コ コ コ コア ア ア ア1 1 1 12 3 4 5

ソケット3
計算科学技術特論A

ハイブリッドMPI実行の注意点（その２） ハイブリッドMPI実行の実行効率を決める要因

ハイブリッドMPI化による通信時間の削減割合
OpenMP等で実現される演算処理のスレッド実行効率

1.
2.

 特に、2は注意が必要。

単純な実装だと、経験的に8スレッド並列を超えると、スレッド実行
時の台数効果が劇的に悪くなる。
 効率の良いスレッド並列化の実装をすると、ハイブリッドMPI実行時に効果がより顕著になる。
 実装の工夫が必要。たとえば
1.
2.
3.

61

ファーストタッチ（すでに説明済み）の適用
メモリ量や演算量を増加させても、スレッドレベルの並列性を増加させる
アンローリングなどの逐次高速化手法を、スレッド数に特化させる
２０１３年度

計算科学技術特論A

ハイブリッドMPI実行の注意点（その３） 通信処理の時間に含まれる、データのコピー時間が、

通信時間よりも大きいことがある
問題空間の配列から送信用の配列にコピーする処理（パッキング） 受信用の配列から問題空間の配列へコピーする処理（アンパッキング） 上記のコピー量が多い場合、コピー操作自体もOpenMP化すると高速化される場合がある。
 反面、ハードウェアによっては、OpenMP化すると遅くなる。
このときは、逐次処理にしないといけない。
 パッキング、アンパッキングをOpenMP化する／しない、もハイブリッドMPI実行では重要なチューニング項目になる
62

２０１３年度

計算科学技術特論A

ハイブリッドMPIの起動方法
 スパコンごとに異なるが、以下の方法が主流（すでに説明済み）。
1.
2.

バッチジョブシステムを通して、MPIの数を指定
実行コマンドで、OMP_NUM_THREADS環境変数でスレッド数を指定

 ccNUMAの場合、MPIプロセスの割り当てを、期待する

物理ソケットに割り当てないと、ハイブリッドMPI実行の効果が無くなるLinuxでは、numactlコマンドで実行時に指定する
 スパコン環境によっては、プロセスを指定する物理コアに割り当てる方法がある。（各スパコンの利用マニュアルを参考）

63

２０１３年度

計算科学技術特論A

数値計算ライブラリとハイブリットMPI実行
 数値計算ライブラリのなかには、ハイブリッドMPI実行をサポートしているものがある


数値計算ライブラリがスレッド並列化されている場合

 特に、密行列用ライブラリのScaLAPACKは、

通常、ハイブリッドMPI実行をサポート
ScaLAPACKは、MPI実行をサポート
 ScaLAPACKは、逐次のLAPACKをもとに構築 LAPACKは基本数値計算ライブラリBLASをもとに構築
 BLASは、スレッド実行をサポート⇒BLASレベルのスレッド実行と、
ScaLAPACKレベルのMPI実行を基にした
ハイブリッドMPI実行が可能

64

２０１３年度

計算科学技術特論A

スレッド並列版BLAS利用の注意
 BLASライブラリは、スレッド並列化がされている
 利用方法は、OpenMPを用いた並列化と同じ


OMP_NUM_THREADSで並列度を指定

 BLASで利用するスレッド数が利用可能なコア数を超えると、動かないか、動いたとしても速度が劇的に低下する
 BLASを呼び出す先がスレッド並列化をしている場合、BLAS内でスレッド並列化をすると、総合的なスレッド数が、利用可能なコア数を超えることがある。このため、速度が劇的に低下する。
 一般的に、逐次実行の演算効率が、スレッド並列の実行効率に比べ良い

上位のループをスレッド並列化し、そのループから逐次BLASを呼び出す実装がよい

65

２０１３年度

計算科学技術特論A

逐次BLASをスレッド並列化して呼び出す例
 通常のBLASの呼び出し

do i=1, Ak
call dgemm(…) ←スレッド並列版BLASを呼び出し(コンパイラオプションで指定)

enddo 上位のループで並列化したBLASの呼び出し!$omp parallel dodo i=1, Ak
call dgemm(…) ←逐次BLASを呼び出し（コンパイラオプションで指定）

enddo!$omp end parallel do
66

２０１３年度

計算科学技術特論A

＜スレッド並列版BLAS＞と＜逐次BLASを上位のループでスレッド並列呼び出し＞する時の性能例
 T2Kオープンスパコン（東大版）

AMD Quad Core Opteron 1ノード（16コア）を利用
 日立製作所によるCコンパイラ（日立最適化C）
 OpenMP並列化を行った
 最適化オプション：“-Os -omp”

BLAS GOTO BLAS ver.1.26（スレッド並列版，および遂次版の双方）
 対象処理高精度行列‐行列積の主計算
 複数の行列‐行列積(dgemm呼び出し)を行う部分
67

２０１３年度

計算科学技術特論A

n=1000での性能（ T2K(1ノード, 16コア））
BLAS内でスレッド並列化する場合に対する速度向上[速度向上]
4.5
4

3.86

3.5

3.80 3.84
3.59

3.79

3.53

3.78
3.50

3.44
3

2.64

2.5

2.59

2.32
1.80

2

1.97

1.77

Speedup to isw=1

1.5
18スレッドを超えると約3.8倍の速度向上！

0.81

0.5

[スレッド数]

0
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

ScaLAPACKにおけるハイブリッドMPI実行の効果の例
 ScaLAPACKの連立一次方程式解法ルーチン

PDGESV 東京大学情報基盤センターのHITACHI SR16000
 IBM Power7 (3.83GHz)
 1ノード４ソケット、１ソケットあたり8コア、合計32コア、980.48GFLOPS／ノード
 SMT利用で、１ノード６４論理スレッドまで利用可能
 ScaLAPACKは、同環境で提供されているIBM社のESSL(Engineering and Scientific Subroutine Library)
ライブラリを利用

69

２０１３年度

計算科学技術特論AScaLAPACKにおけるハイブリッドMPI実行の効果の例
SR16000の2ノードでの実行 (問題サイズN=32,000)
902.2

947.7
GFLOPS

880.5

833.4

788.5
最適な
ハイブリッド
MPI実行

712.9

５％の高速化ピュアMPI
実行

ノードあたりの構成
P（MPIプロセス数）×T（スレッド数）

P64×T1
70

P32×T2

680.3

ノード内は
全て
スレッド実行

ハイブリッドMPI実行の組合せ
P16×T4

P8×T8

２０１３年度

P4×T16

計算科学技術特論A

P2×T32

P1×T64

コンパイラ最適化の影響（その１） MPI化、および、OpenMP化に際して、ループ構造を逐次から変更することになる この時、コンパイラに依存し、コード最適化が並列ループに対して、効かない（遅い）コードを生成することがある 上記の場合、逐次実行での効率に対して、並列実行での効率が低下し、台数効果の向上を制限する
 たとえば、ループ変数に大域変数を記載すると、コンパイラの最適化を阻害することがある

71

特に並列処理制御変数である、全体のMPIプロセス数を管理する変数、自分のランク番号を管理する変数は、大域変数であることが多いので注意。
２０１３年度

計算科学技術特論A

コンパイラ最適化の影響（その２） MPI並列コードで、ループに大域変数を使っている例

C言語の例

Fortran言語の例

ib = n/numprocs;for( j= myid * ib; j<(myid+1) * ib; j++) {y[ j ] = 0.0;
for(i=0; i<n; i++) {y[ j ] += A[ j ][ i ] * x[ i ];}}

ib = n/numprocs
do j = 1 + myid * ib, (myid+1) * iby( j ) = 0.0d0
do i=1, ny( j ) = y( j ) + A( j, i ) * x( i )
enddo
enddo

 上記のmyidは大域変数で、自ランク番号を記憶している変数
 コンパイラがループ特徴を把握できず、最適化を制限
←逐次コードに対して、演算効率が低下し、台数効果を制限

 解決策：局所変数を宣言しmyidを代入。対象を関数化。
72

２０１３年度

計算科学技術特論A

ハイブリッドMPIプログラミングのまとめ ノード数が増えるほど、ピュアMPI実行に対する効果が増加

経験的には、１０００MPIプロセスを超える実行で、ハイブリッドMPI実行が有効となる
 現状での効果はアプリケーションに依存するが、経験的には数倍（２～３倍）高速化される
 現在、多くの実例が研究されている

エクサに向けて10万並列を超える実行では、
おそらく数十倍の効果が期待される

 ノードあたりの問題サイズが小さいほど、

ハイブリッドMPI実行の効果が増大

73

弱スケーリングより強スケーリングのほうがハイブリッドMPI実行の効果がある
２０１３年度

計算科学技術特論A

レポート課題（その１） 問題レベルを以下に設定
問題のレベルに関する記述：•L00: きわめて簡単な問題。
•L10： ちょっと考えればわかる問題。
•L20： 標準的な問題。
•L30： 数時間程度必要とする問題。
•L40： 数週間程度必要とする問題。複雑な実装を必要とする。
•L50： 数か月程度必要とする問題。未解決問題を含む。
※L４０以上は、論文を出版するに値する問題。

 教科書のサンプルプログラムは以下が利用可能


Samples-fx.tar



Mat-Vec-fx.tar



PowM-fx.tar



Mat-Mat-fx.tar



Mat-Mat-d-fx.tar



LU-fx.tar

74

２０１３年度

計算科学技術特論A

レポート課題（その２）[L20] 使える並列計算機環境で、教科書のサンプルプログラムを並列化したうえで、
ピュアＭＰＩ実行、および、ハイブリッドＭＰＩ実行で性能が異なるか、実験環境（たとえば、12ノード、192コア）を駆使して、性能評価せよ。

1.



1ノードあたり、12MPI実行、1MPI+16スレッド実行、2MPI+8スレッド
実行、4MPI+4スレッド実行など、組み合わせが多くある。

75

２０１３年度

計算科学技術特論A

レポート課題（その３）2.
3.

[L10] ハイブリッドMPI実行がピュアMPI実行に対して有効となるアプリケーションを、論文等で調べよ。
[L20～] 自分が持っている問題に対し、ハイブリッドMPI実
行ができるようにプログラムを作成せよ。また、実験環境を用いて、性能評価を行え。

76

２０１３年度

計算科学技術特論A