14

大規模GPUコンピューティングによる天然岩石内の二相流シミュレーション

Two - phase Porous Flow Simulation for Natural Sandstone on GPU Supercomputer

大規模画像データセットの機械学習のための分散コンピューティング

Distributed Computing for Machine Learning on Large -Scale Image Dataset

マルチGPU超並列クラスタシステムを用いた大規模ナノ炭素分子の電子状態計算
Electronic structure calculation of large nano carbon molecules
using multi - GPU massively parallel cluster system

大規模 GPUコンピューティングによる天然岩石内の二 相流シミュレーション
辻健

蒋飛

九州大学カーボンニュートラル・エネルギー国際研究所

GPUスパコン（ TSUBAME 2.5 ）を用いることにより、格 子ボルツマン法（ LBM）による２相の間 隙 流 体シミュレーションを、巨大なデジタル岩石（1000 3 グリッド ）に適 用することが可能となった。計算領域のサイズアップによって、これまで間隙 スケー
ル（μ mスケール）で実 施していた計算を、岩 石スケール（ mmスケール）で実 施できるようになった 。それによりLBMの計算
結 果を、実 験 結 果と比 較し、計算 結 果を検 証することが可能となった 。岩石間 隙 内の２相 流を扱う例として、CO 2 地中貯 留がある。このプロジェクトでは、貯 留層内での水 — CO 2 の２相 流 挙 動が重 要となる。その挙 動 は貯 留層の条 件（ 例えば界面張 力、間 隙 形 状、濡れ性、水 圧など）に強く依 存する。今 回 開 発したアプローチを用いて、様々な貯 留層条 件で CO 2 の挙 動を計算することで、効 果 的にCO 2 を貯 留できる条 件を明らかにした 。さらに貯 留層に圧 入したCO 2 の挙 動を考慮して、CO 2 が鉱物 化する過 程をモデル化し、それが浸 透 率といった水 理 特 性に与える影 響を評 価することも可能となった 。

1

はじめに

しかし多くのLBMシミュレーションでは、計算効率の制限から、
低 解 像 度かつ小さなスケールで実 施されることが多かった。例えば、Ahrenholz et al.[6] は、11 μ mの解像度で 200 3 グリッドの砂岩モデルに対してLBMの計算を実施している。またRamstad et al.[7]も、6μmの解 像 度で 256 3 グリッドの砂 岩などのデジタルロック

多孔質媒質内の流体挙動は、理学から工学にわたる様々な分野でに対して、LBMを適用している。

取り扱われる。特に、岩 石 内 部の２相 流の挙 動は、石 油の増 進 回

近年、Graphics Processing Unit（ GPU）が科学計算にも利用さ

収（EOR）や、二酸化炭素（CO 2 ）の回収・貯留（CCS）、難水溶性汚れ、その大規模並列計算が一般的になってきている。GPUを利用

染物質（NAPL）の除去、地熱開発等で重要となる。ここでは、２相した計 算では、並 列 計 算 用アーキテクチャであるCUDAを利 用す

流シミュレーションを、CCSの諸問題に適用した事例を紹介する。

る。過去数年間に、単相のLBM 計算をGPUに実装した研究が実施

CCSプロジェクトは、CO 2 の大 規 模 排 出 源で CO 2 を回 収し、地され、CPUを用いた場合に比べて、パフォーマンスが大きく向上す

下 深 部の貯 留 層（ 岩 石 ）の 間 隙にCO 2 を注 入することで、永 久 的[1]ることが示されている [8,9,10] 。

にCO 2 を地下に閉じ込める技術である 。CCSは、大気中のCO 2 を

本研 究では、2 相のLBMの計算コードをGPUに実 装し、三次 元多削減できる近未来的な技術として注目されている。貯留層に圧入

孔 質 媒 質（天 然 岩 石のデジタルモデル；図 1）内の多相 流 挙 動を調（ 貯留）したCO 2 の挙動は、多孔質媒質内の 2 相流として扱われる。

べた。デジタル岩石モデルの構築は、高解 像 度のマイクロ CT スキャ

CO 2 の圧入は、非濡れ性流体のCO 2 が、岩石間隙内の水（濡れ性流ナを用いて岩石サンプルの間 隙 構 造の断面 画 像（ 2D ）を取 得し、補体）を押し出す排水過程（Drainage process）と考えることができ

間することによって解 像 度を落とさずにバイナリー三次 元 計算 格 子る。一方で、CO 2 圧入 が 終了した後は、地層水がCO 2 を押し出す浸を生 成した。この巨大 岩石モデルに対して、GPUスーパーコンピュー

潤過程（ Imbibition process）を考える必要がある。このCO 2 貯留

ターを用いた 2 相のLBMシミュレーションを実 施した。

において、間隙スケールでみられるHaines jumpや、スナップオフ、
フィンガリング現象は、貯留層内でのCO 2 の安定度、浸透特性、飽和 度などに関 係する。そのため、間 隙スケールのCO 2 の挙 動を評価することは、CO 2 の貯留能力、漏洩リスクの評価、効率的なCO 2貯留を実現する上で、不可欠であると考えられている。
近年、
格 子ボルツマン 法（ Lattice Boltzmann Method; LBM）が、
多 孔 質内の流 体 挙 動を計算する手法として注目されている [2] 。この手法は、一 般的な数値 力学的手法とは異なり、流 体を衝突と移動を繰り返す多 数の粒子の集合 体と捉え、粒子の動きを計算することによって流 体の運 動を求める数値 計算法である。
既 存の多 相 流モデルとして、Level set 法 [3] 、Volume of fluid（ VOF ）法 [4] 、Phase-field 法 [5] などがあるが、一 般に界 面を取り扱うために複 雑なアルゴリズムを必 要とする 。一 方、LBM では界 面を分 離することが容 易であり、また界 面を操 作する必 要もない 。
さらに LBM では、複 雑な境 界を取り扱うことができる 。このような特 徴から、LBM は複 雑な間 隙 形 状を持つ天 然 岩 石 内の多 相 流
挙 動を計 算する上で、適 当なアプローチといえる 。

02

図1

マイクロX-CT 画像から構築されたべレア砂岩のデジタルロック。黒い部分が間隙であり、CCSプ
ロジェクトでは、この間隙にCO 2 を貯留する。

格子ボルツマン法と GPU 実装

2

大 規 模 計 算を行う際には、メモリを節 約するため、計 算する必要のない固 体メッシュの情 報を保 存しない方 法を用いた。さらにストリームパタンを工 夫し、一つだけの配 列に分 布 関 数を保 存する方法をとった。
マルチGPU 計 算では、はじめに計 算 領 域を分 割する必 要があ格 子ボルツマン法（ LBM）は近 年、多 孔 質 媒 体 内の流 体 挙 動の有る。今回は、流体格子数を均等に割り、z 方向に分割する方法を用

効な計 算 手 法として注 目されている。今 回の計 算では、RK color

いた。分 割した領 域の境 界では、一 層のGhostレイヤーが設けら2 相の LBMモデル [ 11] を用いた。このモデルでは、MRT（ multiple

れ、これは領 域 間のデータ転 送かつカラー勾 配の計 算に用いられ- relaxation - time）モデルなどの改良を加えることで、2 相流の粘る。GPU 間のデータ転送は、2 つのパターンがある（ 図 3 ）。同一の性比の高い条件でも安定して計算することができる。境界条件は、

マザーボード（ ノード ）に装着しているGPUボード間のデータ転送

Non-slip 境界条件をBounce - backスキームで実装した。濡れ性は、はPeerToPeer 方式で行う。一方で、異なるノード間のデータのや固体格子の擬似カラー関数で再現した。

り取りはMessage Passing Interface（ MPI ）に準じて実施する。デー

LBMには、大規模 GPU 並列計算を行う上で、いくつかの利点がタ転 送を効 率するために、CPUホスト側のメモリ保 存は、Pinnedある。まずPoisson 方程式を解かない陽解法であること、かつ局所メモリ方式で行った。このようなコードの最適化によって、2 相流的な計 算がメインであることである。また界 面の追 跡にフェーズ

シミュレーションの計算効率を大幅に向上させることに成功した。

フィールド法の概念を取り込み、並列化した計算が容易である。

以上のようにLBMの計算コードを効率化し、さらにTSUBAME 2.5

マルチGPU の実 装は、CUDAコアレベルと、GPU 間のレベルの

ス パ コ ンの 20ノード（ 60 枚 Tesla-K20 GPU）を利用することで、

2 段 階に分けて実 施した。シングルGPUチップ内 部の並 列 計 算ア1000 3 グリッドという大きな領域で CO 2 圧入シミュレーションを実

ルゴリズム（CUDAコアレベル）は、マルチスレッド処理と呼ばれ現した（図 4 ）。このような巨 大グリッドモデルの計 算を行うことる。CUDAはソフトウェア上で、大量のコアを、スレッドとブロッができれば、長さ数ミリメートルの岩 石 内 部の CO 2 の挙 動を調べ

クの構造で管理している。並列演算が可能な最少単位をスレッド、ることができるようになる。このように計 算 領 域のサイズ アップ

スレッドの集合体をブロックと呼ぶ（図 2）。計算では、各スレッドを行うことで、間隙スケール（μmスケール ）から岩石スケール（ mm

が担 当する格 子を指 定することが重 要である。すなわち、スレッ

スケール）の数値シミュレーションを実 施することが可能となり、ドIDと、更 新される格 子の空 間 座 標を示さなければならない。本

実験室の結果と比較した議論を行うことができるようになった。

研 究では、図 2 a に示すように、z 方 向の一 列の格 子 点の更 新を 1
スレッドで行うこととした。1ブロックが担 当する計 算 領 域は、
blockDimx, blockDimyのパラメータによって決められる。図２aに示すように、ブロックをx 、y 平面にマッピングした。従って、x 、
y 方向のブロック数は、それぞれNx/blockDimx , Ny/blockDimyである（Nx 、Nyは計算領域 x 、y 方向のメッシュ数）。これにより、関
数 内のx 、yループが必 要なくなり、GPUスケジューラが自 動 的に複数のコアを利用して、並列計算を行うことが可能となる。スレッ
ド、ブロックID（ threadID.x, threadID.y, blockID.x, blockID.y）と格子座標（ jx, jy）の関係は、次式で表される。

GPUメモリ内に確 保されている配 列のデータレイアウトも、計算を効 率 的に行う上で重要となる。GPUのメモリバンド幅を無 駄
なく活 用するため、コアレッシングとよばれるデータまとめて転送する方 式の利 用が重 要となる。ここでLBMの分 布 関 数は、以 下のレイアウトで保存される。

またデータアクセス量を減らすため、衝 突 過 程と移 流 過 程を統 合した。これにより、計算効率を、30％も向上することに成功した。

図2

GPU 内部のマルチスレッド化 。

03

大 規 模 GPUコンピューティングによる
天 然 岩 石 内の二 相 流シミュレーション

図3

マルチGPU 間のデータ転送モデル。
黒い四角は、Ghostレイヤーを示している。

適用例

3

図4

2 相のLBMシミュレーション結果（グリッドサ
イズ：1000 3 ）。白い部分が間隙内に圧入した
2 非 濡れ性 流 体 ）である 。ここでは、岩 石
CO（構成粒子と水は透明にしている。右側から左側に向かってCO 2 を圧入している。

挙 動を推 定できている。これらの結 果から、今 回 開 発した 2 相のLBMシミュレーションを用いれば、様々な貯 留 層の条 件における
CO 2 の飽和度や浸透率を定量的に推定できることが示された。

2 相のLBMシミュレーションを使って、CO 2 地中貯留における課題を扱った研究成果を紹介する。ここでは（ 1）貯留層条件が、CO 2 水の挙動に与える影響を調べた研究、（ 2）CO 2 を圧入した岩石に再び水を圧入する湿潤過程を計算して、効果的にCO 2 残留トラップ
や溶解トラップを引き起こす貯留層条件を調べた研究 、（3）
CO 2 の挙動から鉱物化をモデリングし、その鉱物化が流体挙動や浸透率に与える影響を明らかにした研究を紹介する。
３．１

様々な貯留層条件での多相流の挙動

貯留層内の流体挙動をコントロールするパラメータ（例えば温度、界面張力、間隙形状、濡れ性、水圧など）は多い。ここではデジタル岩 石を用いたアプローチにより、貯 留 層 条 件に依 存した流 体 挙
動 特 性（浸 透 率など）を定 量 的に評 価した研 究を紹 介する（図 5）[12,13]。

界 面 張 力を変 化させてCO 2 の挙 動などを調べたところ、界 面 張
力の小さい場合（ 図5a ）のほうが CO 2 は速く流れることが分かった。
つまり、浸透率が高く、圧入効率が良いことになる。実際に、CO 2の挙動から相対浸透率を計算すると、界面張力が増加するにつれて、浸 透 率が減 少することが明らかになった（図 5c）。特に、非 濡
れ性流 体であるCO 2 の浸透率ほうが、界面張力の影響を強く受けることが分かった。一方で、界面張力が大きい場合（図 5b）のほうが、最終的なCO 2 の飽和度は大きくなることが分かった。つまり、
界 面 張 力が大きいほうが、限られた間 隙スペースを有 効に使ってCO 2 を貯 留できることになる。ここでは、界 面 張 力を変 化させた例（図 5）だけを紹介したが、これ以外にも様々な条件下で、CO 2 の04

図5

界面張力（ IFT ）に依存したCO 2 の挙動 [12] 。
（a）低い界面張力の場合。
（b）高い界面張力の場合。
（c）水とCO 2 の飽和度が 50%の場合の相対浸透率と界面張力の関係

３．２

残留・溶解トラップが効果的に働く貯留層条件の推定

貯 留 層 内での CO 2 の残 留トラップ量を定 量 的に調べるため、LBM

鉱 物 化の計 算 結 果から、鉱 物 化が生じる場 所は、流 速だけに依存しているのではなく、不 均 質な間 隙 形 状にも依 存していることが明らかとなった（ 図 7）。鉱物化に伴う浸透率の変化を計算する

シミュレ ーションを用いて界 面 張 力が残 留 CO 2 クラスターに与えと、絶対浸透率は鉱物化によって大きく影響を受けること、CO 2 のる影響を調べた。CO 2 の圧入（ 図 5）が終了した後に、水を圧入する

相 対 浸 透 率は大きく低 下することが分かった。一 方、水の相 対 浸

シミュレーションを実施し、間隙内に残留するCO 2 の量と分布を調

透率は、あまり変化しないことが分かった。これは鉱物化によって、べた（ 図 6 ）[14] 。この残留 CO 2 の量が多いほど、安定的に岩石間隙

間隙（ Pore throat）の径が小さくなり、そこを通過するのに要する内にCO 2 を貯 留できていることになる。また残 留 CO 2 クラスター

Capillary pressureが増加する。それにより、CO 2 といった非濡れのサイズ、長さ、接触面積の分布を計算し、界面張力がこれらの性性 流 体の移 動が制 限され、浸 透 率が大きな影 響を受けると考えら質に対する影響について検討を行った。

れる [16] 。

計算の結果、界面張力が増えるにつれて、残留 CO 2 クラスターの

今 回、デジタル岩 石 内で生じる鉱 物 化や、それが水 理 特 性に与サイズが大きくなり、トラップ量が増 加する傾 向があることが分える影 響を調べることができた。この鉱 物モデリングで得られる

かった（ 図 6c ）。一方で、低い界面張力では、クラスターが小さく、情報を、一般的に実施されるフィールド（km）
スケールの貯留層シ

CO 2 は広い表 面 積を持つことが分かった（ 図 6a）。つまり、この条

ミュレーションに組み込めば、鉱物反応を考慮して、長期的なCO 2

件は、CO 2 の地層水への溶解（溶解トラップ）にとって有利であるとの挙動予測が可能となる。

考えることができる。この方 法によって、異なった貯 留 層 条 件における残 留 CO 2 の状 態を予 測することができる。つまり、安 定したCO 2 地 中 貯 留を実 施する上で、適 当な貯 留 層 条 件を明らかにすることができる。

図7

図6

３．３

残 留した CO 2 の分 布 [14]。左 側から右 側 へ
界 面張力が高くなる。ここでは、岩石構成
粒子と水は透明にしている。

間隙流体挙動を考慮した鉱物化シミュレーションの結 果 [15] 。赤 色 部 はCO 2 が化した炭 酸 塩 鉱 物を示し、白 色 部は間 隙を示す 。左 側は鉱 物 化 前の間
隙 形 状である 。

CO 2 の鉱物化モデリング

圧 入したCO 2 は、最 終 的に鉱 物 化すると考えられている。ここではCO 2 の鉱 物 化をモデル化し、それが間 隙 流 体 挙 動（例えば浸 透率）に与える影響を調べた [15] 。まずLBMで計算された流速場から、
移 流 拡 散 方 程 式を用いることで、各グリッドでの CO 2 濃 度を計 算した。その濃度が、ある閾値を超えた場合に、そのグリッドを鉱物（ 炭酸塩鉱物）に変化させた。その後、鉱物化した岩石モデルに対して、再度、LBMの流体シミュレーションを実施した。この操作を各時間ステップで行い、鉱物化の進行を再現した（図 7）。なお、この鉱 物 沈 殿の計 算 結 果を、実 験 結 果と比 較することよって、鉱 物化モデルに用いるパラメータを正確に決定することができた。

05

大 規 模 GPUコンピューティングによる
天 然 岩 石 内の二 相 流シミュレーション

まとめ

4

[5]

B a d a l a s s i , V. E . , H . D . C e n i c e r o s , S . B a n e r j e e ( 2 0 0 3 ) ,Computation of multiphase systems with phase fieldmodels. J. Comput. Phys. 190(2), 371–397.

[6]

Ahrenholz. B, J. Tölke, P. Lehmann, A. Peters, A. Kaestner, M.
Krafczyk, W. Durner (2008), Prediction of capillary hysteresis

GPU 大 規 模 計 算により、2 相の格 子ボルツマン法を巨 大なデジタ

in a porous material using lattice-Boltzmann methods andル岩石モデルに適用し、間 隙 内の流 体 挙 動を計 算することに成 功

comparison to experimental data and a morphological

した。今 回のようなLBMによる流 体 挙 動シミュレーションでは、

pore network model, Advances in Water Resources, 31 (9),1151–1173.

マルチGPUの利用が効果的である。
今 回 開 発したアプローチを、C O 2 貯 留プロジェクトの課 題に適[7]

Ramstad. T, N. Idowu, C. Nardi, P.E. Øren (2012), Relative

用した結果、次のことが明らかになった。（1）貯留層の状態が、多permeability calculations from two-phase flow simulations

相流の挙動に与える影響を明らかになった。
（2）CCSプロジェクト

directly on digital images of porous rocks, Transport in
Porous Media, 94 (2), 487-504.

において残 留 CO 2 トラップや溶解 CO 2 トラップが効果的に働く条件を明らかにすることができた。
（ 3）C O 2 の 間 隙 内での挙 動から[8]

Tölke, J., M. Krafczyk (2008), TeraFLOP computing on a

C O 2 の鉱 物 化をモデル化し、その鉱 物 化が浸 透 率に与える影 響を

desktop PC with GPUs for 3D CFD. Int. J. Comput. Fluid

調べることができた。また、ここでは紹介しないが、本研究で用いDyn. 22(7), 443–456.

たデジタル岩 石の弾 性定数を、弾性波動シミュレーションで推 定することも可能となってきている [17] 。

[9]Tölke, J. (2010), Implementation of a lattice Boltzmann
kernel using the compute unified device architecture

このようにデジタル岩 石を用いた大 規 模 計 算によって、貯 留 層

developed by nVIDIA. Comput. Vis. Sci. 13(1), 29–39.

内での C O 2 の挙動や反応を再 現することが可 能となってきた。近[10] Kuznik, F., C. Obrecht, G. Rusaouen, J.J. Roux (2010), LBM

い将 来、CO 2 地 中 貯 留をはじめとした地 下 流 体を扱うプロジェク

based flow simulation using GPU computing processor.

トでは、貯 留 層 全 体をデジタル化し、貯 留 層 内で生じている挙 動

Comput. Math. Appl. 59(7), 2380–2392.

や反 応を正 確にモデル化することが可 能になると考えられる。そ[11] Tölke, J., S. Freudiger, M. Krafczyk (2006), An adaptive

れにより、貯 留 層を有 効に利 用または正 確にコントロールできる

scheme using hierarchical grids for lattice Boltzmann multi-

ようになると考えられる。

phase flow simulations. Comput. Fluid 35, 820–830.
[12] Jiang, F., T. Tsuji, and C. Hu (2014), Elucidating the role

謝辞

of interfacial tension for hydrological properties of two-

本 研 究は、日 本 学 術 振 興 会（ JSPS ）の新 学 術 領 域 公 募 研 究（ no.

phase flow in natural sandstone by an improved lattice

15H01143）、基 盤 研 究 A（no. 24246148 ）、スタ ートア ップ（ no.

Boltzmann method, Transport in Porous Media, 104, 1, 205-

26887028）、エジプトMOSR と の 二 国 間 共 同 研 究、JICA/JST の229.

SATREPSプロジェクト、文 部 科 学 省の世 界トップレベル研 究 拠 点（WPI ）プログラムの助成を受けて実施したものである。

[13] Yamabe, H., T. Tsuji, Y. Liang, and T. Matsuoka (2015),
Lattice Boltzmann simulations of supercritical CO 2 - waterdrainage displacement in porous media: CO 2 saturation

参考文献

and displacement mechanism, Environmental Science &[1]

Technology, 49 (1), 537-543.

Metz, B., et al. eds. (2005), IPCC Special Report, Carbon
Dioxide Capture and Storage. Cambridge University Press,[2]

residual CO 2 clusters in porous sandstone, Water Resources

Boek, E.S., M. Venturoli (2010), Lattice -Boltzmann studies

Research, 51, 1710-1722.

of fluid flow in porous media with realistic rock geometries.
[3]

[4]

06

[14] Jiang, F., and T. Tsuji (2015), Impact of interfacial tension on

Cambridge.

[15] Jiang, F., and T. Tsuji (2014), Changes in pore geometry andComput. Math. Appl. 59, 2305–2314.

relative permeability caused by carbonate precipitation in

Sussman, M., A.S. Almgren, J.B. Bell, P. Colella, L.H. Howell,

porous media, Physical Review E 90, 053306.

M.L. Welcome (1999), An adaptive level set approach for[16] Kitamura, K., F. Jiang, A.J. Valocchi, S. Chiyonobu, T. Tsuji,

incompressible two-phase flows. J. Comput. Phys. 148,

and K.T. Christensen (2014), The study of heterogeneous

81–124.

two-phase flow around small-scale heterogeneity in

Hirt, C.W., B.D. Nichols (1981), Volume of fluid ( VOF)

porous sandstone by measured elastic wave velocities andmethod for the dynamics of free boundaries. J. Comput.

lattice Boltzmann method simulation, J. Geophys. Res. 119,

Phys. 39(1), 201–225.

7564-7577.

[17] Yamabe, H., T. Tsuji, Y. Liang, T. Matsuoka (2016), Influenceof fluid displacement patterns on seismic velocity during
supercritical CO 2 injection: Simulation study for evaluationof the relationship between seismic velocity and CO 2
saturation, International Journal of Greenhouse Gas
Control, 46, 197-204.

07

大規模画像データセットの
機械学習のための分散コンピューティング
佐 藤 育郎 *渡邉 隆太郎 ** 株 式 会 社デンソーアイティーラボラトリ

西村 裕 紀 **** 株 式 会 社デンソー

野村 哲弘 ***

松 岡 聡 ****** 東 京 工 業 大 学

ディープニューラルネットワーク（DNN）と称される機 械 学 習手 法を使って、大 規 模データセットの教 師あり学 習を行うことで、非常に高い識 別 能 力を獲 得できることが、最 近の研 究によって明らかとなってきた 。本研 究は、大 規 模データセットにおいて、複 数の DNNモデルを学 習・評 価することで性 能の良いモデルの条 件を定 量 化し、また優 良な複 数のモデルをアンサンブル化することで、高い識 別 能 力を持つ画 像 識 別 器を生 成することを目的としたものである。学 習には独自開 発した非 同 期 型
分 散学習プログラムを使 用し、1モデルあたり最多で 96GPUを利用して高速な学習を行った。学習はTSUBAME グランドチャ
レンジを利 用して行われ、ピーク時で 1146 台の GPUを使い、微 分 計算 部分においては 1GPUあたり約 1TFLOPS の単 精度 浮
動 小 数 計算性 能を引き出した 。

1

はじめに

ソフトウエア2

画 像 識 別は情 報 技 術 分 野における主 要なテーマのひとつである。

DNNの学習アルゴリズムの分散には、一般的にモデル並列とデー

産 業 応 用の例としては、監 視カメラにおける人の識 別、生 産 設 備

タ並列の二種の方針がある。モデル並列とは、1 つのモデルを複数における製品の外観検査、自動車の自動緊急ブレーキ装置におけのGPU（ここでの議 論はGPUの利 用を前 提とする）に分 散させ必る障 害 物 検 知などが挙げられる。特に人の安 全にかかわる分 野に要に応じて変 数を相 互に通 信し同 期を取って処 理することで、コ

おいては、極めて高い識別精度が要求されるため、技術的な革新がスト関数の微分計算を行うものである [1, 4, 5] 。微分計算に必要とさ求められている。

れる全 変 数の容 量が 1 台のGPUのメモリを超える場 合、この方 針

一 般 的に、機 械 学 習を用いた識 別 器の性 能は、使 用するアルゴが取られることになる。つまり比 較 的 大きなネットワークを使 用

リズム、データセットのサイズ、計算機の能力に依存する。アルゴする際に不可欠となる並列化方法である。もう一方のデータ並列

リズムにおいては、あらかじめ決められた手 順で画 像 特 徴を抽 出は、一 度の微 分 計 算に使 用する複 数のデータを、複 数のGPUに分する従 来 的アプローチよりも、学 習によってその手 順を決 定する

散させ処 理するものである [4, 5] 。プログラムにおいては、モデル並

DNNの方が、物体識別のタスクにおいて圧倒的に高い識別能力を

列のみ [1] 、または本 研 究のようにデータ並 列のみを実 現すること有することが、Krizhevskyらによって実証された [1] 。このブレークも可能であるし、またその両方を実現することも可能である [4, 5] 。

[2]

スルーの背 景にデータセットの大 規 模 化がある ことは注 目に値

データ並列には、さらに同期型 [6] のものと、非同期型 [4] のものがする。事 実、小さなデータセットではDNNは従 来 的アプローチに有り得る。同 期 型の場 合、全ワーカー GPUが、それぞれ画 像を読[3]

性 能で劣ることがあることが報 告されている 。DNNの最 適 化にみ込み、その画像に基づいた微分計算を終えたのち、通信によっては多 大な演 算 量が要 求されるため、GPUの利 用が現 在では標 準 的微分を加算し、重みパラメタを更新し、最新の重みを全ワーカーにである。それでもなお、大規模データセットの学習に、単体のGPU

配 布するサイクルを反 復することとなる。利 点は、単 体のGPUでで 1 ヶ月またはそれ以 上の期 間を費やすことはそれほど珍しいこ行う学習に比べて単位時間あたりに処理できる画像枚数が大きくとではない。

出来る点と、収束の保証がある点にある。一方、非同期型の場合、こうした背 景から、我々は自 動 車の走 行 安 全への応 用に向け

計 算 機のリソースを間 断なく利 用するため、同 期 型に比べてさらて、ノード分 散 型の DNN 学 習プログラムを開 発し、単 体の計 算に反復計算の速度を向上させることが可能である利点がある一方、機に比 較して数 十 倍の速 度で反 復 計 算を実 行 可 能であることを「古い」重みパラメタを使った微 分を使うことから理 論 的に収 束

TSUBAME2.5 上で確認した。

の保 証がないといった欠 点が挙げられる。もっとも、後 者の欠 点

本 研 究では、ノード分 散 型の学 習プログラムを大 規 模データには実質的には問題を起こさないというという指摘もある [4] 。

適用し 、（1）DNNのモデル選択とアンサンブル化による高精度画

データ並列には、同期／非同期に関わらず、さらに、重みパラメ

像識別器を、TSUBAME2.5 で比較的短期間に生成した。また、（2）タの更 新のための「パラメタサーバ」と呼ばれる特 殊ノードを置

モデルの学習パラメタ数と識別性能との関係を定量的に明らかにく方 式 [4] と、パラメタサーバの機 能をワーカノードに分 散 化させした。
（ 2）については現 段 階では予 備 的 結 果であるものの、パラた、特 殊ノードを置かない方 式 [5, 6] の 2 種 類がある。これ以上の説

メタ数を広 範 囲に分 布させた結 果、パラメタ数と識 別 性 能の間に明は、本文献の主旨とは若干隔たりがあるため、割愛する。

ある興味深いトレンドが見えてきた。

本 研 究では、特 殊ノードを置かない方 式の非 同 期 型データ並 列のみを実 現した、独 自 開 発のプログラムを使 用した。分 散アルゴ08

リズムの詳 細と学 習 速 度の定 量 評 価については、別の機 会で論ず

めに、全結合層にdropout [10] などの正 則 手 段を適 用することがよる予定である。

く見 受けられる（ 例えば [1, 8] ）。本 検 討では全 結 合 層の使 用は最 後の一層のみとし、また、最後のマップサイズは、フィルタサイズと同じ 3 × 3 とおいた。このようにして、全結合に必要なパラメタ数

データセット3を抑えた。我々の行った予備実験において、表1のようなネットワー
クは、複数の全結合層を持つ似たようなネットワークに比べて、パラメタが少ないにも関わらず高い汎 化 性 能を持つことを確 認している。
最 適 化の詳 細については表 2 に記 載する。分 散 計 算には 1モデ

本 研 究で使 用したILSVRC2012データセット [2] は、ImageNetと呼ルあたり最少で 24GPU 、最多で 96GPUを使用した。本研究で行っばれるデータセットのサブセットであり、128 万 学 習サンプルと 5たモデルの多くは、TSUBAMEのグランドチャレンジを利用して学

万検証サンプルと 10 万枚のテスト画像を含んでいる。これは、公習したものである。モデルあたりの分 散 学 習 日 数はおおむね 7 日開の画像識別データセットとしては最大級のもののひとつである。

である。グランドチャレンジでは、ピーク時で 1146 台のGPUを使ここでいうサンプルとは、画像と正解クラスの組みを指している。

い、13モデルを同時に学習した。特定のDNNの学習部分（ 順伝搬

クラス数は 1000 である。なおクラスあたりの画 像 枚 数は一 律でと逆伝搬）において、1GPU（ NVIDIA Tesla K20X , 単精度浮動小数はない。テスト画 像に対する正 解クラスは開 示されていない。テ

計 算 理 論 性 能 3.95TFLOPS ）あたり 0.98TFLOPSの単 精 度 浮 動 小

ストサーバに識 別 結 果を投 稿することで、テストセット全 体の誤数計算性能を実測した。

識別率が公表される仕組みとなっている。

DNN モデル4

本 検 討におけるモデル設 計の指 針について述べる。パラメタ数と識 別 性 能の関 連を明らかにしたいため、パラメタ数は学 習サンプ
ル数の 1 倍から 30 倍程度の広い範囲で設定した。活性化関数による非線形写像を含む畳み込み層を連続的に反復させることで、ネッ
トワークの表 現 能 力を高められるという近 年の知 見 [7-9] に従って、畳み込み層は 11 層から 48 層の範 囲で設 定した。畳み込みのフィ
ルタサイズは一律 3ピクセル四方とした。各層の特徴マップ（畳み込みなどの操作によって得られる 2 次元信号のこと）の数は、畳み表1

本検討で使用した、畳み込み層 15 層、プーリング層 5 層 、パ ラ メタ数 16M 個 を持つ
ネットワーク構造。
C は畳み込みを、CPは畳み込みとプーリングを、
S はソフトマックス関数を意味する。活性化の記載は省略している。層の番号は、パラメタを持つ層に対して与えられている。

込み層の層 数とパラメタ数の設 計を元に、手で設 定した。本 検 討においては、全 結 合 層は最 終 隠れ層と出 力 層の間の一 層のみとした。入 力 画 像サイズは、目 視による検 査の結 果 200ピクセル以 上に設定することとした。これは、特徴マップの解像度を下げるプー
リング層を畳み込み層と畳み込み層の間に適 当に挟むことで実 現できる。プーリング層の層数は 2 から 5 の範囲である。参考までに、
本 検 討で使 用したモデルをひとつ取り上げ、そのネットワーク構造を表 1 に記載する。
全結合層について補足する。通常、畳み込みネットワークでは、マップの解 像 度を十 分に低 下させたのち、全ニューロンを使った

表 2 最適化の詳細

線形結合の層（ 全結合層 ）を複数層（ 3 層程度 ）設けるのが標準的である（ 例えば [1, 8] ）。これらの層では、ニューロン毎に特定のパラメタが学習されるため、多量のパラメタを必要とし、結果として過学 習を生じさせやすくする傾 向がある。こうした過 学 習を防ぐた09

大 規 模 画 像データセットの機 械 学 習のための分 散コンピューティング

5

定性評価結果中央にいるパフォーマーのポーズなどから、注目物体の種別を推測した結果であると思われる。つまり、DNNが画像全体から、シーンのコンテキストを読み取れていることが分かる。同様の分析は図 1（b）の画像（ 正解ラベルは「パドル」）にも当てはまる。図 1（c ）はシャボ ン玉の識 別に成 功した事 例を示している。従 来、半 透 明 物

適切に識別できている検証用画像の例を図 1 に提示する。図 1（ a）

体の識 別は難しいとされてきた。透けた部 分にある背 景の画 像 特の正 解ラベルは「一 輪 車」である。このサンプルでは、注 目 物 体

徴が混入してしまうことで、識別対象の特徴が「汚染」されてしまである一 輪 車の車 輪 部 分が画 像 外にあり、シャフトがかろうじてうことが識別を困難にしていると考えられる。しかしながら我々の見える程 度という点に難しさがある。しかしながらアルゴリズム

モデルはシャボン玉を高スコアで識別できており、こうした困難さ

は高い確度で識別に成功できている。これは、上を見上げる群集と、をDNNが克服していることが窺える。

図1

識別に成功できたImageNet検証データの例。

（b）paddle（a）unicycle（c）bubble

10

識別に失敗している検証用画像の例を図 2 に提示する。図 2（a）
では、「アラビアンキャメル」が正 解であるところを、「犬ゾリ」と誤 答している。これは、DNNがコンテキストを過 剰に評 価した

6

定量評価結果結 果、画 像にうつるアラビアンキャメルの頭 部の特 徴が過 少 評 価されたものと推 測される。図 2（ b）は「トリケラトプス」が正 解であるが、我々のモデルはこれを識 別できていない。正 面 方 向か

識別性能とパラメタ数の関係を図 3 に記載する。これによると、パら見たトリケラトプスの形 状が識 別できていないことから、汎 化

ラメタ数が 8 ×10 6 〜 40 ×10 6 の範 囲にあるとき識 別 性 能は同 等に限界があることが窺える。

程度に良好だが、パラメタ数が 8 ×10 6 を下回るにつれ識別性能が急 速に劣 化する傾 向が確 認できる。使 用したモデルのネットワー
ク構造は様々であるにも関わらず、このような傾向が確認できたということは、識 別 性 能においてパラメタ数の影 響が支 配 的である
ことを示唆している。分散計算環境における計算速度を考慮に入れると、パラメタ数の少ないモデルの使 用が好ましい。このため、

図2識別に失敗したImageNet検証データの例。

今回の検討からは、8 ×10 6 個程度のパラメタ数を持ったモデルの適用が理に適うこととなる。

（a）Arabian camel

図3

識別性能とパラメタ数の関係。
青線は学習サンプルの個数を示す。

この結果が、機械学習の古典的な考えとはいくぶん乖離した結果であることは指摘に値する。古典的には、モデルのパラメタ数は学
習サンプル数を下回るよう設定することが、過学習による性能劣化を防ぐために重要であると考えられている。しかしながら今回の結果を見る限り、パラメタ数が学習サンプル数を下回ることで良い識別性能が得られるとは推論しがたく、学習サンプル数の６倍程度またはそれ以上のパラメタ数を持つモデルが良好であることが確認できる。DNNがこれほど多数のパラメタ数を必要とする原因はまだ明らかにできていないため、今後の調査課題としたい。
DNNモデルを網羅的に組み合わせて生成したアンサンブル識別器のうち、表 3 に記載した 6モデルからなるものが、ImageNetの検証（b）triceratops

用データにおいて優良な識別性能を示した。ここでのアンサンブル
識別器とは、ある入力に対するソフトマックス出力の対数の期待値

11

大 規 模 画 像データセットの機 械 学 習のための分 散コンピューティングの大きさに従ってクラス分類を決定するものである。アンサンブル化によりおよそ 2 %の性能向上が確認できた。これは、各モデルがそれぞれ少しずつ違った特徴を抽出していることを意味している。

今後

7

以 前は困 難とされてきた、深 いネットワークの学 習を効 率 的に行
うための知 見が積み重なってきたため、今 後はこれらの知 見を活用して、パラメタ数を低く抑えつつ ネットワークを多 層 化する検討を進めたい。パラメタ数の低 減は学 習の高 速 化にとって極めて重要である。現状のDNNに必要とされるパラメタ数は学習サンプ
ル数を大きく上 回っており、パラメタ数を低 減させる余 地はまだ多く残されていると推 測できるため、今 後の研 究 課 題としたい 。
表3

本研究で使用したアンサンブル識別器の構成モデル。

代表的な手法との比較結果について述べる。テストサーバによる識別性能の評価結果を表 4 に記載する。予測された上位 5 個のクラスの中に正解が含まれない割合を評価指標としている。現在時点の最良のスコアは 3.57% [9] であり、これは人で試験した結果 [2] を上回るものである。手法の特色は、複数の層をまたいだ結合を随所に持つネットワークを利用して、下層の信 号との残 差に対する学習を行う点にあり、この工夫によって、これまででは困難とされていた 100
層を超えるような非常に深いネットワークであっても効率的に学習が行えるようになっている。GoogLeNet [7] は、Inception Moduleと呼ばれる異なるサイズの畳み込み演算を並列的に行う層を持つことが特色で、これにより従来に比べ少ないパラメタ数で高いスコアを獲得することに成功している。OxfordNet [8] は、畳み込み層を多層に謝辞
本 研 究 で 使 用 し た、マルチノ ー ド 対 応 の 非 同 期 分 散 学 習 プロ
グラム は、TSUBAME 産 業 利 用 トライアルユ ー ス 制 度 に よ る
TSUBAME2.5 の無 償 利 用を通じて開 発したものである。また、本研 究 成 果は平 成 27 年 度 秋 季 TSUBAMEグランドチャレンジを利用したものである。ここに感謝の意を表する。
参考文献[1]		A. Krizhevsky, et al., “ImageNet Classification with Deep
Convolutional Neural Networks”, NIPS 2012.
[2]		O . R u s s a k o v s k y, e t a l . , “ I m a g e N e t L a r g e S c a l e V i s u a l
Recognition Challenge”, IJCV 2015.
[3]		M. Oquab, et al., "Learning and transferring mid-level image
representions using convolutional neural networks", CVPR

渡って重ねた比較的シンプルな構造となっており、多量のパラメタにより表現力を高める一方で、dropoutによる正則を採用し過学習を抑えている。AlexNet [1] は 2012 年に性能におけるブレークスルーを起こしたモデルとして有名である。本研究で使用したモデルのパラメタ数はGoogLeNetに次いで少なく、層数はAlexNetに次いで少ない。スコアはAlexNetより良好であり、概ね妥当な範囲に収まっていると思われる。

2014.
[4]		J. Dean, et al., “Large Scale Distributed Deep Networks”,
NIPS 2012.
[5] 		R. Wu, et al., “Deep Image: Scaling up Image Recognition”,
arxiv:1501.02876, 2015.
[6]		F. N. Iandola, et al., “FireCaffe: Near-Linear Accelerationof Deep Neural Network Training on Compute Clusters”,
arxiv:1511.00175, 2015.
[7]		C. Szegedy, et al., “Going Deeper with Convolutions”, CVPR
2015.
[8]

K. Simonyan and A. Zisserman, “Very Deep Convolutional[9]

K. He, et al., “Deep Residual Learning for Image Recognition”,

Networks for Large-Scale Image Recognition”, ICLR 2015.
arxiv:1512.03385, 2015.
表4

DNNを使った手法の識別性能の比較評価結果。
内積を含む層の数とパラメタ数は、アンサンブ
ルに使用したモデルのうち最大のものとした。
（* 開示されているネットワーク構造からパラメ
タ数を計算）[10] G. E. Hinton, et al., “Improving neural networks by preventing
co-adaptation of feature detectors”, arxiv:1207.0580, 2012.
[11] V. Nair and G. E. Hinton, “Rectified Linear Units Improve
Restricted Boltzmann Machines”, ICML 2010.
[12] I. Sato, et al., “APAC: Augmented PAttern Classification with
Neural Networks”, arxiv:1505.03229, 2015.

12

マルチGPU
超並列クラスタシステムを用いた大規模ナノ炭素 分子の電子状態計算
河 東田 道 夫 *成 瀬 彰 **

中嶋 隆 人**（国）理 化 学 研 究 所 計 算 科 学 研 究 機 構 量 子 系 分 子 科 学 研 究チーム** NVIDIA Japan本 研 究では、最 近 新 規にマルチGPGPU を施した RI-MP2エネルギー計 算プログラムを用い、大 規 模ナノ炭 素 分 子集合 体の電子状態計算をTSUBAME 2.5 で実施した。本稿ではTSUBAME 2.5 でのマルチGPUジョブの性能評価結果を報告する。
特に、TSUBAME 2.5全ノード規模で実施した最大規模のナノ炭素分子集合体の計算（ナノグラフェン2量体（C 96 H 24 ）2 , 1,349
ノード・4,047 GPU 使 用）の RI-MP2 計 算では、515 TFLOPs の実 測ピーク性 能を達 成した 。本 稿ではさらに、π-πスタッ
キングしたナノグラフェン2 量 体に働く分 子 間 相 互 作 用の解 析 結 果も報 告する 。

はじめに

1

ナノ炭素分子集合体は、有機分子機能材料への応用を視野に入れて実 験 的、理 論 的、計 算 科 学 的アプローチによる研 究が盛んに行われている。有 機 分 子 機 能 材 料の分 子 設 計や機 能 評 価を行うた
めには、分子集合体の局所的な構造や安定性の理解が不可欠である。電子状態計算はこれらの機能を評価する上で極めて有効な手法であると広く認 識されている。Møller - Plesset 2 次 摂 動（ MP2）

図1

法 [1] はファンデルワールス相互作用などの弱い分子間相互作用のNTChemに実装されている主要な機能。

記 述に重 要な役 割を果たす電 子 相 関を正しい振る舞いで取り扱うことが可能な手法であり、ナノ炭素分子集合体の安定性の議論の最近の研究にて、著者らはResolution-of-Identity MP2（ RI-MP2）

際には有 用な手 法となる。しかしながら、MP2 計 算は前 段 階に実法 [4,5] に基づく電 子 相 関エネルギー計 算の超 並 列アルゴリズムと施されるHartree-Fock（HF）
計算よりも高コストな計算となるため、

計 算プログラムの開 発を行い、京 8,911ノードを用いて 360 原 子、

計 算が可 能 分 子のサイズに制 約がある。そのため、京コンピュー

9,840 原 子 軌 道 系の大 規 模 分 子の RI-MP2エネルギー計 算に成 功

タやTSUBAME 2.5 などの大規模計算システムを活用し、大規模分

[6]
している。

子のMP2 計算を実施することが求められている。

また、最 近になり、graphics processing unit（ GPU）などの加我々の研究チームでは、京コンピュータなどのペタスケールスー

速 器を搭 載したヘテロジニアスなシステム構 成のスーパーコン

パーコンピュータを活 用し大 規 模 分 子の電 子 状 態の超 並 列 計 算

ピュータが年々増 加しており、これらのシステムに搭 載されていを実現することを目的として、NTChemプログラム [2,3] の開発を進る加 速 器を活 用した大 規 模 電 子 状 態 計 算のニーズも日 増しに高めている。NTChemには様々な電子状態理論や効率的な計算アル

まっている。NTChemの開 発においても、このニーズに応えるたゴリズムが実装されている（ 図 1 ）。代表的な計算手法として、（ 1）

め、プログラムのgeneral-purpose computing on GPU（GPGPU）

HF 法や密 度 汎 関 数 理 論（ DFT）に基づき、分 子の基 底 状 態と励 起

実装を進めている。最近の研究にて、著者らはRI-MP2エネルギー状 態の電 子 状 態を計 算する機能 、（ 2）より高精度な電子相関 理 論

計算の超並列アルゴリズムの改良とGPGPU 実装を行い [7] 、マルチであるMP2 法、結 合クラスタ法、量 子モンテカルロ法に基づく基

ノード・マルチGPUシステムを活用しMP2エネルギー計算を実現

底状態と励起状態の電子状態計算機能、などが挙げられる。

可能とする計算基盤の整備を行った。

これらの計算方法のうち、HF 法、DFT 法、MP2 法のコードについ

本稿では、TSUBAME 2.5 上で上記実装を用い、ナノ炭素分子複ては京コンピュータでの超並列計算のためにMPI/OpenMPハイブ

合 体の RI-MP2エネルギー計 算をマルチノード・マルチGPU 環 境

リッド並 列 化が施されており、京 数 千ノードから数 万ノード規 模

で実 施した結 果を報 告する。2 章ではNTChemプログラムにおけを用いた超並列計算が可能である。

る RI-MP2 計 算のマルチGPGPU 実 装の概 略を説 明する。3 章および４章では、初めにTSUBAME 2.5 でのマルチGPGPUコードの性能
評 価 結 果を報 告する。5 章では、π- πスタッキングしたナノグラ
フェン2 量体層間に働く分子間相互作用解析を行った結果を報告する。

13

マルチGPU 超 並 列クラスタシステムを用いた
大 規 模ナノ炭 素 分 子の電 子 状 態 計 算

NTChem プログラムにおける
RI-MP2 計算のマルチ GPGPU 実装

2

改善では、最外ループとなる仮想軌道インデックス（ a , b ）の第 1 軸のMPI 並 列 化に加え、このループ内で実 行される積 分 計 算と行 列
演 算に関わる行 列をブロック分 割し第 2 軸の並 列 化を行った。このアルゴリズムとコードの改 良により、前 述の 3,360 原 子、9,840
原 子 軌 道 系の大 規 模 分 子の計 算を最 大で京 80,199ノードを用い

２．１

RI-MP2 法

RI-MP2 法 [4,5] は 4 中心電子反発積 分（ electron repulsion integralsて実施することに成功している。[7]２．３RI-MP2 法のマルチGPGPU 実装（ ERIs））のRI 近似（あるいは密度フィッティング近似 ）を適用する

本 研 究を実 施する以 前の研 究において、2.2 章で説 明したRI-MP2

ことでMP2 計算を高速化することが可能な計算手法である。

エネルギー計算プログラムのCUDA [10] を用いたGPGPU 実装を行った [7] 。GPGPU 実装を行うにあたり、行列 - 行列積計算に必要なデー
タをCPUからGPUに転送し、計算をGPUで実施するオフロードモ
デルを採用した（ 図 2 ）。

ここで、式（1）および（2）に対し原子軌道ペアに対して補助基底
関数（インデックス（ l , n ））の線形展開を行う近似を導入している。
この近 似によりMP2 計算のプレファクタを大きく減らすことが可能なため、結 果として計 算コストと使 用メモリ量を大 幅に減らすことが可 能となる。RI-MP2 計 算において、最もボトルネックとなる処理は式（1）の行列 - 行列積計算であるが、この処理は最適化されたBLAS [8] に実装されているDGEMMを用いることにより、非常に高い実 行 効 率およびスレッド並 列 性 能での処 理が可 能である。
MP2 電子相関エネルギーは式（ 1）、（ 2）から計算された近似積分を用いて式（3）より計算される。

図2

RI-MP2 計算のGPGPU 実装の概要

各 GPU へのタスク割り付けは 2.2 章で言 及した行 列 - 行 列 積 演
算の並 列 化に用いている MPIプロセスを各 GPU に割り当てること
で実 現されている 。行 列 - 行 列 積 演 算は NVIDIA により最 適 化さ
２．２

RI-MP2 計算の超並列実装

している 。この GPU 計 算を行う際には CPU - GPU 間のデータ転 送

本研究を実施する以前の研究において、RI-MP2エネルギー計算プ

が処 理の大きなボトルネックとなる 。この 処 理のオーバーヘッド

ログラム の 超 並 列 アルゴリズム 開 発 と NTChemプログラム へ のを減らすために、異なる仮 想 軌 道のインデックス（ a , b ）の複 数のMPI/OpenMPハイブリット並列実装を行った。[6,7] 本実装では、計

行 列 演 算をまとめて実 施するバッチ処 理を実 装した。このバッチ

算のボトルネックとなる（ 1）式の 4 中心電子反発積分計算の並列処理では、同じバッチに属するデータ転送と演算を 1 度にまとめて化を行う際に仮想分子軌道のインデックス（ a , b ）を分散すること

行い、さらに ペー ジロック
（ピ ンド）
メモリを活用することにより、により、MPI 並列化を行っている。この並列化手法により、占有分

大きなオーバーヘッド削 減が可 能である。ページロックメモリは子軌道のインデックス（ i , j ）をMPI 並列化していた従来の並列アル

CPUのホストメモリとのデータ転 送をする際にGPUからのダイレ

ゴリズム [5,9] と比較して大幅な並列スケーリングの向上が可能であ

クトメモリアクセスを実 現することができるため、CPUを使 用せる。さらに、分散メモリの有効活用を目的として、式（1）と（2）にずにデータ転送が可能であり、性能向上に大きく寄与する。

現れる中 間データをインコアメモリに保 持するアルゴリズムを選本研究で使用したCPUコードのコンパイルにはIntelコンパイラ

択した。本 手 法によりMP2 従 来 法で大きな負 荷となっていたI/O

を、GPUコードはCUDA 5.0コンパイラを使用した。また、数値計

処理を除去することが可能となり、性能の向上に寄与している。

算ライブラリとして、Intel Math Kernelに実装されているBLASおまた、最近になり、数万 MPIプロセスを用いた超並列計算を実行可能とするため、局所的なMPI communicator を活用し 2 次元的なMPI 階層並列化を導入したアルゴリズム改善を実施した。この14れた cuBLAS [ 11] ライブラリcublasDgemm ル ー チンを用いて実 装

よびLAPACK 、cuBLASライブラリを使用した。

TSUBAME2.5 における RI - MP2 法のマルチ GPU 実装の性能評価 :GPU による高速化と並列性能

3

TSUBAME2.5 全ノード規模を用いた
マルチ GPU 大規模 RI - MP2 計算

4

3 章で述べた性 能 評 価に続いて、TSUBAME 2.5 全ノード規 模を用本 研 究では、まず TSUBAME 2.5 におけるRI-MP2 計 算ジョブの実いて大きなナノ炭 素 分 子のRI-MP2エネルギーマルチGPU 計 算を

行 時間、並 列 性 能、GPUによる加速効果に関する性能評価を行っ

2 のRI-MP2実施した。本計算では、ナノグラフェン2 量体（C9 6 H 2 4 ）た。性 能 評 価のための計 算ジョブとして、ナノグラフェン C 96 H 24

計 算を 1,349ノードおよび 4,047 GPUを用いて実 行した。本 計 算のRI-MP2計算をTSUBAME 2.5の64から512ノードを用いて行った。

にはcc-pVTZ基底関数とcc-pVTZ補助基底関数を用いた（240原子、これらの計 算には cc-pVTZ 基 底 関 数 [12] と cc-pVTZ 補 助 基 底 関 数

6,342 原子軌道、600 占有軌道、5,832 仮想軌道、16,992 補助基底[13]を用いた（120 原子、3,216 原子軌道、300 占有軌道、2,916 仮想

関 数、対 称 性なし）。3 章の計 算と同 様にGPU 計 算では 1ノードあ軌道、8,496 補助基底関数、対称性なし）。MP2 計算の前に行われたり 3 MPIプロセスと 4 OpenMPスレッドを、CPU 計算では 1ノーるHF 計算は京コンピュータで実行し、得られた分子軌道データを

ドあたり 1 MPIプロセスと 12 OpenMPスレッドを用いた。

MP2 計 算の入 力データとして用いた。また、GPU 計 算では 1ノー

図 3 にGPUジョブとCPUジョブの実 行 時 間の詳 細プロファイル

ドあたり 3 MPIプロセスと 4 OpenMPスレッドを、CPU 計算では 1を示す。図 3 に示されるように、GPU 利 用により実 行 性 能の大 幅

ノードあたり 1 MPIプロセスと 12 OpenMPスレッドを用いた。表な向 上が達 成されている。特に、4 中 心 積 分を計 算するための行

1 に実行時間、並列性能、GPUによる加速率を示す。

列 - 行 列 積 演 算を GPU へオフロードするにより、この演 算が非 常に短時間で完了したため（ 図 3 の赤色の部分 ）、GPUジョブはCPUジョブと比較して 5.9 倍の高速化を達成した。実行時間はGPUジョブが 419 秒、CPUジョブが 2047 秒であった。また、実 測ピーク性
能はGPUジョブが 514.7 TFLOPsであったのに対して、CPUジョブ
は 87.5 TFLOPs であり、大 幅な実 測ピーク性 能 向 上を達 成した。
その一 方で、GPUジョブの実 測ピーク性 能 / 理 論ピーク性 能 比（10 % ）
はCPUジョブの結果（ 42 % ）よりも遥かに低い結果となった。
これはGPUジョブで最も時 間を要する処 理が 4 中 心 積 分のノード
内 集 団 通 信と 3 中 心 積 積 分のノード間 1 体 1 通 信となったのが原因である。更なる性 能 向 上のためにはこれらの処 理の性 能向上が必要となる。

表1

TSUBAME 2.5 におけるRI - MP2 計算の実行時間、並列加速率、GPUよる加速率（ナノグラフェンC 96 H 24 RI-MP2/cc-pVTZ ）いずれのノード数を用いた場 合でも、GPUジョブはCPUジョブと比較して 4.1- 6.6 倍の高速に計算が可能であることを確認した。
また、GPUジョブの並列性能はCPUジョブと比較して低いものの、
並列性能ノード数に対して良好にスケールすることを確認した。

図32 の RI- MP2 / cc-pVTZ
ナノグラフェン2 量体（ C 96 H 24 ）計 算の GPUジョブと CPUジョブの TSUBAME 2.5 における詳細プロファイル

15

マルチGPU 超 並 列クラスタシステムを用いた
大 規 模ナノ炭 素 分 子の電 子 状 態 計 算

ナノグラフェン 2 量体のπ - π 相互作用解析

5

おわりに

6

本 研 究で実 施した計 算 結 果より、ナノグラフェン2 量 体 間に働本研究では、最大でTSUBAME 2.5 1,349 ノード、4,047GPUを使用しくπ - π 相 互 作 用エネルギーの解 析を行った。モデル分 子にはてナノ炭素分子集合体のRI-MP2エネルギー計算を実施した。初め（ C 24 H 12 ）（ C 54 H 18 ）（ C 96 H 24 ）。モデルを作成す
2、2、2 を用いた（ 図 4 ）にGPUによる計 算の高 速 化の効 果と並 列 性 能の評 価をTSUBAME

る際にはABスタッキング構造（ 図 4（ a））を対象とし、C–C 結合距

2.5 512ノードまで用いて実 施した。GPU 実 装によりRI-MP2 計 算離（ 1.45 Å）と層間距離（ 3.35 Å）はバルクグラファイトの実験値

は大幅（4.1- 6.6 倍）な高速化が達成された。また、GPU 実装の並列を、C–H 結 合 距 離はベンゼンの実 験 値を用いた（ 1.1 Å）。相 互 作

性 能はノード数の増 加に対して良 好なスケーリングであった。次用エネルギーはMP2 法、Spin-component scaled MP2（SCS-MP2）に、本 研 究で実 施した最 大 のナノ炭 素 分 子 ( ナノグラフェン 2 量

法 [ 14] 、経験的分散力補正を施したダブルハイブリッドDFT 法であ体（ C 96 H 24 ）2 、240 原 子、6,342 原 子 軌 道 ）の RI - MP2エネルギ ーるB2PLYP-D3 法 [15] を用いて相互作用エネルギーの計算を行った。

を TSUBAME 2.5 全ノード規 模の 1,349ノードを用いて実 施した。

基 底 関 数と補 助 基 底 関 数は cc-pVTZ を用いた。基 底 関 数 重ね合GPUジョブでは 514.7 TFLOPsの実測ピーク性能を達成し、CPUジョ

わせ誤差はCounter poise 法 [16] で補正した。炭素原子 1 原子あた

ブのピーク性能（87.5 TFLOPs）からの大幅な性能向上を達成した。
さらに、2 層ナノグラフェン2 量 体に対してπ - πスタッキング

りのπ - π相互作用エネルギーの結果を表 2 にまとめる。

相 互 作 用エネルギーの解 析を行った。π - πスタッキング相 互 作
用エネルギーは分子のサイズが大きくなるにつれて、バルクの 2 層
グラファイトの相 互 作 用エネルギーに収 束するのを確 認した。また、SCS-MP2 法、B2PLYP-D3 法の計算結果はバルクグラファイトの解 離エネルギーの実 験 値を良く再 現した。一 方、MP2 法の結 果はこれらの方法よりも大きな相互作用エネルギーとなった。
本研究で使用したコードは 3 中心原子軌道積分計算ルーチンのGPGPU 実装も実施しておらず、またTSUBAME 2.5 のネットワーク
図42 層ナノグラフェン2 量体のモデル分子に対 応するように通 信の最 適 化も実 施していない。今 後の研究では、TSUBAME 2.5 のシステムに適するようにこれらの実 装を進める予定である。
謝辞

本 研 究を実 施するに際して、多くの有 益な研 究 議 論と助 言をいた
だいたNVIDIA Japanの平野幸彦氏に深く謝意を表す。本研究実施に伴うTSUBAME2.5 利用においては、TSUBAMEグランドチャレン
ジプログラムから全 面 的な支 援をいただいた。また、本 研 究で利表2ナノグラフェン2 量体の炭素原子当たりのπ -π 相互作用エネルギー（meV/atom）

用したプログラムの開 発およびテストでは「京」プロジェクト共通基盤研究の支援をいただいた。本研究を実施するにあたって事前に行った予 備 計 算では京コンピュータと自 然 科 学 研 究 機 構 岡

これらの結 果は、π - π 相 互 作 用エネルギーはモデルのサイズが大きくなるについてバルクの 2 層グラファイトの相 互 作 用エ

崎 共 同 利 用 施 設 計 算 科 学 研 究センターを利 用した。記して謝 意を表す。

ネルギーへ収 束していくことを示している。最 大のモデルである（ C 96 H 24 ) 2 のおける相 互 作 用エネルギーの計 算 結 果は、SCS-MP2

参考文献

法では 68.0 meV/atom 、B2PLYP-D3 法では 61.2 meV/atomであ[1]

C. Møller and M. Plesset: Note on an Approximation

り、バルクグラファイトの解 離エネルギーの実 験 値 52 ± 5 meV/Treatment for Many-Electron Systems, Phys. Rev., Vol. 46,

atom [17] をよく再現している。一方で、MP2 法の結果は、これらのpp. 618-622 (1934)

方 法の結 果や実 験 値よりはるかに過 大 評 価した結 果（ 85.7 meV/atom）となった。

[2]

T. Nakajima, M. Katouda, M. Kamiya, and Y. Nakatsuka:
NTChem: A high-performance software package forquantum molecular simulation, Int. J. Quant. Chem., Vol.

16

115, pp. 349-359 (2015)[3][4]

NTChem, http://labs.aics.riken.jp/nakajimat_top/ntchem_

pp. 553-566 (1970)[17] R. Zacharia, H. Ulbricht, and T. Hertel: Interlayer cohesive

e.htmlenergy of graphite from thermal desorption of polyaromatic

M. Feyereisen, G. Fitzgerald, and A. Komornicki: Use ofhydrocarbons, Phys. Rev. B, Vol. 69, pp. 155406 (2004)

approximate integrals in ab initio theory. An applicationin MP2 energy calculations, Chem. Phys. Lett., Vol. 208, pp.
359-363 (1993)[5]

D. E. Bernholdt and R. J. Harrison, Large-scale correlated
electronic structure calculations: the RI-MP2 method on
parallel computers, Chem. Phys. Lett., Vol. 250, pp. 477-484(1996)[6]

M. Katouda and T. Nakajima: MPI/OpenMP hybrid parallel
algorithm of resolution of identity second-order Møller–
Plesset perturbation calculation for massively parallel
multicore supercomputers, J. Chem. Theory Comput., Vol. 9,
pp. 5373-5380 (2013)[7]

M. Katouda, A. Naruse, and T. Nakajima: Massively parallel
algorithm and implementation of RI-MP2 energy calculation
for peta-scale many-core supercomputers, manuscript in
preparation.

[8]

BLAS, http://www.netlib.org/blas/[9]

M. Katouda and S. Nagase: Efficient parallel algorithm of
second-order Møller–Plesset perturbation theory withresolution-of-identity approximation (RI-MP2), Int. J. Quant.
Chem., Vol. 109, pp. 2121-2130 (2009)[10] NVIDIA CUDA ZONE, https://developer.nvidia.com/cudazone/[11] cuBLAS http://docs.nvidia.com/cuda/cublas/[12] T. H. Dunning Jr.: Gaussian basis sets for use in correlated
molecular calculations. I. The atoms boron through neonand hydrogen, J. Chem. Phys., Vol. 90, pp. 1007-1022 (1989)[13] F. Weigend, A. Köhn, and C. Hättig: Efficient use of thecorrelation consistent basis sets in resolution of the identity
MP2 calculations, J. Chem. Phys., Vol. 116, pp. 3175-3183(2002)[14] S . G r i m m e : I m p r o v e d s e c o n d - o r d e r M ø l l e r – P l e s s e t
perturbation theory by separate scaling of parallel-and
antiparallel-spin pair correlation energies, J. Chem. Phys.,Vol. 118, pp. 9095-9102 (2003)[15] L. Goerigk and S. Grimme: Efficient and Accurate DoubleHybrid-Meta-GGA Density Functionals—Evaluation with
the Extended GMTKN30 Database for General Main Group
Thermochemistry, Kinetics, and Noncovalent Interactions,J. Chem. Theory Comput., Vol. 7, pp. 291-309 (2011)[16] S. F. Boys and F. Bernardi: The calculation of small molecularinteractions by the differences of separate total energies.
Some procedures with reduced errors, Mol. Phys., Vol. 19,

17

● TSUBAME e-Science Journal vol.14
2016 年 2 月 29 日 東京工業大学 学術国際情報センター発行 ©
ISSN 2185-6028
デザイン・レイアウト：キックアンドパンチ
編集： TSUBAME e-Science Journal 編集室
青木尊之 渡邊寿雄 佐々木淳 板倉有希
住所： 〒 152-8550 東京都目黒区大岡山 2-12-1-E2-6
電話： 03-5734-2085 FAX：03-5734-3198
E-mail：tsubame_ j@sim.gsic.titech.ac.jp
URL： http://www.gsic.titech.ac.jp/

18

vol.

TSUBAME 共同利用サービ ス『みんなのスパコン』TSUBAME 共同利用サービスは、
ピーク性能 5.7PFlops 、18000CPUコア、4300GPU 搭載
世界トップクラスの東京工業大学 のス パ コン TSUBAME2.5 を東京工業大学以外の皆さまにもご利用いただくための制度です。

TSUBAME 共同利用サービスの利用区分とカテゴリ
TSUBAME 共同利用サービスには「学術利用」と「産業利用」の利用区分があります。
「学術利用」は、HPCI（革新的ハイパフォーマンスコンピューティングインフラ）、
JHPCN（学際大規模情報基盤共同利用・共同研究拠点）の採択により無償でご利用になれる制度と、
東京工業大学学術国際情報センターが実施する有償の制度がご利用になれます。
「産業利用」は、HPCI の採択により無償でご利用になれる制度と、
東京工業大学学術国際情報センターが実施する有償の制度があり、有償の制度には利用目的や利用成果を非公開とする「成果非公開」のカテゴリもあります。

TSUBAME 共同利用サービスの利用料金
TSUBAMEにおける計算資源は口数を課金単位としております。1 口は 3000ノード時間積で、
1 計算ノード（12CPUコア、3GPU 、58GBメモリ搭載）を 3000 時間、
あるいは 300 計算ノードを 10 時間というように、ご利用の用途に合わせ自由にご利用になれます。
TSUBAME 共同利用サービ スの 利 用区分・カテゴリ・利用料金を下 表 に示します。

お問い合わせ●東京工業大学 学術国際情報センター 共同利用推進室●e-mail kyoyo@gsic.titech.ac.jp

Tel. 03-5734-2085

Fax. 03-5734-3198

詳しくは http://www.gsic.titech.ac.jp/tsubame/をご覧ください。

14