5
SC’11 Special Issue

SC’11 特集号
ゴードンベル賞
特別賞
奨励賞

Graph 500 ランキング 3位
テクニカル・ペーパー 3件

http://www.gsic.titech.ac.jp/TSUBAME_ESJ

SC’11 特集号

03

09

ACM ゴードンベル賞・特 別 賞 ＆ SC'11テクニカル・ペーパー

TSUBAME 2.0スパコンにおける樹 枝 状 凝 固 成 長のフェーズフィールド法を用いたペタスケール・シミュレーション
下川辺隆史
遠藤敏夫

青木尊之
丸山直也

高木知弘
松岡 聡

山中晃徳

額田 彰

ACM ゴードンベル賞・奨 励 賞

TSUBAME2 における大 規 模 生 体 流 体 力 学シミュレーション
Massimo Bernaschi
Massimiliano Fatica

Mauro Bisson
松岡 聡

遠藤 敏夫
Simone Melchionna

Sauro Succi

14

18
23

Graph50 0 ランキング 3 位

大 規 模グラフ処 理ベンチマークGraph500 のTSUBAME 2.0 における挑 戦
鈴村 豊太郎

上野 晃司

SC'11テクニカル・ペーパー

Physis: ヘテロジニアススパコン向けステンシル
計 算フレームワーク
丸山 直也

野村 達男

佐藤 賢斗

松岡 聡

SC'11テクニカル・ペーパー（ 最 高 得点 獲 得 ）F T I : ヘテロジニアススパコン向け耐 障 害インタフェース〜 100TFlops超 東北地方太平洋沖地震シミュレーション 〜
Leonardo Bautista-Gomez
Franck Cappello松岡 聡

Dimitri Komatitch
中村 武

丸山 直也

坪井 誠司

TSUBAME 2.0スパコンにおける樹枝状凝固成長のフェーズフィールド法を用いた
ペタスケール・シミュレーション
下川辺隆史 * 青木尊之 **遠藤敏夫 **    丸山直也 **高木知弘 ***
松岡 聡 **

山中晃徳 ****

額田 彰 *** 東京工業大学大学院・総合理工学研究科 ** 東京工業大学学術国際情報センター*** 京都工芸繊維大学大学院・工芸科学研究科 **** 東京工業大学大学院・理工学研究科

軽 量・高強 度な新材料の開発は低 炭素社 会の実現に向けて非常に重要である。材料強 度を左右する材料のミクロな組 織は凝固過程で決定されるが、機械的強度を判定するには数mmというマクロなスケールまでの大規模計算を必要とする。
フェーズフィー
ルド法というメソスケールのモデルを用い、
TSUBAME 2.0 で合 金の樹枝状凝固の大規模シミュレーションを行った。CUDAを用いて有限差分法で離散化された時間発展方程式を解き、領域分割により並列化することで TSUBAME 2.0の4000 GPUを用いて、2.0
PFLOPS（単精度）という高い実行性能を得た。この成果に対し 、2011年 の ACMゴードンベル賞・特別賞（本賞）が与えられた 。

はじめに1模計算する。CUDAを用いて有限差分法で離散化された時間発展方
程式をプログラミングし、MPIを用い領域分割法で並列化することで
TSUBAME 2.0 のほぼ全ての計算資源を使った大規模計算を行う。フェー
ズフィールド法による超大規模計算を行った例は国内外で報告されておらず、本計算結果が材料科学分野に与えるインパクトは大きい。

軽量・高強度な新材料を開発することにより、物資を高効率（低燃費）で輸送することができ、低炭素社会の実現に大きく貢献する。金属材
料の強度はミクロな材料組織に強く依存し、そのミクロ組織は凝固過程において形作られる（図１）。一方、現実に機械的強度を判定するには、

2

フェーズフィールド法

数 mmのマクロなスケールでの解析が必要となる。ミクロな凝固ダイナ
ミクスを解明するために、非平衡統計力学から導出されるフェーズフィー
ルド法 [1] が近年注目されている。導出される方程式は時間空間の偏微

フェーズフィールド法は非平 衡統 計 物理学から導出され、分子スケー

分方程式になっていて、有限差分法や有限要素法などの格子法で解かルとマクロなスケールの中間のメソスケールの現 象を記 述することがれることが多い。フェーズフィールド法は複雑な非線形項を多く含み、スできる。秩 序 変 数 （フェーズフィールド変 数）を導入し、固相部分にテンシル計算で解く場合には1格子点あたりの演算量が多くなる。さらを、液相部分にを設 定する。界面を含む領域はが0かに、非常に狭い界面領域に複数の格子を含む必要があるため、格子サイ

ら1へと急峻かつ滑らかに変 化する拡 散界面として扱われ、ズを小さくする必要があり、格子数が膨大になり時間ステップを小さく取が界面位置となる。フェーズフィールド法では従 来から使われている

らなければならない。このため、CPUで計算すると時間がかかり過ぎるた

界面追 跡法 等の手法が不要となり、領域全体で同一の計算を行うこめ、これまでは主に2 次元計算による解析が行われてきた。

とができる。

本研究では、フェーズフィールド法に基づいて二元合金の一方向凝

本 研 究で対 象とする二元 合 金の樹 枝 状 凝 固 成 長では、フェーズ

固における樹 枝 状（デンドライト）組 織の成長を複 数 GPUにより大 規

フィールド方程 式と溶質の拡散 方程 式を解く。フェーズフィールドに対しては、界面エネルギーの異方性を考慮した次の方程式を解く。ここでは勾 配 係 数 、

ルギー障壁、は融 解エントロピー、

また、関数

と

は、
それぞれ

を用いている。
図 1 材料 組 織のミクロな組 織とデンドライト

03

はモビリティーであり、はエネ

は過冷 却度を表している。
と

SC’11 特集号

一方、溶質の拡散 方程式は、次の方程式を用いる。

フェーズフィールド・モデルの計 算は多くのメモリアクセスを伴うため、効 率よく計 算するためにはグローバルメモリへのアクセス
回 数を低 減することが重 要である。GPUコンピューティングでは、ス
レッドおよびブロックをどのように実 際の計 算に割り当てるかが高 速

ここで、

と

はそれぞれ固相、液相の拡散係 数を表す。

、 はを満たしている。

化のキーポイントである 。1 つの GPU が担 当する計 算 領 域サイズをnx × ny × nz とし 、図 3 のように 64 × 4 × 32 の複 数のピースに分割する。GPUで実 行されるカーネル関 数を 1 つのスレッドブロック

3

単一 GPU への実装

が１つのピースを担当するよう設定し、各ブロックで 64 × 4 × 1 のス
レッドを実行させる。ある点（i, j, k ）を割り当てられたスレッドは、z 方
向へマーチング（ ループカウンタk ）しながら（ i, j, k0）から（i, j, k0+31）の 32 格子点を計算する。ここでk0 は 32 の倍数で 0 、
32 、
64…、nz - 32である。あるスレッドがk + 1 番目の面を計算するのに必要なデータ

3 次 元 直 交 格 子 上で（ 1）式 と（ 2）式を 2 次 精 度 有 限 差 分 法で離 散の一部は k 番目の面を計算した際にそのスレッド自身によってすでに化し 1 次 精 度の時 間 積 分（オイラー法）を行う。CUDAを用いてGPU

使われている。このようなデータはレジスタに一 時 的に保 持し再 利コンピューティングのプログラミングを行い 、フェーズフィールド変数、

用することで、再びグローバルメモリへアクセスすることを回 避する

濃 度 変 数は全てビデオ・メモリ（CUDAではグローバルメモリと呼ことができる。本 研 究で用いるTSUBAME 2.0 に搭 載されたM2050

ばれる）上に確保し、計算の毎時間ステップでのPCI-Express Bus を

GPU はFermi コアのGPU であり、L1/L2キャッシュが利用できるため、介したGPUとCPU 間のデータ通 信を排 除している。3 次 元 空 間の格

これまでのような複 数のスレッド間でデータを共 有する目的（ソフト

子点（ i, j, k）上のフェーズフィールド変数 （図 2（a））と濃度変数

ウェア・マネージド・キャッシュ ）でシェアードメモリを使 用する必（図 2（b））を時 間 発 展するために必 要な近 傍 格 子 点 上の変 数の空

要がなくなった。

間配置を示す。1タイムステップでは 、それぞれの格子点でのびおよを計算するために、隣接 19 格子点のフェーズフィールド変数と隣接 7 格子点の濃度変数をメモリから読み込み、計算結果をメ

モリへ 2 回書き込む。

図 3 CUDA スレッドとブロックの割り当て

複数ノードに搭載されたGPU による計算（ a） フェーズフィールド変数のステンシル4

大 規 模問題を高速に計算するには、複 数のGPUを用いて計算することが必 須となる。GPU 内 部での並列 計算に加えて、上位でのGPU 単
位での並列化が必 要となる。大 規 模 並列計算では 、3 次 元領域分割を行うことにより領 域 間のデータ通 信量を最小にすることができるが、GPU 計算の場 合はx 方 向で分割した場 合の y-z 断 面のメモリア
ドレスが不連 続となり、メモリアクセス性 能の低 下のオーバーヘッドの方が通 信量の低 減より大きくなる。ここでは、y方向 、z 方向に分割する 2 次 元の領 域 分 割 法を用い 、それぞれの分 割 領 域の計 算を各GPUへ割り当てる。複 数 CPU 計算と同 様に隣 接するGPU 間での境
界領域のデータ交 換が必 要になるが、現 状ではGPUは他のGPUのグ（ b） 濃度変数のステンシル
図 2 時間発展に必要な隣接格子点アクセス

ローバルメモリ上のデータに直 接アクセスすることができない。そこ
で GPU 間のデータ転 送はホスト側のメモリを経 由し、次の 3 段 階で構 成される。
（ 1）CUDAランタイムライブラリによるGPUからCPUへ

04

ACM ゴードンベル賞・特別賞 ＆ SC’11 テクニカル・ペーパー

TSUBAME 2.0スパコンにおける樹枝状凝固成長のフェーズフィールド法を用いたペタスケール・シミュレーションの転 送 、（ 2）MPIライブラリによるCPU 間のデータ転 送 （、 3）CUDA

究 [2] では境界も中心領 域もGPUで計算したが、ここでは境界領 域をランタイムライブラリによるCPUからGPUへの転 送を行う。

CPUで計算しGPU 関数を分割することなく通信を隠蔽する 。ただし、

従 来の複 数ノードの CPU 計 算の場 合と比 較すると、GPUの演 算

全計算領域サイズによってCPUによる通 信と境界領域の計算にかか性 能が高いため GPU 内での計算 時 間が短い。GPU 間のデータ通 信

る時間が中心領 域のGPU 計算 時間よりも長いことがあり、CPUがボの時間が無 視できず、大 規 模 計算では大きなオーバーヘッドになる。

トルネックとなることが考えられる。

GPU 間のデータ通 信の時 間をGPU 計算と如 何にオーバーラップす（c ）Hybrid -Y metho d：（ b）と同じようにCPU 計 算も利 用するが、ることで隠ぺいするかが非常に重 要となる。ここでは、（ a）GPU-only

z 方向の境界領 域（ x-y 断面）はカーネル関 数を分割してGPUで計算

method （、 b）Hybrid-YZ method 、（ c）Hybrid-Y methodの 3 つのする。境界領域を担当するGPUカーネル関 数は 、グローバルメモリの実装を行う。

連 続アクセスであるため計算 効 率がよく、また通 信バッファへデータを詰め替える必 要もない。

（a）GPU-only method：比較のために計算と通信のオーバーラップ

TSUBAME 2.0 の各ノードには 3 GPUと 2 CPU sockets (12 CPUを行わず、
上記の通信の 3ステップを順に行う。

co res) が搭 載されている。このため CPU による境 界領 域の計 算に（ b）Hybrid-YZ method： 1つのGPUが担当するサブ領域を、y方向

1 GPUあたり 4 CPU cores を割り当て、OpenMP を用いた並 列 計算の境界領 域 、z 方向の境界領 域 、残りの中心領 域に分割する。y 、z 方をおこなっている。

向の境界領域をCPUで計算し、中心領域をGPUで計算する。先行研

図 4 GPU-only method のダイアグラム

図 5 Hybrid-YZ method のダイアグラム

図 6 Hybrid-Y method のダイアグラム05

SC’11 特集号

大規模 GPU 計算の実行性能

5

GPUコードの実 行 性 能の評 価において、GPUコードの浮 動 小 数
点 演 算 数を直 接 測 定することは困 難である。ここでは全く同じ結果を出すCPUコードに対してPAPI（ Performance API）を用いて実測した浮 動 小 数 点 演 算 数を元にGPUの実 行 時 間から実 行 性 能を評価している。図 8 の強 スケーリング（単精度計算）は、問題を固東 京 工 業 大 学・ 学 術 国 際 情 報 センタ ー の GPUスパコン で あ る

定してGPU数を増やして行くときの実行性能の向上を示している。

TSUBAME 2.0 の複数 GPU（M2050）を用い（a）GPU-only method、（b）

前述の（a）、（ b）、（ c）の 3 つの実装の違いによる特徴が表れている。

Hybrid-YZ method （、c）Hybrid-Y methodの 3 つの実装の実行性

全体の計算領域を 5123 、10243 、20483 とした。図 5 、図 6 のように能について述べる。ここでは、Al-Si合金の凝固成長の計算を行った。

計算と通信のオーバーラップを導入した（b）Hybrid-YZ method、
（c）

図7
（a）
は、大阪大学の安田秀幸教授らが大型放射光施設（SPring-8）

Hybrid-Y methodではGPU 数が少ないときには通信を隠蔽できていを用いて撮影した合金の凝固過程の画像である 。図 7（b）は、本研

る。GPU 数の増加とともにGPUの計算時間が短くなるため 、ある段究のフェーズフィールド法により、同じような 2 次 元 的な形 状で、階で計算時間より通信時間の方が長くなり、もはや通信を隠ぺいす

4096 ×128 × 4096 格 子を使ってシミュレーションを行った結 果

ることができなくなる。
（ c）Hybrid-Y methodでは、GPU 数を増やしである 。合 金 系は異なるが 、成 長 過 程が非 常によく一 致している

た時にCPUの計算時間が隠ぺいの足を引っ張る部分が大幅に改善ことが分かる。

され、広い範囲のGPU 数で最適化を導入していない（a）GPU - only
methodと比較して実行性能が大幅に改善されている。

図 7（a）  SPring-8 を用いて撮 影された合 金の凝 固 過 程の実 験 画 像（安田教 授のご好意による）
図 8    複 数 GPU 計算の実 行 性 能（強スケーリング：単 精度 計算）
弱スケーリングは 1GPUあたりで実行する問題規模を一定にし、GPU 数の増 加とともに問 題サイズも大きくする性 能 評 価である。
この 測 定 で は、1GPU あ た りの 計 算 格 子 サイ ズ を 4096 × 160 ×
128 として、単精度で計算した。図 9 のように非常に良い弱スケー
リングの結 果が得られ 、（ c）Hybrid-Y methodでは、4000 GPUを利 用して 4096 × 6480 × 13000 格 子に対して行った計 算で、2.0
PFLOPS（ GPU:1.975 PFLOPS, CPU: 0.025 PFLOPS）という極めて高い実行性能を達成した。これまでにステンシル計算の実アプリケー
ションとしてPFLOPS を超えた計算例は報告されていない。

図 7（b）フェーズフィールド法を用いてGPUで計算した凝 固 成長

06

ACM ゴードンベル賞・特別賞 ＆ SC’11 テクニカル・ペーパー

TSUBAME 2.0スパコンにおける樹枝状凝固成長のフェーズフィールド法を用いたペタスケール・シミュレーション

図 9 複数 GPU 計算の弱スケーリング（単精度計算）
TSUBAME 2.0 では消 費 電 力を詳 細にモニタリングしている 。
4000 GPUと 16000 CPU を使って計 算し、2.0 PFLOPS（ピーク性
能の 44.5 %）の実 行 性 能を達 成したときの使 用した計 算ノードとネットワークの消 費 電 力は約 1.36 MW であった（図 10）。電 力 性
能は 1468 MFlops/Watt であり、少ないエネルギー消 費で目 的とする計算結果が得られたことになる。

図 11 Al-Si 合 金の樹 枝 状凝 固 過 程の大 規 模
シミュレーション

図 10 TSUBAME 2.0 の消費電力モニター
樹 枝 状 凝 固の 3 次 元 成 長を調べるために、4096 × 1024× 4096
格 子での計 算を行った。初 期に床 面に固 相のシードを置き、柱 状晶が形成される際の相互干渉などによる空間選択のメカニズムを明らかにすることができる（図11）。

おわりに

6

フェーズフィールド法により二 元 合 金の樹 枝 状 凝 固 成 長のプロ
セスを GPUスパコンである TSUBAME 2.0 を用いて大 規 模 計 算を行い 、格子法に基づいたステンシル計算でありながら単精度で 2.0
PFLOPS という極めて高い実 行 性 能を達 成した。これはピーク性能に対して 44.5% であり、GPUスパコンが実 用アプリケーション
に対して十 分 有 効であることを示すことができた。同じような計算でもGPUスパコンを使うことにより大 規 模 計 算が可 能になり、さまざまな分野での研究の発展や開発の進展が期待できる。

07

SC’11 特集号

謝辞
本研究は TSUBAME グランドチャレンジ大規模計算制度を利用して実施させて頂き、一部は科学研究費補助金・基盤 研究（B）課題番
号 23360046「GPU スパコンによる気液二相流と物体の相互作用の超大規模シミュレーション」、科学技術振興機構 CREST「次世代テクノロジーのモデル化・最適化による低消費電力ハイパフォーマン
ス」および「ポストペタスケール高性能計算に資するシステムソフト
ウェア技術の創出」、日本学術振興会（JSPS）グローバル COE プロ
グラム「計算世界観の深化と展開」(CompView) から支援を頂いた。
記して謝意を表す。
参考文献[1]

R . Ko b ay ashi: M o d e lin g a n d num e r i c a l sim ul at i o ns o fdendritic crystal growth. Physica D, Nonlinear Phenomena,63(3-4), 410 - 423 (1993)[2] 小川慧 , 青木尊之 , 山中晃徳 : マルチGPUによるフェーズフィー
ルド相 転 移 計 算のスケーラビリティー ― 40GPUで 5 TFLOPSの実 効性能 , 情 報 処理学会論文 誌コンピューティングシステム
Vol. 3 No. 2 67-75 (2010 June)[3] A . Ya m a n a k a , T. A o k i , S . O g a w a , a n d T. Ta k a k i : G P U accelerated phase-field simulation of dendritic solidification
in a binar y alloy. Journal of Cr ystal Grow th, 318 (1): 4 0 45 (2 011). The 16 th International Conference on Cr ystal
Grow th (ICCG16 )/ The 14 th International Conference on
Vapor Growth and Epitaxy (ICVGE14)[4] T. Shimokawabe, T. Aoki, T. Takaki, A. Yamanaka, A. Nukada, T.
Endo, N. Maruyama, and S. Matsuoka: Peta-scale Phase-Field
Simulation for Dendritic Solidification on the TSUBAME
2 . 0 Supercomputer, in Proceedings of the

2 011 ACM/I E E E I n t e r n a t i o n a l Co n f e r e n c e f o r H i g h Pe r f o r m a n c e
Computing, Networking, Storage and Analysis, SC’11, IEEE
Computer Society, Seattle, WA, USA, Nov. 2011.

08

TSUBAME2 における大規模生体流体力学シミュレーション
Massimo Bernaschi* Mauro Bisson* 遠藤 敏夫 **Massimiliano Fatica*** 松岡 聡 ** Simone Melchionna** CNR-IAC, Istituto Applicazioni Calcolo, Consiglio Nazionale delle Ricerche, Rome, Italy

Sauro Succi*** 東京工業大学学術国際情報センター*** Nvidia Corp. Santa Clara, CA, USA生 体 中の流 体のマルチスケールシミュレーションのためには、複 雑 形 状を流 動する粒子の運 動を計 算する必 要があり、数億 個オーダーの粒子の相 互作用、粒子と流 体との相 互作用を計算する必 要がある。一つの例として、人 間の冠 動 脈の血 流を赤 血 球のサイズの空 間 解 像 度を用いて、生 理 学 的なヘマトクリット値（赤 血 球の体 積 率）の条 件で計 算する。本 稿ではTSUBAME2スーパーコンピュータ上で高いスケーラビリティを持つ血行力学シミュレーションの手 法を提 案し、600TFlops 以
上の高性 能を実 現した 。この結果は、新 規 数 学モデル、計算アルゴリズム、コンピュータアーキテクチャ、チューニング技 術の統 合により、臨 床 的に意 義のある生体 内流れの予 測が可能であることを示している。

1

はじめに

TSUBAME2 の 4000GPU に均等に割り当てる必要がある。この領域
分割問題は流体計算部分だけであっても困難であるが，さらにこのマルチスケール法においては MD 部分の負荷のバランスも良好に保つ必要がある。
このような複雑かつ大規模形状のシミュレーションは既存研究で輸送現象は生体における基本現象であり、筋肉の収縮、消化現象、は稀である。本研究では、大規模並列アーキテクチャの利用によって、細胞への栄養物の輸送、血液循環はその例である。血液は代表的

過去最大規模のシミュレーションを、形状の規則性に依存せずに可な生体流体であり、代謝、免疫反応、組織修復などの基本的な生能であることを示す。

理的機能に必要な、生物学的物質を供給する。
生体内の現実的な血液と血管を計算することは困難な課題である。なぜならその計算モデルは、複雑な形状の血管の中を流れ、心拍によって不規則に流速や圧力が変化する流体の動きと、赤血球、白血球や他の粒子の挙動を統合する必要があるからである。
この数年で、大規模血行力学シミュレーションは大きな進歩を遂
げた [1-3] が、現実的な形状やサイズの血管に対して、流体と赤血球などの粒子挙動を連立するシミュレーションには至っていなかった。
流体の圧力による非局所的相関のために 、大域的な形状が、特に動脈の壁面せん断応力に大きな影響を与える 。壁面せん断応力は 、動脈硬化につながる複雑な生体力学的変化を誘発するとされている。
正確で信頼できる血行力学（壁面せん断応力）のシミュレーションは、
循環器系疾患の進行を予測するための、新しい非侵襲的な手法となることが期待される。

図 1 シミュレートされる冠 動脈の形 状、Lattice Boltzmannメッシュの中に赤血球（RBC）
が存在する。

本稿では、CT 血管造影から再構築された人間の冠動脈の 、初のマルチスケールシミュレーションについて述べる。冠動脈は心筋に血液を供給し心臓全体に渡るネットワークを成す。シミュレーション
規模は 5cm 程度であり、解像度は赤血球の直径の約 8 μ m に対しマルチスケール生体流体力学

10 μ m 程度である。シミュレーション対象は図 1 に示すような 、約2

10 億ノードから成る流 体と 1000 万～ 4 億 5000 万の粒子である。
シミュレーションは、流体のための Lattice Boltzmann（LB）法と粒生物は主に二つの要素から成り立つ：水溶性の溶媒、血漿 、サイト

子のための特化された Molecular Dynacmics（MD）法を統合した、ソルなどの液体と、溶媒の中を流れる生体粒子、細胞 、タンパク質、[4,5]

MUPHY ソフトウェア（Multi PHYsics/multiscale）に より行 う。

このシミュレーションは TSUBAME2 スーパーコンピュータの 4000
GPU を用いて 500TFlops 以上の平均性能、90%以上の並列化効率を達成した。

09

DNA などの懸濁物質である。このような条件をシミュレートするた
めに、MUPHY [4,5] ソフトウェアを開発してきた。
以前のバージョンの MUPHY と比べ 、現在のフレームワークでは
流 体と粒子の挙 動の数値的解法を扱うことができる。これにより

本 研 究 は、 高性 能 計算 技 術と物 理 / 計算モデリングの双 方の

MUPHY は 、溶媒と溶質のライブラリを基にして、生体内流れ向けの側 面 に お いて 特 徴 を持 つ。 非 常 に 複 雑 な不 規 則 形 状 の 領 域 を計算基盤として利用することができる。機能の一部は以下の通りで

SC’11 特集号

ある：ニュートン / 非ニュートンレオロジー特性の選択、分子流体力点が生じる。実際、LB 法の計算コストは格子点数 M について線形学から確率的摩擦動力学の再現のためのカーネルの選択、確率的変である。ここで O（M）=O（N/c）であり、c = N/M は溶質の密度の定動の効果 、分離した粒子または重合体と分子の形状の溶質の選択、数である。粒子間の直接力の計算のために link-cell 法を用いること

多様な流体力学形状の懸濁物質 、異なる粒子間力の選択 、規模に対により、粒子力学の計算量は O（N）となる。この性質のために、溶媒応する粒子・流体間結合機構 、不規則境界の物質と組織への対応。

を介した粒子間作用は局所的で陽的となり、LB-MD 結合のコストはこのような選択肢により、生体流体は様々な手法でモデル化でき、格子数と懸濁物質数に比例する。しかしながらこの計算は 、O（100）

生物学的・生理学的システムのマルチスケール・マルチフィジックス

の格子点と不連続なメモリアクセスの高コストのために、最も時間をの挙動が理解できる。

要する部分である。

MUPHY は 2 つの計算エンジンを持つ。一つは流体運動学上の流体の挙動を Lattice Boltzmann（LB）法により扱う：衝突による手法
により連続体中の流体の挙動を再現する（Navier-Stokes 方程式に基づく直接的巨視的な流体の挙動と対象的に）。二つ目のエンジン

MUPHY

は Molecular Dynamics（MD）と技術的側面が共通する手法に基づ

3

き、Lagrangian 体の挙動を扱う。しかしながら懸濁物質の性質のためには通常の MD では対応できず、本質的な拡張が必要である。

MUPHY（MUlti PHYsics/multiscale）ソフトウェアは元々 Fortran 90

最後に、流体と粒子の結合のために運動モデルに基づいて設計さで記述され、並列化のために MPI を用いていた [4]。不規則形状を柔れたカーネルが存在する。これは巨視的流体力学に基づく stresslet

軟でかつ効率的に扱うために、MUPHY は間接アドレッシング法をや境界積分法などと大きく異なる。これらの計算環境により、現実用いる 。これによりメモリ必要量を最小に抑えるだけでなく、良好な的な条件における生体流体の研究に利点をもたらす、時間に陽なシ

負荷分散と、プラットフォームに応じた最良な通信パターンの選択がミュレーション技術を提供する。特に 、1）冠動脈のような不規則な可能となる。当初 MUPHY は 、比較的低クロックの数千の PowerPC

境界形状に対応する 、2）粒子の動力学を用いるため、複雑な流動プロセッサを高速専用ネットワークで結合した IBM BlueGene アー

学を描写できる 、3）メッシュ以下の空間解像度の流体学上の相互作

キテクチャ [7] 上で開発された。BlueGene/P 上で、我々が利用可能な用を可能とする 、4）流体と粒子間の複雑な界面の記述を避ける。

最大構成（294,912 コア）で良好なスケーラビリティをすでに示して後の二つはシミュレーション性能の加速と、特に物体が流体に及ぼいる [11]。この時数十 TFlops の性能を達成したが、PowerPC アーキ

す剛体力に関連した数値安定性の保証のために必須である。フレー

テクチャの SIMD 的命令を利用できないなどの制限があった（これらムワークは、剛体力についてシングルもしくはマルチ時間ステップアの命令は連続アクセスを必要とするが、MUPHY ではデータアクセス

ルゴリズムを採用する。ここで強調したいことは、LB 法が、力の急激パターンが不連続である）。一方で近年の GPU の計算性能の向上をな変化に耐性があり、1000 のオーダーのレイノルズ数のほふく流に受け、我々は GPU クラスタを対象とした MUPHY を開発した [11]。

よる濃厚サスペンジョンをシミュレート可能なことである。この利点

複数 GPU を用いるソフトウェアは、GPU 内と GPU 間の、二レベにより現実的な生物学的条件下において広範な物理現象に対応し、ルの並列性の記述を必要とする。通常のマルチコアシステム上では構造的関係に支配される連続体に同化できない、非局所的な流動学

ハイブリッド方式（OpenMP+MPI など）もしくは単純な分散メモリ方

的反応の再現が可能となる。

式（MPI ライブラリが共有メモリでも効率的に動作することを期待し

2.1 Lattice Boltzmann 法と Molecular Dynamics 法て）にて実装する一方、複数 GPU システムにおいてはハイブリッド
方式が唯一の選択肢である。考慮すべき点はそれだけではない：通LB 法 [10] は、離散速度で動く「流体粒子」の格子点 x、時間 t にお常のマルチコア環境では並列度は数スレッド、ハイエンドのものでもける存在確率の singlet 分散の発展に基づく。流体粒子は物理的粒数十スレッドである。GPU ではハードウェアを十分に稼働させるため子の一団 (population と呼ばれる ) の集団的動 作を表現する。流には数百スレッドが必要であり、はるかに細粒度の並列性が必要で体と個体の結合は回転並進カーネルに基づき、それらは流体中を動ある。さらに、GPU 上のグローバルメモリのデータを同マシンの別く剛体を、貫通できない物質、柔軟な小胞、またはその組み合わせ

GPU とやりとりする際には、CPU を介する必要がある（ただし最新として扱う。

の NVIDIA GPU と CUDAドライバにより、条件によっては可能となる）。

流 体 と固 体 の 流 体 力 学 的 作 用 は、i 番 目の 粒 子 を 中 心 とした 伝

GPU 間の通信のために CPU を介することによるオーバヘッドは

達 関 数に従 い 形 成される。 関 数は 球もしくは 楕 円 体 対 称 であり発生するものの、CUDA の stream の概念と非同期メモリコピーによcompact support である。流体力学的形状はメッシュサイズより小り、GPU と CPU 間のデータ転送とカーネル（GPU 上で実行される関

さくなりうる [6]。流体と粒子の結合は N 個の懸濁物質ごとに格子点数 ) の実行をオーバラップすることが可能である）。さらに CPU 上の{x} にまたがる畳み込み演算を必要とする。

関数（MPI 関数など）の実行は GPU とは並行に行われる。結果として、基本的に流体力学的作用の計算量は系のサイズの三乗で増加す

CPU は MPI 用のコプロセッサのような役割を果たす。

るが、流体と固体が並行で時間発展するため、アルゴリズム上の利

10ACM ゴードンベル賞・奨励賞

TSUBAME2 における大規模生体流体力学シミュレーションも課題となったのは、部分領域の境界部分の隣に位置し、領域内粒

3.1 領域分割

子とも領域間粒子とも作用するような粒子の特定である。我々はセ本研究のシミュレーション対象の領域は、図 1 に示すように非常にルの概念 [8] によりこれを解決する。ここでのセルは、作用のカットオ

不規則であるため、計算資源間で領域を分割すること自身が大きフ距離以上の辺の長さを持つ立方体であり、計算される不規則形状

な課題となる。以前の実験では、計算負荷を均等に分散させるためを埋め尽くすものである（図 2 ）。この概念により、プロセッサが領域[12]の並列版である。PT-

内 / 領域間粒子ペアを効率的に探索可能であり、領域境界に存在す

SCOTCH を用いた。このツールはグラフバイセクションアルゴリズム

る粒子と領域間を移動する粒子の上位集合を交換することにより通に基づき、計算領域の形状にかかわらず分割を行うものである。し信量を削減可能である。

に、SCOTCH グラフ / メッシュ分割ツール

かしながら形状の知識が無いと、分割数が増え、部分領域がより不

MUPHY のもう一つの構成要素は流体 - 粒子間相互作用の計算部規則になると、分割の性質が悪くなると分かった。境界面積が増え分である。各懸濁粒子は、流体力学的な力と、流体の巨視的な速通信オーバヘッドが増加する。より適した解は、グラフに基づく分度と、4 × 4 × 4 サイズの格子に渡る渦度の影響を受ける。同様に

割と flooding に基づく手法（graph-growing 法とも呼ばれる）を、

格子点は、粒子の運動量輸送の影響を受ける。これらの非局所的な

以下のように統合した場合に得られると分かった：格子はまず PT-

処理は、ある粒子を持つ各プロセッサが流体力学的数量を周りの領

SCOTCH により固定数（256 ）の部分領域に分割される。そして各部域と交換するために、複数の通信ステップを必要とする。

分領域は flooding アルゴリズムによりさらに分割される。

流体 - 粒子間結合は、粒子と隣接領域の境界セル中の格子ノード

MUPHY では通信パターンは実行時に判明する。各タスクは、そのの間の作用を含む非局所的な処理である。このため結合計算においタスクがシミュレーション中にアクセスする必要のある格子点を持っても領域間力の計算と同様 、セルの概念を用いる。流体が粒子に及た、隣接タスクを決定する。シミュレーション中に、LB アルゴリズ

ぼす力やトルクの計算のために、各プロセッサは境界セル中の粒子

ムにおいては非局所的な population の流れが存在し、分子のシミの情報を近隣と交換する。最後に外部の粒子に関連する力やトルクュレーションにおいては粒子の移送やドメイン間力の計算があるこは隣接プロセッサと交換され、受信側では外部からの影響が境界粒とに注意が必要である。前処理段階においては主に MPI 集団通信

子に伝えられる。この手法は他の手法より効率的である。たとえばが用いられ、それ以降の実行においてはほとんどの通信が以下のよ

粒子の代わりに境界セル中の格子点を交換した場合、一ステップにうな形の一対一通信である：ノンブロッキング通信関数が用いられ、おける通信回数を削減することができる。しかしこの手法では界面受信処理が先に起動されてから送信処理が起動される。その後各に近いすべての格子ノードの交換のために通信オーバヘッドは大きく

タスクは MPI_wait 関数により、ノンブロッキング通信の終了を待つ。

なる 。

Lattice Boltzmann 部 分 にお いては 大 域 数 値（x, y, z 方 向 運 動

最後に、粒子が流体へ及ぼす運動量輸送については、境界上の粒子量など）の計算にの み MPI 集 団 通 信（reduction）が用いられる。

の情報が交換される。この点は流体から粒子への影響の場合と同様

Molecular Dynamics 部分においてもほとんどの通信が一対一でだが、二回目のデータ交換は不要である。

ある。ただし不規則的形状であることにより領域分割が例外となり、その点は次節で述べる。
3.2 並列 Molecular Dynamics
多くの並列 Molecular Dynamics アプリケーションにおいてはシミュ
レーション形状は規則的であり、各タスクがほぼ同数の粒子を持つ
ように単純な Cartesian 分割が用いられる。不規則的形状の場合にこのような方法を用いると、LB のための分割と MD のための分割という、二つの異なる分割が生じる。その結果一つの部分領域が二つ以上のプロセッサに所属することとなり、粒子と流体の相互作用が複雑な通信パターンを持つ非局所的処理となってしまう。そのため、
我々は MD と LB の間で並列分割を一致させることにした。これにより各タスクは LB の演算、MD の演算 、さらに粒子 - 流体間作用の演算をほぼ局所的に実行できる。本手法では LB 格子は 、各粒子が所属する部分領域を特定するためにも用いられる：場所 R にある粒子
は、そのベクトルにもっとも近い整数ベクトルに対応する格子を持つ

図 2 境界セルと内部セルを不規則領域に敷き詰めた様子。

部分領域に所属する。粒子は領域間にほぼ均一に存在すると期待で隣接領域とのデータ交換を限定することができる

きるため、MD 部分の負荷分散も良好に保たれる。
不規則形状の領域に適した新規並列化戦略を開発した。なかで11

SC’11 特集号

実験結果

4

通信に費やしている 。この結果は非同期通信方式が良好に動作し、GPU が計算にほとんどの時間を費やす一方で、CPU がその MPI コプロ
セッサとして働いていることを示している。
TSUBAME2 の 1334 ノード（4000MPI タスク）を用いた場合の、（重みつき）平均性能は 、600TFlops 弱となる。図 5 にその内訳を示す。ここれまでに述べたシミュレーションについて、全体実行時間と、計算の成果を実用にもたらすとすると、この TSUBAME2 全体を用いることと通信の内訳について測定する。シミュレーション対象は以下の通により、完全な心拍の 、マイクロ秒の時間解像度における、赤血球をりである：流体については約 10 億の格子ノードから成る領域であり、全て含んだシミュレーションを、48 時間で行うことが可能である。

約 3000 億ノードの bounding box に収まる。粒子は 4 億 5000 万の赤血球（RBC）である 。流体 - 粒子結合のパラメータについては、
以下のように設定した：各赤血球は小球の最小・最大主要方向について 4 μ m および 8 μ m の流体力学的形状を持つとし、ヘマトクリッ
トは 58%とした。
多くの演算は単精度浮動小数で行われ、一部の reduction 処理
は倍精度で行われた。LB 計算においては各 GPU スレッドは複数格
子ノードの更新処理を行い、そのノード数は利用した GPU 数に依存する（スレッド設定は各 GPU 固定とした）。たとえば 512GPU の場
合各スレッドは 8 格子ノードを担当する。MD 計算においては、粒子
間作用の演算は粒子ごとに行われる。ここでスレッドグリッドは粒子の配列にそのまま対応するスレッドはそのグローバル ID に応じて粒子に対応し、作用ペアの探索は各スレッド独立に行われる。各スレッドは隣接セルを調べ、各作用ペアについて力への寄与を計算する。
この実験により、生理学的なレベルのヘマトクリットにおけるコードの信頼性の基礎試験を行う。
図 3（上部）では、シミュレーションステップごとの実行時間と、LB
部分、MD 部分の内訳時間を示す。いくつかの知見を述べる：MUPHYの Lattice Boltzmann 部分の 1GPU 上の性能は、高度に最適化され

図3 上部: 時間ステップあたりの経過時間。下部: 並列化効率。

た他の CUDA LB カーネルと一致する。GPU 数を増やした場合には実
行時間は大幅に削減されており、4000GPU の場合には 256GPU の12.5 倍の性能向上が見られ、並列化効率は約 80%に相当する。MD
における直接力の計算は効率 95%に相当する。LB 部分は非常に効率的であり、負荷は常に 4%以下である。
またこの結果は非同期通信と通信・計算のオーバラップの組み合わせにより GPU 間の直接データ交換ができない問題を緩和できていることを示す。
図 3（下部）は全体の並列化効率を示す。その基準は、利用可能なメモリの範囲でシミュレート可能なヘマトクリットレベルの場合には 256GPU であり、それより高レベルの場合は 512GPU である。
1024GPU までは線形を超えた性能向上が起こっており、また利用した最大 GPU 数の場合の並列化効率は 80%程度であると分かる。この高いスケーラビリティは、複数の枝を持つ冠動脈向けに提案したグラフ分割と flooding 法のハイブリッド方式により実現できたと考えられる。1200 領域を用いるテストケースでは赤血球の分布には広い分散があったが、MD/LB 部分の実行時間は小さかった。これは部分領
域間の境界が小さく計算と通信の共有が最適化されたためである。
図 4 は 4000GPU を用いた場合の通信時間の分散を示す。平均的

図 4 4000の各タスクが通信に費やした時間の割合の分布にはほとんどのタスク（4000 中約 3000）が、全体時間の 50%程度を

12

ACM ゴードンベル賞・奨励賞

TSUBAME2 における大規模生体流体力学シミュレーション

謝辞
本研究は TSUBAME グランドチャレンジ大規模計算制度で実施された。
E. Kaxiras, C.L. Feldman, A.U. Coskun, F.J. Rybicki, A.G. Whitmore,
G. Amati, F. Pozzati, F. Schifano の各氏との議論に感謝する。
参考文献
[1]

D.A. Vorp, D.A. Steinman, C.R. Ethier, Comput. Sci. Eng., pp.
51 (2001).

[2] A . Q u a r t e r o n i , A . Ve n e z i a n i , P. Z u n i n o , S I A M J . N u m .
Analysis, 39, 1488 (2002).
[3] L. Grinberg, T. Anor, E. Cheever, et al., Phil. Trans. Royal Soc. A,
367 1896 2371 (2009).
[4] M. Bernaschi, S. Melchionna, S. Succi et al., Comp.Phys.
Comm., 180, 1495, (2009).
[5] S. Melchionna, M . B ernaschi, S. Succi et al, Comp. Phys.
Comm., 181, 462, (2010).
[6] S. M el chio nna, Ma cro m o l. T h e o r y Sim ., D O I: 1 0 .1 0 0 2 /mats.201100012 (2011).
図5 各計算コンポーネントの計算時間の割合とそれぞれの性能[7] M . B e r nas chi, M . Fat i c a , S . M e l chi o nna , S . Su cci a n d E .
K a x ir a s , Co n c u r re n c y a n d Co m p u t at i o n: Pr a c t i ce a n dExperience, DOI: 10.1002/cpe.1466 (2009).

おわりに

5

[8] M. Bisson, M.Bernaschi, S.Melchionna, Commun. Comput.
Phys., 10, 1071 (2011).
[9] S. Melchionna, J. Comput. Phys. 230, 3966 (2011).
[10] S . S u c c i , T h e L a t t i c e B o l t z m a n n E q u a t i o n f o r F l u i d
Dynamics and Beyond, Oxford University Press, USA (2001).

世界トップクラスの GPU クラスタを用いた、心臓規模の冠動脈の計
算生体流体力学シミュレーションを初めて行った 。シミュレートされた血行力学系は赤血球レベルの空間解像度を持ち、人間の冠動
脈の複雑な形状から成る 。計算された環境は 10億の流体ノードと、並行に流れる数億の赤血球を含む 。TSUBAME2 スーパーコンピュータの持つ 4000GPU を用いて、約 600TFlops の性能と 90%以上の並列化効率を実現し、4 億 5000 万の赤血球の存在下で約 2 兆回 /秒の lattice 更新を行った。
この成果は、高性能計算技術と物理的 / 計算モデルの双方における新規開発によるものであり、我々の知る限り、理想化されない形
状に対応した初の実装である。この研究は、現実の生体流体力学の研究に対するコンピュータシミュレーションの可能性を大幅に進歩させるものであり、心臓血管の臨床における応用を可能にするもの
である。

13

[11] A. Peters et al., Proceedings of Supercomputing 2010, New
Orleans, 2010.
[12] http://www.labri.fr/perso/pelegrin/scotch

大規模グラフ処理ベンチマークGraph500 のTSUBAME 2.0 における挑戦
鈴村 豊太郎 */**

上野 晃司 *

* 東京工業大学大学院・情報理工学研究科**IBM Research – Tokyo

Graph 500とは、スーパーコンピュータのグラフ処理性能を測定する新しいベンチマークである。
スパコンのベンチマークでは、数値 計算性 能を測るLinpackによるTop 500 が有名だが、近 年、大 規 模グラフ処 理が、重 要性を増しており、Graph500ベンチマークが広がりを見せている。Graph500 のリファレンス実 装は、使 用されているアルゴリズムの問題により、分 散メモリ環 境で大 規模にスケールさせることができなかった。そこで、大 規模にスケール可能な２次 元 分割に注目した。本論文では、２次 元分割をTSUBAME2.0 上に実装し、1366ノードで頂 点数 2^36
（ 68.7 billion）、エッジ数 2^40（ 1.1
trillion）のグラフ（Graph500 のScale 36）のBFS( 幅優 先探索 )を10.955 秒で計算した。
TEPS 値は100.366 GE/sであり、2011年 11月に発表されたランキングでは世界 3 位を獲得した。

はじめに

1

大 規 模グラフ処 理はWeb ページのリンク解 析、タンパク質 間の相互作用解 析、サイバーセキュリティ、VLSIのレイアウトや道 路 網、送電 網の最 適 化など様々な応 用 分 野あり、近 年 盛んに研 究されている。従 来、スーパーコンピュータは物 理シミュレーションなどの数値計算に、主に使われてきたが、大 規 模グラフ処理も重 要なアプリケー
ションとなりつつあり、そのような中、スパコンのグラフ処 理 性 能を測る、Graph500 [1] という新しいベンチマーク登場し、注目を集めている。Graph500 は、スパコンの通 信 性 能や、グラフデータを格 納
するメモリの大きさや、メモリへのランダムアクセス性 能を測るという、データインテンシブなベンチマークであり、数値 計算性 能を測る
Top 500ベンチマークとは計 測する性 能が全く異なる。本 論文では、
Graph500 の概 要、我々が提 案するスケーラブルな最 適 化 実 装とTSUBAME2.0 における性 能評 価について述べる。

Graph500 の概要

2

本 章では，Graph500ベンチマークの概 要と、分 散 BFSアルゴリズムについて述べる。
２-1

Graph500ベンチマークの概 要

図 1 クロネッカーグラフ[4]Graph500 は、大 規 模なグラフに対してBFSによる探 索を実 行する

ベンチマークを実 行するプログラムは、（ a）
グラフデータの生 成、（ b）

ベンチマークである。単位時間に処理できたエッジ数と、扱える最 大計算するのに最 適なデータ構 造への変 換 （、c）
BFSによる探 索、（ d）

問題サイズが評価指標となる。計算インテンシブなTop500ベンチマー

計算 結果の検 証の４つの部分から成る。ベンチマークの実 行順は次

クと違い、Graph500ベンチマークは、
データインテンシブなベンチマーのようになっている、最初に
（a）
（
、b）
によりグラフデータを構築し、グクである。扱える問題サイズは、グラフの頂 点数 =2SCALEであるよう

ラフから始 点を 64 個 選ぶ 。次に、64 個の始 点それぞれに対して順なSCALE値で表す。
単位時間に処理できたエッジ数は、TEPS (Traversed

番に、
（ c）
BFSによる探 索と、（ d）
計算 結果の検 証を行う。複 数の始点

Edges Per Second) 値で表す。例えば、100 万 TEPSとは、100 万個のからの探 索を同 時に行うことはできない。時 間を計 測しベンチマー

枝を持つ連結グラフのBFSが 1秒で完了した場合の性能である。

クとする部分は、（ b）のグラフデータ構 造への変換（ Kernel 1）と、
（ c）

14

Graph500 ランキング 3位

大規模グラフ処理ベンチマークGraph500 のTSUBAME 2.0 における挑戦のBFSによる探 索（Kernel 2）のみである。
（a）では、枝 数が頂 点 数の点γが入っていたとすると、11、12 行目でγの隣接頂 点が訪問済みかど

16 倍となるようなクロネッカーグラフ[4] を生 成する。枝はすべて重みうかチェックされ、まだ訪問していない頂点はNQに格納される。次のなし、無向 辺である。ここで生 成されるデータは規 則 性のない順 番

レベルでCQにはこれらの頂 点が格納されていることになる。9 行目ので並んだ、枝のリストである。
（b）
では
（a）
で生 成された枝リストから、

forと、11行目のforは並列化が可能なループである 。リファレンス実装

隣 接 行列のCSR (Compressed Sparse Row)や、CSC (Compressedの４つのMPI 実装は、基本的にはLevel-synchronized BFSを実装してSparseColumn)などのグラフデータ構 造に変 換する。
（ c）のBFSは、

いるが、グラフデータの分 散 方法などに違いがある。４つのリファレン

BFSで辿った頂 点の軌 跡であるBFS 木を出力する。
（d）では、このBFS

スMPI 実装の処理の仕方の違いや TSUBAME 2.0 上における性能特性

木が正しいかどうかチェックする。このチェックでは、BFS 木にループの結果は我々の先行研究 [2] を参照して頂きたい。

がないこと、枝の張っている頂 点同士の深さの差が 1以下であること、などの 5つのルールを満たしていることをチェックする。
2-2 分 散 BFS アルゴリズム

2 次元分割によるスケーラブルな実装

Graph500 のリファレンス実 装には、OpenMPで書かれた実 装や、3

MPIで書かれた実 装、Crayの共 有メモリ型プログラミング環 境 用の
実 装、など複 数の種 類が用 意されている。TSUBAME2.0 で分 散 実

リファレンス実 装は全て 1 次 元 分 割を使っているが 、1 次 元 分 割は

行するには， MPIで書かれた実 装を使 用する。MPIで書かれた実 装スケールさせることが難しい [2] 。そこで、隣 接 行列を 2 次 元に分割すには、さらにアルゴリズムや実 装 方法の異なる 4 種 類の実 装が用意

[3]
るアルゴリズム（ 2 次 元 分割）を実装したプロセッサをP = R x Cの2されている。これらの実装は、対 象としているプログラミング環 境や次 元メッシュ（ mesh）に配 置する。このメッシュの行を「プロセッサ

分 散 方法などは異なるが、全てベースとなるアルゴリズムとしてLevel-

行」、列を「プロセッサ列」と呼ぶことにする。隣 接 行 列を図 4 のよsynchronized BFSを使っている。このアルゴリズムは、各レベル（深うにR*C 個の行とC 個の列に分割し 、
プロセッサ（i, j）
は、隣 接 行列のさ）について、そのレベルの頂 点をすべて処理してから、次のレベルに

A _(i,j)^((1) )~A _(i,j)^((C) )のCブロックを担当する。頂 点は、R x C

進むというアルゴリズムである。

個のブロックに分割し、
プロセッサ（I, j）は、j*R+I 番目のブロックを担当する。
１レベルにつき、expandとfoldの２段階の通 信を行う。各プ
ロセッサは自分の担当する頂 点ブロックのCQを同じプロセッサ列の他のプロセッサに送 信する。これをExpandという。Expandは１次
元 分割の縦分割と同じように、
CQをコピーする通 信であるが、隣 接 行
列は横にもC 個に分割されているので、通 信は、同じプロセッサ列の他のプロセッサとだけ行う。次に、各プロセッサはCQと各プロセッサが持っている部分 隣 接 行列から、
CQの隣 接頂 点を探す。PRED やNQを更新するため、CQの隣 接頂 点を、その頂 点の担当プロセッサに送信する。この通 信をFoldという。PREDを更新するのに、親の頂 点が必 要なので、Foldでは、CQの隣 接頂 点と、親 頂 点（CQの頂 点）の組
みを送 信することになる。Foldは１次 元 分割の横 分割と同じように、
CQの隣 接頂 点を担当プロセッサに送 信する通 信だ。しかし、２次 元
分割では、隣 接 行列の分割方法から、Foldの通 信を行う必 要のある
相手は、同じプロセッサ行の他のプロセッサのみとなる。
２次元分割の利点は、通信で絡むプロセッサ数が少ないことである。

アルゴリズム１はLevel-Synchronized BFS の擬 似コードである。

１次 元 分割では、２種 類の分割方法のどちらも、全 対全の通 信が必

ま ず、BFS 木 を 格 納 するPREDと、頂 点 が 訪 問 済 み かどうか を 格要だったのに対し、２次 元 分割の場 合、Expandでは同じ列のノード

納 するVISITEDを 初 期 化 する。PRED [v] は 頂 点 vのBFS 木 に お ける

(R-1)プロセッサと、foldでは同じ行のノード(C-1)プロセッサとしか通

親 頂 点 を 表 す、初 期 値 -1はBFS 木 にまだ 入って い な いことを 表

信を行わない。よって、通信するプロセッサ数を少なくすることができ、して い る。VISITED [v] は 頂 点 v が 訪 問 済 み かどうか を 表 す。初 期

大 規 模に分 散可能になる。

値 0 はまだ 訪 問して い な いことを 表して い る 。次 に、BFSの 始 点と な る 頂 点 を CQ（Current Queue）に 入 れ 、探 索 を 開 始 す る。
探索においては 、7~16 行の１ループが１レベルに相当する。このルー
プ中で、CQは現在のレベルで訪問する頂点 、NQ（Next Queue）
は次のレベルで訪問する頂点が格納されている。例えば、レベル1でCQに頂15

SC’11 特集号

装は、
リファレンス実装のreplicatedと比べると、性能が低い。これはreplicatedのアルゴリズムはノード数が小さい場 合には通 信データ量を小さくすることができ、有 利だからである。図 3 から分かるように、
replicatedの優 位性もノード数が増えるにしたがって急 激に低下し、通 信データ量は 512ノードで２次 元 分割と逆 転する。実 際、図 3 から
replicated-cscはノード数 128 で既に性能の限界が見え始めている。
また、図 4 は 1ノードあたりの Weak Scaling によるスケーラビリティの評 価だが、1024ノードまでノード数を増 加させても性 能が向 上し、十 分なスケーラビリティが得られていることがわかる。

図 2 隣 接 行 列の２次 元 分割

性能評価

4
図 3 ２次 元 分割と参 照 実 装の比 較

TSUBAME2.0 上での性能評 価の結果について述べる。TSUBAME2.0
は、1400 以 上 の ノード が Fat-Tree に よ るフルバイセクシ ョン のInfinibandネットワークで接 続されている。各ノードには、Intel CPU
Xeon 5670 2.93GHz（Westmere EP、6コア、256-KB L2 キャッシュ、
12-MB L3）が２つ、NVIDIA M2050 GPU（ Fermi）が３つ、48GBのメモリが搭載されている。通信は、各ノードはInfiniband QDRが２リン
ク使 用可能で、合計 80Gbpsの通信バンド幅を備えている。
最 大 1024ノードまで使 用して実 験した。なお、TSUBAME2.0 はGPUメインのスパコンだが。GPUは使 用していない、TSUBAME2.0
は 1ノードあたり物理コア12 個だが、SMTを有効にすると仮 想的に 24コアになる。1ノード24コアとして、各プロセスに均 等に割り振った。
gcc 4.3.4（ OpenMP 2.5）、MVAPICH2 1.6 [4] 。比 較するリファレンス
実装は、執 筆時点で最 新のversion 2.1.4 である。
図 3 は２次 元 分割とリファレンス実装の比 較である。横軸はノード

図 4 1024ノードまでのスケーラビリティ

数で縦軸はTEPS（GE/s）である。リファレンス実装のreplicated-csr 、
replicated-cscと最 適化実装は、１ノードあたり2MPIプロセスで実行
し、Simpleは１ノードあたり16MPIプロセスで実行した。図 7 は図 6 の性能をノード数で割り、１ノードあたりの性能を算出したグラフである。
リファレンス実装のsimpleを参 考に掲載した。リファレンス実装で、データがない部分はエラーなどで計測できなったところである。
２次 元 分割の実 装は、リファレンス実 装のsimpleの２倍 程 度の速度が出ている。これは、送 信処理と受 信処理の並列化や、OpenMPによるプロセス内の並列 化の効 果によるものである。２次 元 分割の実

16

Graph500 ランキング 3位

大規模グラフ処理ベンチマークGraph500 のTSUBAME 2.0 における挑戦

まとめと今後の展望

5

本 論文では大 規 模分 散環 境でスケールさせるため２次 元 分割によるBFSを実装した。2011年11月のGraph500におけるスコアは、1366ノードで頂 点 数 2^36
（68.7 billion）、エッジ数 2^40
（1.1 trillion）のグラフ（Graph500 のScale 36）のBFS（幅優 先探索）を10.955 秒で計算した。
TEPS 値は 100.366 GE/sであり、2011年 11月に発 表されたランキン
グでは世界 3 位を獲 得した。我々は、２次 元 分割の他に、通信データの圧 縮や頂 点の並び替えなどによる最 適化も行なっている。
それらの成 果は、また別の機 会に発 表する。
謝辞
本研 究の成 果は、TSUBAME2.0グランドチャレンジ制 度 、科 学 技 術
振 興 機 構 CREST「ポストペタスケール高性 能 計算に資するシステム
ソフトウェア技 術の創出」から支 援を頂いた。
参考文献
[1]

Graph500 : http://www.graph500.org/.

[2]

Toyotaro Suzumura, Koji Ueno, Hitoshi Sato, Katsuki Fujisawa
and Satoshi Matsuoka, "Performance Evaluation of Graph500
on Large-Scale Distributed Environment", IEEE IISWC 2011 (IEEE International Symposium on Workload Characterization) ,2011/11, Austin, TX, US[3]

Andy Yoo, et al, A Scalable Distributed Parallel Breadth-First[4]

J. Leskovec, D. Chakrabarti, J. Kleinberg, and C. Faloutsos, “Realistic,

Search Algorithm on BlueGene/L. SC 2005.
mathematically tractable graph generation and evolution, using kronecker multiplication,” in Conf. on Principles and Practice of Knowledge Discovery in Databases, 2005.

17

Physis: ヘテロジニアススパコン向け
ステンシル計算フレームワーク
丸山 直也 *

野村 達男 **

* 東京工業大学・学術国際情報センター

佐藤 賢斗 *****Google, Inc.

松岡 聡 **** 東京工業大学・情報理工学研究科

ステンシル計算を対 象とし、TSUBAME2.0 のような大 規 模 GPUスーパーコンピュータを簡 便に利 用可能とするフレームワーク
を提 案する。ステンシルを表 現する関 数を記 述し、それを基にフレームワークが自動 的にGPU 実行コードに変 換する。また通信と計算の最 適化などの多数の計算ノード上の複 数 GPUを効率良く使うための種々の最 適化を自動的に施す。
本 稿では同フレームワークの概 要を報告し、TSUBAME2.0 を用いた評 価 結果より良 好な性能を達 成できることを示す。

はじめに

1

が可能である。一 般に抽 象 度を上げることによって性 能 上のコストが発 生するが、それによって生 産性が大幅に改善されることも多くの場 合重要である。
本 稿ではそのような高い抽 象 度を有したプログラミングモデルの例として、ステンシル計算に特 化したアプリケーションフレームワー

通常のCPUに加えてGPUアクセラレータを共 用した計算がその高い

クPhysisを提 案する。格 子上のステンシル計算では各格 子点につい性能および電力効率より注目をあびている。Tsubameに搭載されててその隣接点を参照することで値を更新する。このような計算パター

いるGPUは 515GFLOPSの性能を持ち、同じくTsubameで用いられてンは数値シミュレーション、特に偏 微 分方 程 式に基づいた計算においるCPUの性能の数 枚倍高速である。また計算速 度だけでなくメモいて頻 繁に出現する。ステンシル計算は典 型的には性 能はシステム

リバンド幅にも優れており、カードあたり150GB/sの速 度を達 成してのメモリバンド幅に律 速される。これはステンシル計 算のB/F 値がいる。これにより計算速 度律速なアプリケーションだけでなくメモリ

GPU 等を含む今日の高性 能計算システムのB/F 値を上回る場 合が多バンド幅 律速なアプリケーションにおいても大幅な性 能向上が可能

いためであり、従ってメモリアクセスの最 適 化が重 要になる。例えば

であり、実際にTsubameにおいて実証されている [2] 。

GPUにおいてはメモリアクセスレイテンシを隠蔽するために非常に多しかしながらそのような異 種プロセッサから構成されるシステム上

くのスレッドの発行やデータのアライメント、キャッシュブロッキングにおけるプログラミングは均 質システムに比べて困 難である。既 存などが知られている [3] 。これらに加えて複 数ノード上のGPUを高効率のプログラミングモデルは低レベルかつ個々の機 種固有な場 合が多に用いるためには通 信と計算のオーバーラップなどの最 適 化が重 要

く、現 状では複 数のプログラミングAPI 、言 語などを共 用したハイブになる。

リッドプログラミングが必 要とされている。例えば、複 数のノード上

Physisフレームワークは実 行 環 境 独 立な可搬 性を有したフレームのGPUを利用するためにはGPU用プログラミングAPIおよび複数ノーワークであり、かつ上 述の環 境 固 有な最 適 化を透 過 的に実 現する。

ド利用のためのメッセージバッシングAPIなどを協調させたハイブリッ

直 交格 子上のステンシル計算をプログラムするために多次 元 格 子のドプログラミングが一 般的である。これによって個々のシステムコン

生 成、データ移動、ステンシルの適 用等を宣言的に表現可能な構文

ポーネントの並列 性を生かした高 効 率なアプリケーションの開 発がを提 供する。直交格子は大域的メモリ空間上に生 成され、その操 作可能になるが、その一方でプログラミングの複 雑さが大幅に増加し、は実 際のメモリシステムの構 成によらず単一のメモリ空 間 上 操 作とプログラム開 発 行 程の増 加につながる。また、正しく動 作するハイ

してプログラム可能である。これにより高い生 産 性および実 行 環 境

ブリッドプログラムを開発させることに加えて、実際に性 能向上を達

独 立性を実現する。また高い抽象 度の宣言的なモデルにより分 散メ

成するためには種々の高度な最 適 化を施す必 要がある。キャッシュ

モリ環 境 上の自動並列化などの高度なコンパイル技 術の適 用を可能

ブロッキングなどの個々の最 適化技 術の多くはこれまでによく知られし、さらに自動チューニングや自動チェックポイントなどのソフトウェた技 術であるが、汎 用プログラミング言 語のコンパイラ等によって自ア技 術の適 用も可能な設 計となっている。

動 的に適 用される場 合は限 定的であり、プログラムを手動で変 更す

本 稿ではPhysisフレームワークのC 言 語をベースにした実装についる必 要がある。GPUクラスタなどの高性 能計算 環 境において高い性て報告する。本実 装では標 準 的なC 言 語を基 盤とし、それに対して能を達 成するためには、複 数のプログラミングモデルを共 用し、さらフレームワークを実 現するための小 規 模な拡 張 構 文を導入したドメにそのような複 雑なプログラム上で種々の最 適 化を適 用する必 要がイン特 化 型言 語（ DSL）
を定義する [4] 。拡張構文を用いたプログラム

あり、生 産性が大幅に低下する。

はPhysisフレームワークによって実 行 環 境 向けのプログラムに自動

この問 題を解決するためにはより抽 象 度の高い統 合されたプログ

変 換される。具体 的にはPhysisプログラムをGPUクラスタ上で実 行

ラミングモデルが必 要である。高い抽象度により生産性を向上させ、する場 合はフレームワークによってMPIおよび CUDAを用いたソース

かつ実 行 環 境によらない可搬 性のあるプログラミングモデルの実現

コードに自動 変 換される。またこの際に通 信と計算のオーバーラップ

18

SC’11 テクニカル・ペーパー

Physis: ヘテロジニアススパコン向け
ステンシル計算フレームワークなどの最 適化も自動的に適 用される。
本フレームワーク実 装の有 効 性を評 価するためにステンシル計算
をフレームワーク上に実装し、その性能をTSUBAME2.0 の256GPUをプログラミングモデル

用いた評 価を行う。本 稿ではその結果の一 部について報告し、良 好3

な性 能およびスケーラビリティが達 成可能であることを示す。より詳細な結果については文 献 [1] を参照されたい。

Physis DSLは標 準的なC 言 語にステンシル計算用のデータ型とイントリンシックスを追 加したものである。本 拡 張を用いて記 述されたステンシル計算フレームワーク

2プログラムはソースコード変換 器によって実行 環 境向けのソースコー
ドに変換される。
Physisでは浮 動 小 数 点値を持った直 交のカルテシアン多次 元 格
子を利 用可能である。多次 元 格 子を表 現するために格 子の型と次元に基づいたデータ型として PSSGrid3DFloatやPSGrid2DDouble

ステンシル計算において高い生 産 性の実 現を狙ったフレームワークなどを導入する。前者は float 型を格 子データとして持つ３次 元 格 子を設 計する。同フレームワークはドメイン特化型プログラミング言 語であり、後者はdouble 型の２次 元 格子である。これらの型の内部 構（DSL）およびアーキテクチャ固有ランタイムから構成される。DSLは造はユーザプログラマからは隠蔽されており、実 態へのハンドルとし宣言的かつ機種独立にステンシル計算を記述可能であり、ソースコーて機能する。

ド変換 器によって実行アーキテクチャ向けコードに変換される。ランPhysisでは次 元および格 子点の型によって異なる複 数の格 子 型がタイムは多次 元 格 子データを簡 便かつ機 種 独 立に操 作するための提 供されているが、以下では説明の簡便 化のために格 子の型の名前

抽 象 化を提 供するものである。本 節ではフレームワーク設 計におけとしてPSGridを用いる。

る主な目標を述べる。

PSGridFloat3D 型の格子は PSGridFloat3DNewおよびPSGridFree
により作成、破 棄される:
PSGrid3DFloat PSGrid3DFloatNew(

自動並列化 : Physis DSLは分 散メモリ環 境を含む様々な環 境においsize_t dimx, size_t dimy,て自動並列 化が可能なものとして設 計する。一 般には自動並列 化、

size_t dimz,

特にデータの局所 性を考慮した並列 化は非常に困 難であり、汎 用言

enum PS_GRID_ATTRIBUTE attr)

語では分散メモリ環境上などでは現実的ではないとみなされている。

void PSGridFree(PSGrid g)

我々はDSLを特定の計算パターンに限 定し、必 要な制約を導入するこ

格 子の各 要 素には以下のイントリンシックスにより一 括および各点とで本 DSLで記 述されたプログラムの自動並列 化が可能となるよう毎にアクセス可能である:にDSLの設 計を行う。
組み込み型 DSL 設 計の採 用 : 全く新しい言 語を設 計することで、対 象

const void *src)の問 題に対して高度に最 適 化されたプログラミングモデルを提 供す

void PSGridCopyout(PSGrid g, void *dst)

ることが可能である。しかしながら、現 実的には既存プログラミング

PSGridGet(PSGrid g, size_t i,

言 語と大きく乖離することは普及の妨げにもなるため、我々はPhysis

size_t j, size_t k)

におけるDSLを広く使われている既存汎 用言 語に対してステンシル計

void PSGridSet(PSGrid g,

算用拡 張を施した組み込み型 DSLとして設 計する。本 稿ではそのよ

size_t i, size_t j, size_t k, T v)

うな広く普及した言 語としてC 言 語を基盤とする。

void PSGridEmit(PSGrid g, T v)

宣 言 的かつ高い記 述 性を持った言 語 : 高い生産性を実現するために

PSGridGet およびPSGridSetのsize_t 型のパラメータは同じくパラ

高い抽象度を持った宣言的プログラミングモデルを採用する。宣言的メータとして与えられる格 子における位 置を表し、格 子の次 元と同とすることで命令的モデルと比較して記述量を削減可能である。例え数の 個 数のパラメータを持つ。PSGridGet の 返り値やPSGridSet,

ばPhysis DSLではプログラマは各格子点の計算について記述するが、PSGridEmitのパラメータvは格子点の型と同一である。

格子全体がどのように計算されるかはフレームワークによって実行環
境毎に最適となるように決定される。高い抽象度のプログラミング言
語はその実装一方で抽象度を必要以上に高めると複雑な現実のシミュ
レーションの記述が困難になりうる。我々は生産性と性能を両立するために対象ドメインのアプリケーションの記述に必要な構文を吟味した設計を行う。

19

void PSGridCopyin(PSGrid g,

SC’11 特集号

３-1

イクセションバンド幅ネットワークを実現するファットツリートポロジ

例を構成する。ソフトウェアとしてはOSとしてSUSE Linux Enterprise

図 1は２次 元 格子における９点ステンシルを実装したステンシル関数である。ステンシル関 数の実際の格 子への適 用はPSStencilMapおよびPSStencilRunを用いて図 2 のようにプログラム可能である

Server 11 SP1を用い、CUDA v3.2 、gcc/g++ v4.1.2 を用いた。
４-1 ウィークスケーリング
図３はdif fusion のウィークスケーリングにおける性 能を示したものである。赤 線および青 線 がそれぞれGPUあたり256x128x128 、
512x256x256 の 計 算を行った場 合であり、Y 軸 および Z 軸 の２次
元 分 割 を行 った。結 果 よりGPUあたりの デ ー タサイズ が大 き い
512x256x256 の場 合はほぼ 256GPUまで性 能がスケールすることがわかる。一方で、256x128x128 の場 合でも 256GPUまで良 好に性 能向上していることがわかる。

図 1 9 点ステンシルの例

図 2 ステンシル関 数の格 子への適 用 例

4

性能評価

図 3 Dif fusion のウィークスケーリング性 能

図４はSeismicのウィークスケーリング評 価 結 果である。各 GPUは提 案フレームワークの有 効 性を評 価するためにプロトタイプ実 装を2563 の領域を計算する。Diffusionとは異なりX 軸および Y 軸にて領行った。実装にはROSEコンパイラフレームワーク を用い、ソースコー

域 分 割を行った。そのためストライドアクセスをともなう境 界領 域

ド変換器およびランタイムから構成される。変換はCPU実行向けコー

交換が発 生し、その結果としてDiffusionと比 較してスケーラビリティ

ドの場合はCコードを生成し、GPUの場合はCUDAコードを生成する。

が劣ることがわかる。また 64GPUを用いた場 合に性 能低下が見られ

また分 散メモリ環 境の場 合はMPIにより並列化を行う。

るがその原因についてはより詳 細な解 析が必 要である。

[5]

性能評 価のために以下の2 つのステンシルをPhysis DSL 上で実装した。
・Diffusion: 3 次 元 7点ステンシル・Seismic: 3 次 元 地震波伝 播シミュレーション
Dif fusionは単一のステンシル関 数から構成される比 較 的小規 模なプログラムとなっているが、seismicは 27 個のステンシル関数およびスタッガード格子を用いて構成される。27 個の内、6 個は 3 次 元 格子の2 次 元 境界平面の計算であり、PSDomainオブジェクトによって格子 全 体に対して部分 的に適 用される。プログラムは 2 つとも単 精度
浮動小 数 点を用いる。
評 価 環 境として TSUBAME2.0 を用いる。同マシンは 1408 台の計算ノードから構成され、各ノードはIntel Xenn Westmere-EP 2.9GHz

図 4 Seismicのウィークスケーリング性 能のCPUを 2 基、NVIDIA M2050 GPUを 3 基 搭載する。メモリはシステ
ムメモリとして 52GB 、GPUメモリとして各 GPUあたり3GB 搭載している。計算ノードは 2 本のQDR Infinibandによって接 続され、フルバ

20

SC’11 テクニカル・ペーパー

Physis: ヘテロジニアススパコン向け
ステンシル計算フレームワーク

ムワークが対 象とする構 造格 子に対してListzを適 用することも可能

４-2 ストロングスケーリングである。しかしながらListzは非 構 造格 子に最 適 化されているため

図 5 は問題サイズ512x512x4096 のDiffusionの性能を示したもので構 造格子において高い効率を達 成するための最 適化が可能かどうかある。128GPUまで用い、1次 元、2 次 元および 3 次 元 分割の場 合をは明らかではない。

比 較した。
1次 元 分割ではZ 軸 方 向に均 等に分割し、2 次 元 分割ではY 軸お
よび Z 軸、3 次 元ではすべての軸で均 等分割を行った。その結果、利用 GPU 数が増えるにつれて多次 元 分割が性能 上優れていることがわおわりにかる。このような多次 元 分割の優 位性自体はすでに広く知られたも6のであるが、本フレームワークを用いた場 合は明 示的に分割を行うプ
ログラミングをせずにフレームワークによって自動的に処理される。

TSUBAME2.0 のような大 規模 GPUスーパーコンピュータにおけるプロ
グラミングの生 産 性を向 上させるためにステンシル計算 向けフレー
ムワークPhysisを提 案した。C 言 語に基づいたPhysis DSLによりス
テンシル計算を簡潔かつ実 行 環 境 独 立に記 述することを実現した。
DSL 変 換 器によって対 象 実行 環 境向けにコード生 成を行い、通 信と計算のオーバーラップなどの最 適 化を自動 的に適 用する。本 稿ではプロトタイプ実 装を行い、TSUBAME2.0 の256GPUを用いた評 価では良 好な性 能を達 成できることを示した。今後はさらなる性 能向上を目的としたフレームワークの拡張を施し、またアプリケーションケー
ススタディを通じて本フレームワークの有効性の実証を行う予定である。
謝辞
本 研 究の一 部は科 学 技 術 振 興 機 構 CREST「高 生 産 性・高性 能ア

図 5 512x512x4096 の Dif fusion の性 能

関連研究

プリケーションフレームワークよるポストペタスケール高性 能 計算の5実 現」、科 学 研 究 費 補 助 金（ 22700047）、JST-ANR FP3C 、NVIDIA
CUDA Center of Excellenceから支 援を受けた。
参考文献
[1]

“ Phy sis: an imp licitl y p ar all el p ro gr ammin g m o d el fo r

れており、高速 化だけでなく電 力 効 率の向上が可能であることが示

s ten cil co mp u t atio ns o n larg e - s c al e GPU -a ccel er ate d

されている。しかしながらアクセラレータを使うことによりプログラs u p e r c o m p u t e r s , ” I n P r o c e e d i n g s o f 2 0 11 AC M / I E E E

ミングの複 雑さが増し、生 産 性の大 幅な低下が大きな問 題となってInternational Conference for High Performance Computing,いる。この問 題に対して生 産 性向上を目的とした取り組みが行われ

Networking, Storage and Analysis (SC'11), Seattle, WA, USA (2011)

ている。

[2]

T.Shimokawabe, T.Aoki, C.Muroi, J.Ishida, K.Kawano, T.Endo,

Mintは指 示 文によるステンシル計算フレームワークである [6] 。ス

A.Nukada, N.Maruyama, S.Matsuoka, “An8 0 -fold speedup,

テンシルのループに対して指 示 文を追 加することでループをGPU上

15.0 TFlops full GPU acceleration of non-hydrostatic weather

で実 行するコードに変 換する。YpnosはHaskell 言 語に基づいたス

model ASUCA production code” in Proceedings of the 2010

テンシル用 DSLであり、GPU 用に自動並列 化を行う [7] 。MintおよびACM/IEEE International Conference for High Performance

Ypnosともに我々と同じく自動並列化等の高生 産性のための機 能をComputing, Networking, Storage and Analysis, (SC’10), New

実現しているが、双 方とも単一 GPUのみを対 象としている点が我々と異なる

21

N . M a r u y a m a , T. N o m u r a , K . S a t o , a n d S . M a t s u o k a ,GPUアクセラレータを使った科 学 技 術計算に関する研 究が多く行わ

Orleans, LA, USA (2010)[3]

S. Ryoo, C. I. Rodrigues, S. S. Baghsorkhi, S. S. Stone, D. B. Kirk, and W. M.

Listzは非構造格子計算のためのDSLである [8] 。我々のフレームワー

W. Hwu, “Optimization principles and application performance

クと同 様に高い抽 象 度を持ったDSLであり、実 際の実 行 環 境 上でのevaluation of a multi- threaded GPU using CUDA,” In Proceedings

実 装 方 式はユーザプログラマからは隠 蔽されている。これによってof the 13th ACM SIGPLAN Symposium on Principles and practice

対象ドメインに特化した最適化を適用可能にしている。我々のフレー

of parallel programming, PPoPP ’08, zpages 73–82, (2008)

SC’11 特集号[4]

M. Fowler. Domain-Specific Languages. Addison- Wesley
Professional (2010)[5] M. Schordan and D. Quinlan, “A source-to-source architecture
for user-defined optimizations,” In Modular Programming
Languages, volume 2 7 8 9 of Lec ture Notes in Computer
Science, pages 214–223. Springer Berlin / Heidelberg (2003)[6]

D. Unat , X . C a i, a n d S . B . B a d e n . “M int : rea l i z in g CUDA
per formance in 3 D stencil methods with annotatedc ,” In Pro ce e din gs o f th e Inter nati o nal Co nferen ce o n
Supercomputing (ICS’11), ICS ’11, pages 214–224 (2011)[7] D. A . O rchard, M . B olingbroke, and A . Mycrof t. “ Ypnos:
declarative, parallel structured grid programming,” In DAMP ’10:
Proceedings of the 5th ACM SIGPLAN workshop on Declarative
aspects of multicore programming, pages 15–24 (2010)[8]

H . Chaf i, Z . D eV i to, A . M o o r s , T. R o mp f, A . K . Su - j e eth,
P. Ha nr a ha n, M . O d e r sk y, a n d K . O lu kot u n, “ L a n gua g e
vir tualization for heterogeneous parallel computing,” In
Proceedings of the ACM international conference on Object
oriented programming systems languages and applications,OOPSLA ’10, pages 835–847 (2010)

22

FTI : ヘテロジニアススパコン向け
耐障害インタフェース〜100TFlops 超 東北地方太平洋沖地震シミュレーション〜
Leonardo Bautista-Gomez* Dimitri Komatitch**
Franck Cappello† 松岡 聡‡ 中村 武 ***丸山 直也 *

坪井 誠司 ***

* 東京工業大学・情報理工学研究科 ** University of Toulouse, Observatoire Midi-Pyrénées
*** JAMSTEC † University of Illinois, INRIA ‡東京工業大学・学術国際情報センター

TSUBAME2.0 のようなペタフロップ級の巨大なシステムは、計算資 源の障害も無視することはできない。
長時間にも及ぶ大 規模 科 学シミュレーションを行うためには、たとえ障害が発 生したとしても計算を続ける必 要がある。
我々は、地 震 波シミュレーションにおいて、リード • ソロモン • エンコーディングを利用することにより、システムの信 頼性を保 障しながら100 TFlops 以 上の性能をTSUBAME2.0 で実現した。

はじめに

1

リード • ソロモン符号を利用した
チェックポイント2

スーパーコンピュータなどのHPCシステムは、本来高信 頼に設 計されこのため、以前よりマルチレベル•チェックポイント/リスートという手るが、
近 年、システムを構成する機 器の増大や複 雑 化により、障害発法が提 案されている。これは、チェックポイント先として並列ファイ生率が増えてきており、現在、ペタスケールシステムでは、その平均 故

ルシステムだけでなく、計算ノード上のメモリ領域やローカルストレー

[11]

障間隔
（MTBF: Mean Time Between Failures）
が数日となっている 。

ジを活用する手法である。実際に、チェックポイントを保 存する際にそのため、長時間に及ぶ大 規 模 科 学 計算を行えるよう、耐障害性のは、並列ファイルシステムより高いI/Oバンド幅をもつメモリ領域やロー

高いシステムの構 築がより一層重 要となっている 。チェックポイント/

カルストレージに対して、高い頻 度でチェックポイントを保 存し、低い

リスタート
（CR）
は、耐障害 手法として、広く利用されている。これは、バンド幅ではあるが、信頼性の高い並列ファイルシステムに対しては、

実 行中のアプリケーションの状 態をチェックポイントと呼ばれるファ

それより低い頻 度でチェックポイントを保 存することにより、信 頼性

イルとして並列ファイルシステムなどの高 信 頼ストレージに保 存してと性能を両立させる手法である。実際、年々増加するプロセッサの高

おき、障 害が発 生したときに、最 新のチェックポイントから計算を再い集 積 度に起因するソフトエラー率（ SER: Soft Error Rate）が増える

開させる手法である。

と報告されている [4,10] が、このような、軽 度な障害は計算ノード上のメ

大 規 模 科 学アプリケーションを対 象とした場 合、その状 態を保 存

モリ領 域やローカルストレージに保 存されたチェックポイントから計するためには、
数 万オーダーに及ぶ数のプロセスが 、それぞれ数 GBの算をリスタートさせることが可能である。不 揮発 性メモリである相変チェックポイントデータを書き込むため、合計 数十 TBのデータが書き

[2]
化メモリ
（ PRAM : Phase-change RAM）
やTSUBAME2.0 にも搭載

込まれることとなる。しかし、並列ファイルシステムの低いI/Oバンド[6]
されているSSD（Solid State Drive）
が、このような用途に最 適な保幅を考えると、
チェックポインティングに要する時間が爆発的に増える。

存領域となる。 一方、チェックポイントの複製やパリティを他の計算

実際に、現在のペタスケールシステムにおいて、大規模科学アプリケー

ノードに分 散させて保 存することにより、計算ノード単体の故 障が起ションを持 続 的に実 行できる頻 度でチェックポイントを取った場 合、きた時に、故障ノードのチェックポイントを新しく割り当てられた計算

その25%の時間をチェックポイントに費やしてしまう事例もある [9] 。

ノード上へ転 送することにより、その計算ノード上で計算をリスタート

さらに、将 来のポストペタ • エクサスケールシステムでは、構成するすることが可能である [1, 7] 。しかし、一 般にマルチレベルチェックポ

機 器の数がさらに増 加すると予見される。それに伴い 、システム全イント•リスタートでは、複 数の計算ノードが一度に故障する場 合に備

体の障 害 発 生率が上 昇し、MTBFが数 分足らずになると予 想されてえて、やはりある程 度の頻 度でチェックポイントを並列ファイルシス

いる。このため、より高い頻 度でチェックポイントを取る必 要がある

テムに保 存する必 要がある。しかし、先に述べたようにそれに起因すが、チェックポインティングに費やす時間が膨 大となってしまうため 、るオーバーヘッドは大きいという問題が残る。

ポストペタ • エクサスケールシステムではアプリケーションの計算が実 質的に進まなくなると危惧されている。

そこで 我々は、トポロジー 情 報 から計 算ノード をグル ーピング
し 、各 々のグル ープ 内 で、
リード • ソロモン • エンコーディング（ RS
Encoding: Reed-Solomon Encoding）
を用いたチェックポイントの符 号 化を適 用することにより、並列ファイルシステムへチェックポイン
トを保 存することなく、
複 数の計算ノードの障害にも耐え得る耐 障害
インタフェース（FTI: Fault Tolerance Interface）
を提 案する。しかし、RS Encodingの複 雑性により、XORによるパリティ符 号などの既存の23

SC’11 特集号

手法よりエンコーディングに多くの時 間を要する。一方、多くのGPU

TFlops（ 384 GPUs）へストロングスケールしており、FTIを使 用したアプリケーションでは、ほとんどの計算をGPU上で行うため、いくつ

場 合でも、
わずか 4%ほどのオーバーヘッドでストロングスケールする

かのCPUコアは使われていない。そこで、この使われていないCPUコ

ことを確 認した。一方、Lustre へチェックポインティングでは、I/Oバ

アをエンコーディング専用のコアとして使 用することにより、エンコー

ンド幅がボトルネックとなり、スケールしない。Flops 値は 5 回の計 測

ディングの時間を隠 蔽する。我々は、性 能モデルに基づき、チェックで得られた実行時間の平均とPAPI [14] から得られた浮動 小 数 点 演 算

ポイントサイズ、グループサイズ違いによる、エンコーディング/デコー数から求めた。

ディング時 間への影 響やスケーラビリティの検 証を行い、このFTIを次に、1000 GPU までウェークスケールさせ、より大 規 模な問題サ

2011年 3 月11日東 北 地 方 太平洋 沖 地 震のシミュレーションに適 用

イズで評 価を行った。Fermi GPU [5] では 12.5%のメモリ容量は ECCした。この計算において1000 個のGPUを使 用し100 TFlops 以 上のに使 用されるため、１GPU で使 用可能な 2.6GBのオンボードメモリ

性能を達 成し、実アプリケーションにおいて、低オーバーヘッドかつ高のうち、2.1GB をアプリケーションの実 行に使 用した。また、チェッ

信 頼なシミュレーションを実現した。

クポイントの間 隔は、Youngの性 能モデル から計 算される最 適な
時 間 間 隔 6 分に固 定した。この評 価では、チェックポイントを行わ

100TFlops 超 高信頼 SPECFEM3D

3ない場 合（ No ckpt.）、ローカルの SSD にチェックポイントを保 存す
る場 合（L1）、L1 に加え、2 回に１度 FTIを使 用しRS Encodingを行う
場 合（FTI-L1,L2）、FTI-L1,L2 に加え、6 回に1度、チェックポイントを
Lustre へ保 存する場 合（FTI-L1,L2,L3）、最後にBLCRを用いてLustre
へ保存する場合（BLCR+Lustre）の５つの場合を比較した。BLCRでは、

スケーラビリティ及び効 率性を検 証するために、FTIをSPECFEM3D

GPUアプリケーションのチェックポインティングはサポートしていないに組み込み、評 価を行った。SPECFEM3Dは地 震 波の伝 搬をシミュため、各々のプロセスが 2.1GBのデータを書き込むようにし、擬 似 的レートするアプリケーションであり、世 界 300 以 上の研 究グループにチェックポイントを行わせた。また、BLCRはカーネルレベルのチェッ

で使 用されている。このアプリケーションは、2002 年 時点で、世界クポイントを取るため、アプリケーションが使 用するメモリ領 域 全て１位であった地 球シミュレータ上で CPU 1944 個を使 用した計算でがチェックポイントとして保 存される。SPECFEM3Dにおいて、これ

5 TFlops を出し、その後、2003 年 Gordon Bell Super Computing

は、アプリケーションレベルで保 存されるチェックポイントデータの 5Award for Best Performance [8] を受賞したアプリケーションである。

倍のサイズにあたる。図 2 にその結果を示す。Flops 値は 3 回の計 測

さらに、SPECFEM3DはTSUBAME2.0 のようなCPU/GPUハイブリッで得られた結果を平均した値である。SPECFEM3Dでは、43 TFlops

ト型システムの性 能を十 分 発 揮できるように、CUDA [12, 13] による、（384GPUs）
から117TFlops（1152GPUs）へウィークスケールしており、GPUへ移植がなされた。一 般に、このような地 震 波伝 搬シミュレー

L1、FTI-L1、L2 の場 合でも、およそ 8%のオーバーヘッドでウィークス

ションでは、単 精度 浮動小 数 点が用いられており、また、このGPU 版

ケールすることを確 認した。FTIでは、アプリケーションで使 用されてSPECFEM3Dは、有限 差 分 法や有限 要素 法の計算のようにメモリバ

いないCPUコアをエンコーディングに使 用しているので、エンコーディ

ウンドな計算特性を持つため、性能がメモリバンド幅律速になること

ングによるオーバーヘッドを低く抑えることができる。FTI-L1、L2 、L3を強調しておく。

では、さらに 3 % のオーバーヘドが生じた。これは 、Lustreへのチェック

TSUBAME2.0 を使 用した評 価を行い、ストロング/ウィークスケール

ポインティングに起 因する。一方、BLCR+Lustreではそのチェックポ

させたときの、SPECFEM3D の 性 能を検 証する。TSUBAME2.0 のイントサイズとLustre へのチェックポインティングが原 因でウェーク

アーキテクチャ及び性 能を表 1に記す。各評 価では、SPECFEM3Dのスケールしていない。また、この評 価では、6 分間隔でチェックポイン

30~40 分実行時間を要した。

トを行ったが、このような短い時 間 間 隔でチェックポインティングを

まず、ストロングスケーリングの結 果を図 1に示す。ストロングス
ケールでは、固 定された問 題サイズに対し、GPUの数を増 加させる

行った場 合でも、実アプリケーションにおいて、100 TFlops 以 上の性
能を達 成している。

ため、１GPUあたりのチェックポイントサイズは減少する。このため、
チェックポイントサイズに応じて、チェックポインティングの間隔を短くし、障害が発 生後にリスタートする際、障害発 生時点までの復旧時間を短くする。この評 価では、チェックポイントを行わない場 合（ No
ckpt.）、FTIを用いてチェックポイントを行う場合（FTI ckpt.）、並列ファ
イルシステム Lustre 上にチェックポイントを保 存する場 合（ Lustre
ckpt）の 3 つの場 合を比 較した。また、FTIはアプリケーションレベル
でチェックポイントを行うため、リスタートに必 要なデータのみ保 存する。そのデータ量はアプリケーションで使 用するメモリ量の20%に相当する。図１より、SPECFEM3D では、5 TFlops（48 GPUs）から 23

24

SC’11 テクニカル・ペーパー

FTI : ヘテロジニアススパコン向け耐障害インタフェース〜 100TFlops 超 東北地方太平洋沖地震シミュレーション〜

表 1 TSUBAME2.0アーキテクチャ

図 1 SPECFEM3D per formance (Strong scaling)25

図 2 SPECFEM3D per formance (Weak scaling)SC’11 特集号

東北地方太平洋沖地震
シミュレーションへの適用

4

我々は、FTIを組み込んだSPECFEM3Dを用いて、東 北 地 方 太平洋
沖 地 震 のシミュレーションを試 みた。このシミュレーションでも
TSUBAME2.0 を活用する。地 震記象（地 震計に記 録された地 震 動の
波 形 記 録）は、震 央 距 離が 30°~100°の範囲のIRIS GSN 及びIFREE
OHP 広域 地 震計で記 録された波 形データを使 用し、その地 震記 象に対しバンドパスフィルターをかけ 0.002Hz~1Hz の周 波 数 帯の波形を取り出す、その後、中村ら [16] と同様の方法で、震 源における滑り
分布を得るために、非負の拘束 条 件下で波 形インバージョンを適 用した。震 源や観 測点付近の波 動 場の計算には、IASP91の伝 達 速 度
構 造モデルを使った 。また、CMT（ Centroid Moment Tensor）モデ

図 3  地 殻 変 動（北 方向）

ルに基づき、走向をN201°E 、断層角を 9°と設 定し、断層長を余 震の発 生分布から 440kmと仮 定する。
上 記の条 件下で、波 形インバージョンを行った結 果、震 源から北
東 及び南 西 100km 付近に 40mほどの逆 断層によるずれあり、また北 東と南 西に向かって、バイラテラル破 壊（中 央から外側へ進む断層破壊）が発 生したことがわかった。また、同様の条 件下で、スペクトル要 素 法 [15, 16] を用いて福島県の広 野 町の地 震 計における、理 論
地 震記 象を計算した。その理 論 地 震記 象を図 3 、4 、5 に示す。
これらの図から、今 回の地 震により、この地 震 計 付近では、南へ 2m 、東へ
2m 、鉛 直 方 向に 1m 静 的 変 位したことがわかる。これは観 測された地殻 変 動の結果とも概ね一 致している。

図 4  地 殻 変 動（東 方向）

図 5  地 殻 変 動（鉛直 方向）

26

SC’11 テクニカル・ペーパー

FTI : ヘテロジニアススパコン向け耐障害インタフェース〜 100TFlops 超 東北地方太平洋沖地震シミュレーション〜

おわりに

5

simulation on the Earth Simulator, Proceedings of the ACM /IEEE Supercomputing SC’2003 conference, November 2003.
[9]

G. Grider, J. Loncaric, and D. Limpart, Roadrunner System
Management Report, Los Alamos National Laboratory, Tech.
Rep. LA-UR-07-7405, 2007.

リード • ソロモン • エンコーディングを用いた耐 障 害インタフェース[10] S. Y. Borkar, Designing Reliable Systems from Unreliable

FTIを紹 介した。エンコーディング専 用のスレッドを割り当て、エン

Components: The Challenges of Transistor Variability andコーディング時 間を隠 蔽することによりSPECFEM3Dにおいて 100

Degradation, IEEE Micro, vol. 25, no. 6, pp. 10-16, 2005.

TFlops 以 上の性 能を達 成した 。また、FTIを組み込んだSPECFEM3Dを用いて、東北 地方太平洋沖地 震のシミュレーションを行い、実アプ
リケーションの高信 頼計算を実現した。

NVIDIA graphics cards using CUDA, Journal of Parallel and

本研 究の一 部は JSPS 、ANR/JST FP3Cプロジェクト、INRIA-Illinois
Joing Laborator y for Petascale Computing から支 援を頂いた。
記して謝 意を表す。
参考文献

X. Dong, N. Muralimanohar, N. Jouppi, R. Kaufmann, Y. Xie.

method, Beowulf computing, and global seismology, Science
298, 1737-1742, 2002.

Leveraging 3D PCRAM Technologies to Reduce CheckpointO v e r h e a d f o r Fu t u r e E x a s c a l e Sy s t e m s . I n AC M / I E E E

process of the 2008 Wenchuan, China earthquake inferredInternational Conference for High Performance Computing,from teleseismic waveform inversion and forward modeling

Networking, Storage and Analysis, Portland, 2009.

of broadband seismic waves, Tectonophysics, vol. 491, 72-

Z. Cheng, J. Dongarra, A scalable Checkpoint Encoding

84, 2010.

Algorithm for Diskless Checkpointing. Proceedings of the 11th

[17] S. Tsuboi, D. Komatitsch, C. Ji, J. Tromp, Broadband modellingIEEE High Assurance Systems Engineering Symposium, HASE

of the 2002 Denali fault earthquake on the Earth Simulator,

2008, Nanjing, China, December, 2008.

Phys. Earth Planet. Inter. 139, 305-312, 2003.

B. Schroeder, E. Pinheiro, W. Weber. DR AM errors in themodeling of computer systems (SIGMETRICS), ACM, New
York, NY, USA, 2009.
[5] NVIDIA Corporation http://www.nvidia.com/object/fermi_
architecture.html
A. Moody, G. Bronevetsky, Scalable I/O Systems via NodeLocal Storage: Approaching 1 TB/sec File I/O. DOE technicalreport, 2009[7] W. D. Gropp, R. Ross, and N. Miller. Providing ef ficient I/O
redundancy in MPI environments. Lecture Notes in Computer
Science, 3241:7786, September 2004
D. Komatitsch, S. Tsuboi, C. Ji and J. Tromp, A 14 . 6 billion
degrees of freedom, 5 teraflops, 2 . 5 teraby te earthquake

27

icl.cs.utk.edu/papi/[15] D. Komatitsch, J. Ritsema, J. Tromp, The spectral-element[16] T. Nakamura, S. Tsuboi, Y. Kaneda, Y. Yamanaka, Rupture

11th international joint conference on Measurement and[8]

MPI on a large GPU cluster, Journal of Computational Physics,[14] PAPI: Performance Application Programming Interface, http://

w il d: A L arg e -S c al e Fi el d St u d y. In Pro ce e din gs o f th e

[6]order finite-element seismic wave propagation modeling with

Design, Modeling, and Evaluation of a Scalable Multi-level

Analysis, New Orleans, 2010

[4]j.jpdc.2009.01.006, 2009.
[13] D. Komatitsch, G. Erlebacher, D. Göddeke, D. Michéa, High-

vol. 229(20), p. 7692-7714, doi:10.1016/j.jcp.2010.06.024, 2010.

for High Performance Computing, Networking, Storage and[3]Distributed Computing, vol. 69(5), p. 451-460, doi: 10.1016/

A . Mo o d y, G . Bronevet sk y, K . Mohror, B. R . de Supinsk i,
Checkpointing System. In ACM/IEEE International Conference[2]

Director’s Colloquium, LANL, May 2004.
[12] D. Komatitsch, D. Michéa, G. Erlebacher, Por ting a highorder finite-element earthquake modeling application to

謝辞

[1]

[11] D. Reed, High - End Computing: T he Challenge of Scale,SC’11 特集号● TSUBAME e-Science Journal No.5
2012 年 2 月28 日 東京工業大学 学術国際情報センター発行 ©
ISSN 2185-6028
デザイン・レイアウト：キックアンドパンチ
編集： TSUBAME e-Science Journal 編集室
青木尊之 ピパットポンサー・ティラポン
渡邊寿雄 佐々木淳 仲川愛理
住所： 〒 152-8550 東京都目黒区大岡山 2-12-1-E2-1
電話： 03-5734-2087 FAX：03-5734-3198
E-mail： tsubame_j@sim.gsic.titech.ac.jp
URL： http://www.gsic.titech.ac.jp/28TSUBAME 共同利用サービ ス『みんなのスパコン』TSUBAME 共同利用サービスは、
ピーク性能 2.4PFlops 、18000CPUコア、4300GPU 搭載
世界トップクラスの東工大のスパコンTSUBAME2.0 を
東工大以外の皆さまにご利用いただくための枠 組みです。

  課 題 公募する利 用 区 分とカテゴリ
共同利用サービスには、「学術利用」「産業利用」、「社会貢献利用」、の3つの利用区分があり、さらに「成 果公開」と「成 果非公開」のカテゴリがあります。
ご利用をご検 討の際には、下記までお問い合わせください。

他大学や公的研究機関の研究者の 学術利用［有償利用］

TSUBAME

民間企業の方の 産業利用［有償・無償利用］

共同利用とは…

その他の組織による社会的貢献のための 社会貢献利用［有償利用］

  共同利用にて提 供する計算資源
共同利用サービスの利用区分・カテゴリ別の利用課金表を下記に示します。TSUBAME 2.0
における計算 機 資 源の割 振りは口 数を単 位としており、1 口は標 準 1ノード
（12 CPUコア,
3GPU, 55.82GBメモリ搭載）の 3000 時間分（≒約 4 ヵ月）相当の計算 機資 源です。
1000 CPUコアを1.5日利用する使い方や、100 GPUを3.75日利用する使い方も可能です。
利用区分

利用者

制度や利用規定等

カテゴリ

利用課金

学術利用

他大学または研究機関等共同利用の利用規定に基づく

成果公開

1口：100,000 円

産業利用

民間企業を中心としたグループ「先端研究施設共用促進事業」に基づく非営利団体、公共団体等共同利用の利用規定に基づく

社会貢献利用

成果公開

トライアルユース（無償利用）
1口：100,000 円

成果非公開

1口：400,000 円

成果公開

1口：100,000 円

成果非公開

1口：400,000 円

  産業利用トライアルユース制度（先端研究施設共用促進事業）
東 工 大のスパコンTSUBAMEを、より多くの企 業の皆さまにご利 用いただくため、初めて
TSUBAMEをご利 用いただく際に、無 償にてご試 用いただける制 度です。
（文部 科 学省 先 端 研 究 施 設 共 用促 進事業による助成）

詳しくは、下 記までお問い合わせください。

お問い合わせ●東京工業大学 学術国際情報センター 共同利用推進室●e-mail kyodo@gsic.titech.ac.jp Tel. 03-5734-2085 Fax. 03-5734-3198
詳しくは http://www.gsic.titech.ac.jp/tsubame/をご覧ください。