64096 GPUを用いた 4096 3 規 模の一 様 等 方 性 乱 流の渦 法 解 析Turbulence Simulation Using 40963
Vortex Particles on 4096 GPUs

次 世 代シークエンサーから得られる
大 量メタゲノム情 報の解 析のための
超高速パイプライン
An Ultra-fast Computing Pipeline for Metagenome Analysis
with Next-Generation DNA Sequencers

大 規 模 並 列 GPU 計 算による
地震波伝播シミュレーション
GPU-Accelerated Large-Scale
Simulation of Seismic-Wave Propagation

34096 GPUを用いた4096 規模の一様等方性乱流の渦法解析
横田 理 央 *Lorena Barba***King Abdullah University of Science and Technology**** 慶應義塾大学 理工学部 機械工学科

成見 哲 *****ボストン大学

泰岡 顕治 ****

機械工学科*** 電気通信大学

情報理工学部

情報・通信工学科

高レイノルズ数の乱 流 解 析はハイパフォーマンスコンピューティングの困 難なアプリケーションの一つである。ここでは、40963 の渦粒子を用いた渦法による一様等方性乱流の解析を行い同等の格子点数を用いた擬スペクトル法との比較を行った。
二つの異なる手法を同 等のアプリーケーションについて比 較することで渦法の計算 精度の定 量的な検 証を行うことができた 。
また、4096 GPUを用いたときの高速多重 極 展 開 法を用いた渦法とFFT を用いた擬スペクトル法の並列 化 効 率の定 量 的な比較を行った 。

1

はじめにす。ただし、スケーラビリティはパフォーマンスとは直 接関 係のない
指 標であることに注 意されたい。寧ろ計算の遅いコードほど通 信を隠蔽することができるため、良 好なスケーラビリティを得やすい。このため、本 計算ではスケーラビリティだけでなく計算 時間も比 較し、
次 世代のハードウェアとアルゴリズムの協 調 設 計に役 立つような指

乱 流は我々の身の回りに広く存 在し、エネルギーや環 境に関する標を提 供するよう心がけた。

様々な工学上の問題に現れる。乱流 解 析は支配 方程 式の非 線 形性
により空間的にも時間的にも高い解 像 度を要するため大 規 模な計算
資 源が必 要となる。近 年、計算 機の発 達によりNavier-Stokes方程式の直接 数値 解 析が多くの場面で利用可能になってきている。
多くの工学的アプリケーションは複 雑な境界条 件や様々な物理 的

渦法

過 程との連 成が必 要となる一方で、このような複 雑な現 象から孤 立

2した素 過程を抽出し解 析することで乱流の普 遍的な性 質が明らかになることもある。一様 等方性乱流はこのような素 過程の一つであり、「渦法」と称される数値 解 析手法には様々なものがある [4] 。離 散

平均せん断 流や壁面の影 響を含まない高レイノルズ数乱流の基 本的化手法で大 別すれば、格子を一切用いないLagrange 型のものと、格な分析が可能となる。

子と粒子を併用するsemi-Lagrange 型の手法がある。支 配 方 程 式

一様 等 方 性 乱 流の直 接 数値 解 析にはこれまで擬スペクトル法がの定 式化で分けると渦度−流れ関 数によるものと渦度−速 度による用いられてきた 。最 大 規 模のものでは地 球シミュレータ上で 2002

ものが挙げられる。いずれの場 合も渦度 方程 式を解くことになるが、年に行われた 40963 格 子点、の計 算が挙げられる [2] 。石このとき粘性 拡散 項の粒子上での扱いにも様々な手法が提 案されて

原ら はこの大 規 模 直 接 数値 解 析 結果から、高 次の乱 流 統 計 量はいる。全ての「渦法」に共 通の特長としては、対 流 項がLagrange 的

高レイノルズ数で異なった挙 動を示すが、構 造関数のスケール指 数はに扱われるため、数値 拡 散や数値分 散の問題から解 放されるという

普 遍 的である可能 性が高いことを示した。このような結 論は最高性点が挙げられる。さらに、圧 力でなく渦度をベースとする定 式化を行能のスーパーコンピュータを利用して初めて得られるものであり、ハイうことで計算領域を渦度がある場所に限 定できる。これは渦輪や翼パフォーマンスコンピューテイングの有用性を象徴しているといえる。

端 渦のような渦 運 動が支 配 的な外 部流において計算点の大きな節

地 球シミュレータから京コンピュータにかけてLINPACKベンチ

約につながるため、演算とメモリ使 用量を大きく削減することができ

[1]

[3]

マークにおける演 算性 能が約 300 倍 向 上したにも関わらず、上 述のる。

4096 格 子点の計算 規 模の記 録は現 在でも破られていない 。これ

本 計算では完全メッシュフリーの離散化手法と渦 度−速 度によるは以下の二通りの解 釈が可能である。a）FFTはトーラスネットワーク

定 式化を採 用した [5] 。このためにまず渦度 場をGauss 分布を基底関を用いた分 散メモリ型の大 規 模 計算 機で性 能を発 揮できない。b）数にもつ渦 粒子の重ね合わせによって表す。渦 度−速 度による定 式

LINPACKベンチマークは実アプリケーションの性 能を予 測する指 標

化では速 度のPoisson方程 式

31

としては適当でない。前者の解 釈では、次世代計算 機で性能を発揮できるアルゴリズムの開発が必 要であるということになる。後者の見方では、LINPACKに変わる新たなベンチマークの採 用を促 進し、ハードウェアとソフトウェアの協調設 計に用いることが求められる。
本 研 究では一 様 等 方 性 乱 流の計算を行うための新たな手 法を提案し、TSUBAME2.0 上で良 好なスケーラビリティが得られることを示

01

[1]

81923 点を用いたスペクトル法のベンチマークの報告はされているが乱流の計算を長時間行ったものは著者らの知る限りまだない。

と渦度 方程 式

FMMの負荷 分 散の問 題を解 決する巧 妙な手 法として前のステッ
プの負荷のばらつきを元に現ステップの領 域 分割を更 新するものが挙げられる。このような手法は９０年代前半から共 有メモリ[8] 、分 散
メモリ[9] の負荷 分 散にそれぞれ用いられてきた。前のステップの情を交 互に解く。速 度のPoisson方程 式はGreen 関数を用いた積分方報を元に負荷分 散を微 調整する発 想自体はOrthogonal recursive

程式に変換するとBiot-Savart 式になる。これは全ての粒子間の相互bisection（ ORB）[7] 、Morton/Hilbert keyの分割法 [9] 、グラフ理論

作用を伴うN 体問題に帰着し、高速多重 極 展開法（ FMM）を用いてを用いた領 域 分割 法 [10] などの基 本 的な負荷 分 散の手法と併用する解くことができる。渦 度 方 程 式の各 項は個 別に解くが、更 新する変ことができる。本 計算に用いるFMMコードはORB 、Morton keyの数が各 項で異なるため部分 段 階 法には相当しない。対 流 項は粒子分割の両方から最 適なものを選ぶことができる仕 組みになっている。

の座標を更 新することで正 確に解くことができる。伸 張 項の計算は上 述の方法で粒子の領域分割が行われると、それぞれの領域で木速 度にBiot-Savart 式を代 入することで N 体 問 題に変 換できる。拡構 造が生 成される。次に、この局所的な木構 造の一 部を各プロセス

散 項は基底関数であるGauss 分布の標 準偏差σを増大させることでが相 互に送 信しながら各部分 領 域にある粒子が必 要な大 域 的な木考慮した。Lagrange 型の離散化手法が収 束するためには基底関数

構造（Local Essential Tree）を構築する必 要がある。このLETの構築が十 分重 複する必 要があるため、本 計算では数ステップに一回渦 粒に伴う通 信は大 域 的な通 信になるため、大 規 模な分 散メモリの計算

子の座標を等間隔に再配置し、基底関数のσを再 初期 化し、渦粒子ではここがスケーラビリティのボトルネックとなる。LETの通 信の最の強 度は元の渦 度 場を再 現するように調 節した。このとき渦強 度を適 化手 法としてORBの二分木の各 階 層を次 元に見 立てた超 立 方体

調節する計算は動径基底関数を用いた補間によって行った。ただし、に沿った通 信パターンを用いるものがある。これは treecode [11] にも動径基底関数には渦要素の基底関数に合わせてGauss分布を用いた。

FMM [12] にも適 用することができる。本 計算で用いたFMMコードでは上記の手法をプロセス数が 2 の冪乗でない場 合に拡 張したものを用いた。さらに、周期境界 FMMに上記の通信 方法が使えるようにも

3

高速多重極展開法 (FMM)

拡張した。
3.2 Dual Tree Traversal
粒子を領 域 分割し、各プロセスで局所 的に木 構 造を生 成し、LETを通 信するところまではtreecodeとFMMの手順は全く同様である。

Biot-Savart 式と渦度 方程式の伸張項の計算は全ての粒子同士の2

これら二手法 間の差異はLETを走 査する方法にある。ただし、以下

相 互作用の計算になるため直 接 解くとO(N ) の計算量になる。高 速の説 明では粒子間 相 互作用の作用される側をターゲット、作用を及

多重 極 展 開 法（ FMM）を用いることでこのときの計算量をO(N) に低

ぼす側をソースと定義する。また、木構 造の末端にあるセルを葉セル減することができる。

と呼ぶこととする。 木構 造の生 成 時には、葉セルあたりの粒子 数がBiot-Savart 式はPoisson方程式から導出できるためLaplace 型のできるだけ均 一になるようにする。Treecodeではターゲットの木 構

FMMカーネル関 数を用いることができる。Biot-Savart 式の右辺と造の葉セルに関するループをまわし、各葉セルに対してソースの木 構

渦 度 方程 式の伸 張項にはそれぞれ回転と勾配の空間微 分演 算子が造を走査する。一方、FMMでは通常木構造の走査は行わず、ターゲッ

存 在するためLaplaceカーネルの二 回 微 分までの項が必 要となる。

トセルに関するループをまわして各ターゲットセルについて作用する

FMMで用いる多重 極 展 開と局所展 開にはこれらの微 分が既に含ま

ソースセルのリストを逐 次 生 成する。粒子の分布が不均 一で木 構 造

れているため、この情 報をそのまま用いることができる。

に偏りがある場 合にはこのソースセルのリストの生 成が煩 雑になる。

3.1 負荷分散

このため、FMMでは通常 multipole acceptance criterion（ MAC）の概 念は用いず、影 響距 離の固定されたソースのリストを使 用する。

FMMを並列化する際に問題となるのはN 体問題に内在するデータ

Dual tree traversal[13] を用いることで上記のソースのリストの生 成の大 域 的な依 存性と粒子のダイナミックな性 質である。格 子法とはを簡略 化・最 適化することができる。これは、FMMでよく用いられる

異なり、N 体 問 題では領 域 分割 法を用いて粒子群を分割したとして

U,V,W,Xリスト[12] にMACの概念を付加したものと考えることもできる。

も問 題を局在 化することはできない。これは、それぞれの部分 領 域

図１にdual tree traversalの概 念図を示す。基 本的な手順はターゲッにある粒子が必 要とする情 報が結局全 領 域にまたがっているためで

トとソースの双 方の木 構 造を同 時に走 査するだけである。この際にある。FMMでは大 域 的な木 構 造を用いた領 域 分割 法が良く用いら生じるターゲットとソースのセルのペアに対してMACを適 用し、セルれるが、粒子の分布が非一様で木構 造に偏りがある場 合には領域を同士の計算を行うのに十 分な距 離にあるかを判定する。十 分でなけ

等 分割すると逆に負荷は偏ってしまう。さらに、粒子が毎ステップ移れば大きい方のセルを分割し、再びターゲットとソースのセルのペア動するため木 構 造を構 築し直す必 要があり、高 速に木 構 造を生 成、に対してMACによる判定を行う。これにより双 対な木構 造の走 査を

負荷分 散する機構が不可欠となる。

行い、ターゲットとソースの両方が葉セルに達した時点で距 離が未だ02

4096 GPUを用いた40963 規模の一様等方性乱流の渦法解析に十 分 遠くなければ直接計算を行う。Dual tree traversalはターゲッ

ことは容 易ではない。さらに、多重 極 展開や局所展開自体にも様々

トの親セルから子セルへソースの候 補を継 承することで、各階層で排な基底のものが存 在し、計算 機のアーキテクチャによってそのカーネ

他的なMACベースのソースリストを構築する最 適な手法を実現していルの性能は異なる。
これらのカーネルの中から最 適なものを選ぶ簡単な方法として、各るといえる。

カーネルの処 理 時 間を予め計 測し実 行 時にその情 報をもとに自動最 適化を行うものがある。Dual tree traversal 中にMACを満たすセルの各ペアについてP2P（ direct）、M2P（ treecode）、M2L（ FMM）のカーネルから最 適なものを選ぶことでいずれの単 独の手 法よりも常に高速な実 装となる。このことは本コードを用いた性 能 試 験で実証している [15] 。

スケーラビリティ

4

本 計 算は全て TSUBAME 2.0 上で行った。TSUBAME 2.0 の thinnode はIntel Xeon5670 を２ソケット、NVIDIA M2050 を３カード、54GBのRAM、120 GBのローカルSSDを有するノード1408台からなる。
ノード間はフルバイセクション・ノンブロッキングのfat-treeネットワー
クにより接 続されており、各ノードはdual linkのQDR Infinibandでリンクされている。
図 ２ にGPU あ た りN=40962 粒 子 を 割 り 当 て た と き の 1 か ら
4096GPUまでのweak scalingとそのときの計算時間の内訳を示す。
本 計算ではGPUあたりMPIプロセスを1つ割り当てたため、MPIプロ
セス数とGPU 数は等しい。計算時間はBiot-Savart 式と渦度 方程 式の伸張項の両方を含んでいる。 凡例の「 P2P evaluation」は近傍場のP2Pカーネルの計算時間、「FMM evaluation」はその他のFMMカー
ネルの合計 値、
「MPI communication」は全ての通信時間の合計 値、「GPU buffering」はGPU 周りのデータのバッファリングとGPUとの送受 信の合計、「Tree construction」は木構 造の生 成と領域分割、負荷分散の時間を表している。ただし、MPIの通信はローカルなP2Pカーネルの計算とオーバーラップしており棒グラフからはこの重 複 部分は引かれている。このため、棒グラフの合 計の高さは全 体の実 行 時 間と一 致している。最も大きい 4096GPUを用いた計算ではN=40963の粒子の計算を100 秒程 度で行い1.01 Pflop/sを達 成した。
図１ Dual tree traversalの概 念図

図２を見ると「GPU buffering」に時間がかかっていることが分かる。
ただし、これはFMMのカーネルをGPU上で効率的に実行するために必 要であり、実行時間の合計 値を低 減する効果があることが分かっ

3.3 CPU/GPU 上での自動最適化
TreecodeやFMMなどの階 層的なN 体 問 題の高 速 化 手 法では近体のスケーラビリティには全く影 響しない。FMMのスケーラビリティ

傍場のparticle-particle（ P2P）カーネルと遠 方場のmultipole-to-に効いてくるのは「 MPI communication」と「 Tree construction」

local（ M2L）やmultipole-to-particle（ M2P）カーネルが計算時間である。
「Tree construction」の中にもMPI 通信があり、これがスケーの大半を占める。また、この際の近傍場のP2Pと遠 方場のM2L/M2P

ラビリティを低下させている。3.1節で述べた超 立方体に沿った階層の計算時間のバランスが重要である。CPUとGPUが混 在する環 境で的なLETの通信を行う手法は、TSUBAME 2.0 のネットワーク上では単はこれらのカーネルの加速率にばらつきがあるため03ている。また、この部分は非常に良くスケーリングするため、FMM 全[14]

バランスをとる

純なMPI_Alltoallvに比べて低 速であった。このため、図２ではMPI_

Alltoallvを用いたときの通信時間を示している。
擬スペクトル法に関しても同様なweak scalingの試 験を行い、図３にこのときの並列 化 効 率をFMMと比 較したものを示す。ただし、本

一様等方性乱流の解析計算を行った時点ではCPU上のFFTWとGPU上のCUFFTの計算速 度

5に（ CPU-GPU 間の通 信時間を含めると）大きな差異が見受けられなかったため、擬スペクトル法の計算はノードあたり３プロセスを用いてここでは前述のFMMコードにより高速化された渦法とFFT WによCPU 上で行った。4096プロセスを用いたときの並列 化 効 率はFMMる擬スペクトル法を用いたでは 74%でスペクトル法では 14%であった。計算 時 間は両手 法とも果を示す。計算領 域は[-π,π ]3 の周期 的な立 方体で格 子点およびに１タイムステップあたり100 秒程 度であった。

粒子の数は 40963 であった。FMMに関しては 273 の周期 鏡 像を置くの一 様 等 方 性 乱 流の解 析 結ことで周期 境界 条 件を考慮した。また、FMMの級 数 展 開の次 数は
p=14 とした。初 期 条 件は指 定したエネルギースペクトルにランダムな位 相を持たせたものを波 数 空間の速 度 場として生 成した。擬スペ
クトル法はこの情 報をそのまま初 期 条 件として用いた。渦法では波数 空間の速 度 場をまず実 空間へと変 換し各セルの中心点（スタガー
ト格子における圧 力の計算点）における渦度を計算した。次に、渦粒子を同じくセルの中心点に配置し、動 径 基底関 数を用いた補 間法によって渦 度 場を再 現するように渦 粒子の強 度を計算した。渦 粒子の渦核半径はオーバラップ比が１となるようにσ =Δ xとした。
図４にt/ T=2 における渦法の計算 結果から得られた速 度勾配テン
ソルの第二不変 量の等 値 面を示す。ただし、Tは大 規 模 渦の回転時間である。微小な渦 構 造が多 数 見受けられるが、大 規 模はコヒーレ
ント渦構造は計算時間が不十分なため確認することができなかった。
計算時間はTSUBAME 2.0 のフルノードを占有できる時間の制約により僅かなものに限られた。ただし、渦法の計算精度やスケーラビリティ
2

図 2 GPUあたりN=4096 粒子を用いた
4096GPUまでのFMMの weak scalingの違いを定 量化するためにはこのような短時間の計算で十 分であるといえる。
図５にt/ T=2 における擬スペクトル法と渦法の速 度場から得られるエネルギースペクトルを示す。ただしここでは、擬スペクトル法の結果が正しいものとして渦法の計算 精度を検 証することを目的とする。
高波 数 成分において僅かな差異はみられるものの両手法間のスペク
トルは定 量 的に一 致していることが見てとれる。過 去の一様 等 方性
乱流の渦法による解 析 例 [5][16] と比べて高い精度が得られている要因
として以下のことが挙げられる。今 回の計算ではFMMの級 数 展 開の次 数が高く（ p=14）、周期 鏡像の数も多い（ 273 ）。また、再 初 期 化の頻 度が高く（5ステップに１回）、再 初期 化の際の動 径基底関数の補間の収 束判 定を厳しく設 定した（ |L|2 =1e-5）。もちろん、膨 大な数の渦粒子（ 40963 ）を用いていることも高い精度を実現できた要因の一つである。逆に、これらの条 件を全て満たさない限り高 精度な渦法計算は実現できないことが分かった。

図 3 FMMとスペクトル法の4096プロセスまでの並 列 化 効 率044096 GPUを用いた40963 規模の一様等方性乱流の渦法解析

おわりに

6

本研 究では一様 等方性 乱流を例にとって、二つの異なる計算手法
による解 析を同等の計算条 件の元で行い、計算 精度と計算速 度、ス
ケーラビリティに関する比 較を行った。FMMをベースとする渦法による解 析とFFTをベースとする擬スペクトル法の解 析を 40963 の粒子・格子点を用いて行いTSUBAME 2.0 上で最 大 4096 GPUまで用いた計算を行った。このときFMMの並列 化 効 率は 74%であったのに対してFFTの並列 化 効 率は 14%であった。両手 法とも計算 時 間は１ステッ
プあたり100 秒 程 度であり、40963 を用いた渦法計算は 1.01Pflop/sの持 続 演 算性 能に相当する。渦法と擬スペクトル法のエネルギース
ペクトルは定 量 的に一 致し、本手法の計算 精度を検 証することができた。ただし、TSUBAME 2.0 全 体を占有できる時間が限られていた
図 4 渦法の計算 結 果から得られた
速 度勾 配テンソルの第二不 変 量の等 値 面

ため、高次の乱流 統 計 量を得るのに十 分な計算時間を確保することができなかった。今後、同 様の計算資 源が安 価に利用できるようになれば、FMMを用いた乱流 解 析の有用性が明らかになると期待される。
謝辞

本研究の計算 機 環 境はTSUBAME2.0 の平成 23 年 秋 期グランドチャ
レンジ大 規 模 計算制 度によるものである。本研 究の一 部は科 学 技術振 興 機構 CREST「ポストペタスケール高性 能 計算に資するシステ
ムソフトウェア技 術の創出」から支 援を頂いた。ここに感 謝の意を
表する。
参考文献
[1] S. A. Orszag and G. S. J. Patterson: Numerical Simulation ofThree Dimensional Homogeneous Isotropic Turbulence. Phys.
Rev. Lett. Vol. 28, pp. 76-79 (1972)[2] M. Yokokawa, K. Itakura, A. Uno, T. Ishihara, and Y. Kaneda:
16.4-Tflops Direct Numerical Simulation of Turbulence by a
Fourier Spectral Method on the Earth Simulator, Proc. ACM/
図５ 擬スペクトル法と渦法の t/T=2 におけるエネルギースペクトル

IEEE Supercomput. Conf., Washington DC, pp. 50-66 (2002)[3] T. Ishihara, T., Gotoh, and Y. Kaneda: Study of High-Reynolds
Number Isotropic Turbulence by Direct Numerical Simulation,
Annu. Rev. Fluid Mech., Vol. 41, pp. 165-180 (2009)[4] G. H. Cottet and P. Koumoutsakos: Vortex Methods, Cambridge
University Press (2000)[5] R . Yokota, T. K . Sheel, and S. Obi: Calculation of Isotropic
Turbulence Using a Pure Lagrangian Vortex Method, J. Comput.
Phys., Vol. 226, pp. 1589-1606 (2007)[6] H. Cheng, L . Gre engard, and V. Rok hlin: A Fas t Adaptive
Multipole Algorithm in Three Dimensions, Vol. 155, pp. 468-498(1999)

05

[7] M. S. Warren and J. K. Salmon: Astrophysical N-body Simulation
Using Hierarchical Tree Data Struc tures, Proc. ACM/IEEE
Supercomput. Conf., pp. 570-576 (1992)[8] J. P. Singh, C. Holt, J. L. Hennessy, and A. Gupta: A Parallel
Adaptive Fast Multipole Method, Proc. ACM/IEEE Supercomput.
Conf., pp.54-65 (1993)[9] M. S. Warren and J. K . Salmon: A Parallel Hashed Oc t-tree
N-body Algorithm. Proc. ACM/IEEE Conf. Supercomput., pp. 1221 (1993)[10] S.-H. Teng. Provably Good Partitioning and Load Balancing
Algorithms for Parallel Adaptive N-body Simulation. SIAM J. Sci.
Comput., Vol. 19, pp. 635-656 (1998)[11] T. Hamada, K. Nitadori, K. Benkrid, Y. Ohno, G. Morimoto, T.
Masada, Y. Shibata, K. Oguri, and M. Taiji: A Novel Multiple-walkParallel Algorithm for the Barnes-Hut Treecode on GPUs –Towards Cost Effective, High Performance N-body Simulation,
Comput. Sci. – Res. Dev., Vol. 24, pp. 21-31 (2009)[12] I. Lashuk, A. Chandramowlishwaran, H. Langston, T.-A. Nguyen,
R. Sampath, A. Shringarpure, R. Vuduc, L. Ying, D. Zorin, and
G. Biros: A Massively Parallel Adaptive Fast Multipole Method
on Heterogeneous Architectures, Proc. Conf. High Perform.
Comput. Netw. Stor. Anal., (2009)[13] W. Dehnen: A Hierarchical O(N) Force Calculation Algorithm, J.
Comput. Phys., Vol. 179, pp. 27-42 (2002)[14] R. Yokota and L. A. Barba: Treecode and Fast Multipole Method
for N-body Simulation with CUDA, GPU Computing Gems,
Morgan Kaufmann (2011)[15] R. Yokota and L. A. Barba: Hierarchical N-body Simulations with
Auto-tuning for Heterogenous Systems. Comput. Sci. Eng., Vol.
14, pp. 30-39 (2012)[16] R. Yokota, T. Narumi, R. Sakamaki, S. Kameoka, S. Obi, and K.
Yasuoka: Fast Multipole Methods on a Cluster of GPUs for the
Meshless Simulation of Turbulence. Comput. Phys. Comm., Vol.
180, pp. 2066-2078 (2009)06

次世代シークエンサーから得られる
大量メタゲノム情報の解析のための
超高速パイプライン
石田 貴士
東京工業大学鈴木 脩司

大学院情報理工学研究科

秋山 泰
計算工学専攻近 年注目を集めているメタゲノム解 析は未 知の微 生物に関するゲノム情 報が得られるだけでなく、その環 境 中の共 生 系の理解や環 境 汚染の監 視等に有用であるが、解 析を進める上で計算量の大きな配 列 相 同 性 検 索 処 理がボトルネックの一つとなっている。この問 題に対処するため、我々はGPUを用いた高 速な配 列 相 同 性 検 索プログラムを開 発し、TSUBAME2 上に大 規
模なメタゲノム解 析向けの超高速パイプラインの構 築をおこなった 。

1

はじめに

メタゲノムマッピング

2

従 来のゲノム解 析は培 養された単 一の種のゲノム情 報を明らかに

次 世 代シークエンサーと呼ばれる現 在のシークエンサーは非 常にするものであったが、近 年シークエンサーの性 能 向 上に伴い、土壌、
海

高いスループットを持ち、１度の読み取り実 行で数 千 億 塩 基 以 上が解洋、
ヒト体 内 等の環 境 中に生 息する微 生 物のゲノムを分 離 培 養せず

読 可能となっている。
しかし、その出 力は 10 0 塩 基（bp）程 度の短い断にそのまま丸ごとシークエンスして解 析することが可能となってきた。こ

片 配 列であるため、そこから意 味のある情 報を得るにはその出 力に対の解 析 手 法はメタゲノム解 析と呼ばれ、未 知の微 生 物に関してのゲノ

してアセンブリやマッピングといった処 理を計 算 機 上で行う必 要があ

ム情 報が得られるだけでなく、その環 境 中の共 生 系の理 解や環 境 汚

る。

[1]

染の監 視 等に有用であり、大きな注目を集めている 。さらに近 年でマッピング処 理は既 知のゲノム配 列に対して断片 配 列を貼り付け、は次 世 代シークエンサーと呼ばれる新 型のシークエンサーの登 場にその一 致する位 置を同 定する処 理である。従 来 行われてきたようなより、膨 大な量のゲノム情 報が短 時 間のうちに入手可 能となっており、単一の生物 種に対する解 析ではリファレンスゲノムが存 在するため、
多その大 規 模な情 報を用いることでメタゲノム研 究が更に進 展すること

くの不 一 致やギャップを許 容する必 要がなく、通常の文 字 列 検 索に近が期 待されている。

い処 理によって解 析が可能であった。現 在では BWA [5] や Bow tie [6]としかし、このメタゲノムの解 析では、サンプルに含まれる多くの種のゲ

いった高 速なプログラムが開 発されており、数 台のワークステーション

ノム情 報がデータベースに登 録されていないため、遠 縁の種の配 列があれば次 世代シークエンサーの出 力に対しても十 分に対処が可能

データとの間で比 較が可 能となる高 感 度な検 索 手 法が必 要となる。

となっている。

配 列 相 同 性 検 索と呼ばれるこの検 索 処 理は多くの計算を必 要とする

その一方、メタゲノム解 析では環 境 中に含まれる微 生物のすべての処 理であり、
その結 果、
メタゲノム解 析を進める上でのボトルネックの一ゲノム配 列が既 知であることは稀であるため、マッピングにおいては近つとなってしまっている [2] 。

縁の種のゲノム情 報を参 照する必 要があり、曖 昧な一 致まで検 知 可

そこで、
本 研 究では次 世代シークエンサーによる大 量のメタゲノム情

能となるような高い感 度の検 索が必 要となる。そのような多くの不 一報を現 実 的な時 間 内に解 析することを目的とし、TSUBAME2 の膨 大致やギャップを許 容する検 索は一 般に配 列 相 同 性 検 索と呼ばれ、従な計 算 能 力を利 用 可 能とする大 規 模な全自動 解 析パイプラインを来のゲノムマッピングに比べて非 常に多くの計 算を必 要とする。さらに構 築した。また、そのパイプライン中で行われる配 列 相 同 性 検 索そのメタゲノム解 析ではより高い検 索 感 度を得るため、A，T，G，C の 4 文ものを高 速 化するため、従 来のメタゲノム研 究で標 準 的に利 用されて字で表 現されるDNA 配 列のままで検 索を行なわず、それらをコドン表[3]

きた配 列 相 同 性 検 索プログラムBLAST X に加え、BLAST Xと同 等に従い 20 種 類のアミノ酸からなるタンパク質 配 列に翻 訳してから解の検 索 感 度を持つオリジナルの GPUプログラムであるGHOSTM [4] を析を行う。DNA 配 列の比 較では塩 基の差は一 致か不 一 致の 2 状 態

開 発し、利 用 可 能とした。これによって次 世 代シークエンサーの一 度

でしか区 別しないことが多いのに対し、
タンパク質 配 列ではアミノ酸の読み取り実 行によって得られるゲノム情 報を数 時 間の内に処 理す

間で性 質の類 似 度の違いから、アミノ酸 置 換 毎に異なるスコアを用ることが可 能となり、今 後このパイプラインによって次 世代シークエンいるため、計算量は更に大きなものとなる。

サーによるメタゲノム解 析が促 進されることが期 待される。

このようなンパク質 配 列 間での曖 昧な一 致も含めた検 索を行うため、現 在では比 較 的 高 速かつ高 感 度で配 列 比 較が可 能な近 似 的
手 法 BLAST [3] が利 用されてきた [ 7 ] 。
しかし、次 世 代シークエンサーは
膨 大な量の DNA 断片 配 列を出 力するため、すべての DNA 断片 配 列

07

をアミノ酸 配 列に翻 訳してからマッピングを行うのは BLAST Xプログ

感 度を有しており、図２で示すように動 的 計 画 法によりアライメントをラムを用いても長い計 算 時 間が必 要となる。現 在、最 新の Illumina

行うSSEARCH [10] の結 果を正 解としたテストでは BLAST には劣るもの

社の Hiseq 20 0 0 DNA シークエンサーの出 力は 1 度の読み取りで合の、B L ATよりも高 い 検 索 感 度を示し、特 に実 際 に 利 用される Bit

計 60 0G 塩 基にも達し、そのデータを解 析するには約 25,0 0 0CPU 日

スコア50 以 上の領 域では BLASTとほぼ同 程 度の検 索 感 度を示してが必 要となってしまう。 そのため、現 在ではメタゲノムの解 析においており、BLASTの代 替として利 用可能な事が示されている。

は、このマッピング処 理が大きなボトルネックとなっており、その高 速 化
が必 要とされている。

また、表１は 約 75bp のリード10 0,0 0 0 本をクエリとして、KEGG のgenes.pepタンパク質 配 列データベース（約 2.5GB）に対して相 同 性
検 索を実 行した際の実 行 時 間であるが、GHOSTM の検 索 速 度は１
GPUを用 いた場 合、
1CPUコア上で 実 行され た BLAST X に 対して約GPU 配列相同性検索プログラム
GHOSTM3130 倍、4GPUを用いた際には約 4 0 0 倍の高 速 化を達 成している。これは検 索 感 度 落とすことで高 速 化を実 現しているBLATの高 速 化 率
約 4 0 倍に比べてもより高 速であり、GHOSTM は高 感 度と高 速 化を同時に実 現している。

我々はメタゲノム解 析において処 理のボトルネックとなっている配 列
相 同 性 検 索を高 速 化するため、GPUを利 用した高 速な配 列 相 同 性
解 析プログラムであるGHOSTM [4] を開 発した。GHOSTM のアルゴリ
ズムは BLASTと類 似したものであるが、GPU 上での実 行により適したものとなっており、メタゲノムの解 析に必 要 十 分な検 索 感 度を持つように設 計されている。GHOSTM は NVIDIA 社の CUDAを用いて実 装されており、
CUDA バージョン2.2 以 上が利 用 可 能なNVIDIA 社 製のGPU が 搭 載された計 算 機で 利 用 可 能である。一 部 の処 理を簡 略
化し、処 理の大 部 分をGPU 化することで大 幅な高 速 化を達 成しており、詳 細は後 述するが１GPUを用いた場 合、
1CPUコア上で実 行された
BLAST Xに対して約 130 倍の高 速 化を達 成している。
3.1 アルゴリズム

図 1 GHOSTM の処理の流れ

図１に示すように、GHOSTMでは BLASTと同 様にまず K 文 字の部
分文 字 列の一 致を探 索する事でアラインメントの候 補となる部 位を同 定し、その後に各アライメント候 補についてSmith-Watermanアルゴリズム [8] により局 所 的アラインメントを行う事で詳 細なアラインメン
トとそのスコアを計 算している。そして、最 終 的に全てのアラインメントのスコアがソートされ、そのスコア上位のヒットが検 索 結 果として出 力される。
GHOSTMではこれらの処 理のうち、計算の大 部分を占めるアライン
メント候 補 探 索と局所 的アラインメントの双 方をGPU上で処 理している。GPU は多くのコアを搭 載しており、処 理を適 切に並 列 化することで大きな高 速 化が可 能であるが、GHOSTMでは大 量の断片 配 列を処理する必 要があるため、
まず、
アラインメント候 補の探 索の際には各 断
片 配 列の処 理を各 GPUスレッドで行い、局 所 的アラインメントの際には各アラインメント候 補に対する処 理を GPUスレッドに割り当てて処理することで並 列 化を行っている。
3.2 検索速度と検索感度

図 2 GHOSTM の検索感度

メタゲノム解 析では高い検 索 感 度が必 要となるため、BLAT [9] のようなBLASTより高 速であるがに検 索 感 度の低いアルゴリズムを利 用する事は困 難であった。一方 GHOSTM はメタゲノム解 析に十 分な検 索08次世代シークエンサーから得られる
大量メタゲノム情報の解析のための
超高速パイプライン

million の DNA リードの処 理を実 現した。GHOSTMと2,520 GPUを用いた場 合では 1,260GPUを用いた場 合に対して速 度の向 上が 2 割
程 度に留まっているが、これは計 算資 源に対して今 回 用いたデータのサイズが小さくなりすぎてしまったため、負荷 分 散 等に失 敗したことが
原 因であり、最 新のシークエンサーの出 力や、複 数 回のシークエンシ
ングの結 果を同 時に処 理すればGPU 数に対する速 度 向 上はほぼ線形となると考えられる。

表 1 実行速度の比較

大規模全自動解析パイプライン

4

我々はメタゲノム 解 析 に お けるボトルネックである配 列 相 同 性
解 析を高 速 化 するため、GPUを 利 用した 高 速 なプログラムである
GHOSTMを作 成した。しかし、一 台の計 算 機では、
この GHOSTMを用いても次 世代シークエンサーの出力を解 析するには不十 分である。
そのため、我々はTSUBAME2 の大 量の計 算ノードを解 析に利 用することで、次 世 代シークエンサーによるメタゲノム解 析を現 実 的な時 間

図 3 BLASTX を用いた際のCPU コア数に対する処理速度の向上

で行うことを目指し、
TSUBAME2 上に全自動のメタゲノム解 析パイプ
ラインを構築した。
多くの計算ノードを利用する大 規模な計算ではデー
タベースのコピー、計 算 結 果の書き込みといったファイルの入 出 力に関する処 理が問 題となるが、このパイプラインではノード間でコピーを
2 分 木 状に行う工夫を行う事で入 出 力に依 存するボトルネックを解 消している。
4.1 大規模実行
我々は本パイプラインについてTSUBAMEグランドチャレンジ（超 大規 模アプリケーション）制 度を利 用し、次 世代シークエンサーによって得られた大 規 模なデータに対する解 析 速 度の検 証を行った。解 析に用いたデータは汚 染 土 壌に関するメタゲノムデータである。メタゲノム
原データは各 75bp の DNAリード224 million 本であるが、実 際に配列 相 同 性 検 索が行われたのは、
ここから低 品 質なデータを除 去 後した 71million 本の DNAリードであり、これらがクエリとして4.2GB のタン
パク質 配 列 情 報が含まれたNCBI nrデータベースに対する検 索が行われた。パイプラインの相 同 性 検 索エンジンとしては CPU 上で動 作
するBLAST Xと我々の開 発した GPU 上で動 作するGHOSTM のそれぞれについて実 効 性 能の比 較を行った。
パイプラインは計算コア数に対してほぼ線 形の速 度 向 上を示し、
図３のようにBLAST Xを相 同 性 検 索に用いた場 合、
TSUBAME2 の CPU
16,0 08コア（1,334ノード）用いた際に 1 時 間あたり約 24 million のDNAリードの処 理を実 現した。また、
GHOSTMを用いた場 合では図４のように、2,520 GPU（8 4 0ノード）を用いることで 1 時 間あたり約 60

09

図 4 GHOSTM を用いた際のGPU 枚数に対する処理速度の向上

おわりに

5

[8] T. F. Smith and M. S. Waterman, “Identification of common
molecular subsequences.,” Journal of molecular biology,
vol. 147, no. 1, pp. 195-7, Mar. 1981.
[9] W. J. Kent, “BLAT--the BLAST-like alignment tool.,” Genome
research, vol. 12, no. 4, pp. 656-64, Apr. 2002.

我々は次 世代シークエンサーから得られるメタゲノム情 報の解 析に[10] W. R . Pear s o n, “ S earchin g p rotein s e qu en ce lib r ar i e s:

おいて現 在ボトルネックとなっている相 同 性 検 索を高 速 化することを目comparison of the sensitivity and selectivity of the Smith-

的とし、GPUによる効率的な相同性検索プログラムGHOSTMを開発し、Waterman and FASTA algorithms.,” Genomics, vol. 11, no. 3,

またTSUBAME2 の大 量の計算 機 資 源を利 用可能とする大 規 模な解

pp. 635-50, Nov. 1991.

析パイプラインの構 築を行った。解 析パイプラインは使 用したCPUコ
ア数、GPU 数に対してほぼ線 形の速 度 向 上を示し、TSUBAME2 のほぼ全 体を利 用した場 合、
1 時 間あたり約 60 millionのDNA リードの処理が可能であった。これは次 世代シークエンサーの 1 度の読み取り実行から得られる出力に対する解 析が数 時 間 以内に可能となる事を示しており、今後我々のパイプラインによって次 世代シークエンサーによる
メタゲノム解 析が促 進されると考えている。
謝

辞

本 研 究 の 一 部 は HPCI 戦 略プログラム「予 測 する生 命 科 学・医 療および創 薬 基 盤」、国 立がんセンター研 究 所がん研 究 開 発 費、及び
NVIDIA 社 CUDA COEプログラムの支 援を受けて行われたものであ
る。
参考文献
[1]

M . A r u m u g a m e t a l . , “ E n te r o t y p e s o f t h e h u m a n g u t
microbiome.,” Nature, vol. 473, no. 73 4 6, pp. 174 - 8 0, May
2011.

[2] J. C . Wo o ley, A . G o d zik , and I. Frie db erg, “A primer on
metagenomics.,” PLoS computational biology, vol. 6, no. 2, p.
e1000667, Jan. 2010.
[3] S. F. Altschul, W. Gish, W. Miller, E. W. Myers, and D. J. Lipman,“Basic local alignment search tool.,” Journal of molecularbiology, vol. 215, no. 3, pp. 403-10, Oct. 1990.
[4] S . S u z u k i , T. I s h i d a , K . K u r o k a w a , a n d Y. A k i y a m a ,“GHOSTM: A GPU-Accelerated Homology Search Tool for
Metagenomics,” PLoS ONE, vol. 7, no. 5, p. e36060, May 2012.
[5] H. Li and R. Durbin, “Fast and accurate short read alignment
with Burrows-Wheeler transform.,” Bioinformatics (Oxford,
England), vol. 25, no. 14, pp. 1754-60, Jul. 2009.
[6] B. Langmead, C. Trapnell, M. Pop, and S. L. Salzberg, “Ultrafast
and memory-efficient alignment of short DNA sequences tothe human genome.,” Genome biology, vol. 10, no. 3, p. R25,
Jan. 2009.
[7] P. J. Turnbaugh, R . E. Ley, M. a Mahowald, V. Magrini, E.
R . Mardis, and J. I. Gordon, “An obesit y-associate d gut
microbiome with increased capacity for energy harvest.,”Nature, vol. 444, no. 7122, pp. 1027-31, Dec. 2006.

10

大 規模 並列 GPU 計算による
地震波伝 播シミュレーション
岡元太郎 *

竹中博士 *** 東京工業大学理工学研究科

中村武史 ***** 九州大学理学研究院

青木尊之 ******* 海洋研究開発機構**** 東京工業大学学術国際情報センター

2011 年の東北 地 方太平洋沖 地 震は、強い地 震 動と巨 大な津 波によって東日本 地 域に計り知れないほどの被害をもたらした 。
この地 震について、巨 大 地 震の震 源がどのようなものであったか、また巨 大 津 波や地 震 動の励 起はどのようになされたのかを探ることが、地 球 科 学においてきわめて重 要な課 題となっている。そのためには大 規 模な地 震 波 計算によって得られる理 論
波 形を用いた定 量 的な解 析が必 要となる。我々はTSUBAMEの GPUを用いた高速・大 規 模 計算によってこの問 題に取り組んでいる。本 稿では、我々が採 用したGPU 計算手法や計算性能について報告する。また、東北 地方太平洋沖地震に関するシミュ
レーション例を紹介する。

はじめに

1

差 分 法にもとづいたものである。この手 法は地 震 波 伝 播シミュレー
ションにおける時 間 領 域 解 法の標 準 的な手 法の一つとして広く用いられているものである。プログラムの主 要 部 分の開 発には東 京 工 業
大 学 学 術 国 際 情 報センターのTSUBAME-1.2（2010年10 月まで）
、および TSUBAME-2.0（ 2010 年 11月以 降）を利 用させていただいた。こ2011 年 3 月11日に発 生した東 北 地 方 太平洋 沖 地 震（マグニチュー

こではTSUBAME-2.0 による結 果について報 告する。

ド９、図１）は、強い地 震 動と巨 大な津 波によって東日本 地 域に計り知なお本 稿は地 震の研 究をテーマとしているが、 地 震 波 計 算は地れないほどの被 害をもたらした。
しかし残 念ながら、この地 域ではこれ

球・惑 星 内 部 構 造を推 定するうえでも非 常に重 要な手 法であること

ほどの巨 大 地 震の発 生は想 定されていなかった。そのため、巨 大 地を付 言しておきたい。特に資 源 探 査の分 野では地 震 波 伝 播を高 速

震の発 生 条 件や強 震 動・巨 大 津 波 励 起の仕 組みを探ることと、今に計 算する手 法や高 性 能ハードウェアへの要 請が非 常に高く、GPU後の地 震防災に向けた研 究を進めることが地 球 科 学においてきわめ

計 算を応 用した地 震 波 計 算や内 部 構 造 推 定の研 究が世界 的に進て重 要な課 題となっている。

められている。

我々は地 震 波を解 析することによって、この地 震がどのようなものであったのかを推 定し、地 震 時に発 生した強 震 動や津 波 生 成との関 連を探る研 究を進めている。 その研 究では、地 球 内 部 不 均 質 性や地形・海 水 層などの効 果を考 慮した大 規 模な地 震 波の理 論 計 算を繰り返し実 行しなければならない。そのため、大 規 模で高 速な計 算 資
源が必 要となっている。そこで我々は GPU の演 算性 能に着目し、大 規
模 並 列 GPU 計算による地 震 波 伝 播 計算プログラムの開 発を進めている。
GPU（Graphics Processing Unit）は、その名 称の通り、画 像 処 理を高 速に実 行するためのデバイスである。GPU の特 徴は、非 常に多くの演 算 用コアを内 蔵していることと、極めて高い数 値 演 算性 能を持つ
ことである
（図２）
。最 近の GPU は一 千 個 以 上の演 算コアを内 蔵し、単精 度 演 算では 2TFlops（テラ・フロップス）
を越える理 論ピーク演 算性
能を持っている。また、絶 対 性 能と同 時に電 力あたりの性 能も高い。
さらに一 般 消 費 者 向けに大 量 生 産される製 品であることから、高い
演 算性 能にも関わらず安 価にかつ容 易に入手できるのも GPU の大きな特 徴である。
GPU にはもう一つ、
メモリ転 送 速 度が非 常に高 速であるという重 要な特 徴がある（図２）。そのため、多 量のメモリ読み出し・書き込みが発生するメモリ集 中 型（memor y intensive）の計算 問 題でも高い性 能を発 揮する。地 震 波 伝 播のシミュレーションもメモリ集 中 型の問 題であり、GPU を応 用することには大きな利 点がある。
本 稿では、我々が開発してきたGPUプログラムについて紹介する[1-3] 。
我々のプログラムは格 子 計算 型の手 法、
すなわち食い違い格 子 型の11

図1
東北地方太平洋沖地震（2011 年 3 月 11 日）の余
震 分 布（2011 年 3 月 9 日 か ら 2011 年 4 月 11 日、
MJMA ≧ 4.0）。マグニチュード 7 以上の地震の震央を星印で示した。震源データは気象庁による。

GPU 計算3

メモリ最 適 化
地 震 波の差 分 法 計算のようなメモリ集 中 型の問 題では、
GPU 内 部
図 2 本研究で利用した GPU（NVIDIA M5020）の構成の概念図。

のメモリ（グローバルメモリ）に差 分 法 領 域の変 数 全てを保 存し、その変 数に対して演 算を行うことになる（図４a）。このグローバルメモリと演算ユニット間のメモリ帯域 幅は、通常の CPU の帯域 幅に比べると非 常に大きい。
しかしそれでもグローバルメモリから演 算ユニットへのデー

計算スキームについて

2

本 稿 で 紹 介 するプログラムで は、時 間 領 域 の 差 分 法（FDTD:
finite-dif ference time domain）を用いている [4] 。この差 分 法では媒 質の粒子 速 度（v i ; i=x, y, z）と応 力（τij ; i,j=x, y, z）とを変 数として、
図３に示す食い違い格 子を用いて計 算 領 域を離 散 化する。このような食い違い格 子は地 震 波 計算のほか、
電 磁 界シミュレーションでも広く利 用されている [5] 。なお、本 稿で紹 介するプログラムでは、数値 分 散性を向 上させるために空 間 差 分 精 度は 4 次 精 度とした。時 間 差 分 精
度は 2 次 精 度である。 完 全 弾 性 体の場 合には１単 位セルあたりのデータは粒 子 速 度・応 力が計 9 個、物 性パラメータが 3 個（τij と同じ

タ転 送には 4 0 0 - 60 0サイクルの遅 延が発 生する。そ のため、演 算ユ
ニット内 部にある高 速な共 有メモリとレジスタを「キャッシュメモリ」として利 用することが性 能 向 上のためには重 要である。
先に述べたように地 震 波 計 算では単 位セルあたりの変 数が多く、その一方で共 有メモリのサイズは小さいため（M2050 の場 合 16 kB
または 4 8 kB）、共 有メモリに 3 次 元のブロックを確 保することは困 難である。そこで我々のプログラムでは共 有メモリとレジスタとを併 用する（図４a）。すなわち横 方 向への差 分 操 作が必 要になる 2 次 元 平 面
上のセルのデータは共 有メモリに配 置し、
同 時にこの 2 次 元平 面 上の各セルに一つのスレッドを割りあてる。こうすると図の縦 方 向に位 置す
るセルのデータについては縦 方 向への差 分のみが発 生しスレッド間で共 有する必 要がないので、これらのデータは各スレッドのレジスタに配 置すればよい（図４b）。

格 子点に置く）、
合わせて12 変 数となる。非 弾 性 減 衰を含めた場 合には変 数の個 数がこの数 倍に増 加する。
（本 稿で紹 介する結 果は全て
弾 性 体 版のものである。
）
このスキームでは食い違い格 子を用いるため、粒子 速 度 変 数と応力変数の時刻が時間刻みの半刻みだけずれることになる。計算上は、
粒子 速 度と応 力を片 方ずつ更 新していくことになる。

図 3 食い違い格子と変数の配置。
図示した格子は単位セルに相当する。

図4
共有メモリとレジスタを併用する方法の概念図。(a）副領域の変数をすべて GPU のグローバルメモリに収める。さらに、その中の小さな 2 次元ブロック（赤色）を共有メモリにコピーする。青矢印は計算が進行する
方向を示す。（b）差分法計算に関係するデータ配置の概念図。赤色の格子点のデータは共有メモリ、灰色の格子点のデータはレジスタにそれぞれコピーする。

12大 規 模 並列 GPU 計算による
地 震 波伝 播シミュレーション

そして、縦 方 向へ計算を進めていく際に共 有メモリまたはレジスタに
一 度 読み込んだデータを、一 種のパイプラインのようにレジスタと共有メモリとの間で交 換することによって、グローバルメモリへのアクセスを低 減することができる。このようなレジスタと共 有メモリの利 用 法は
地 震 波シミュレーション [1,6 -9] や気 象シミュレーション [10] で利 用されて効果をもたらしている。
この最 適 化のほかに、物 性パラメータに関して参 照テーブル方 式を採 用することや、単 位セル中 心 以 外の格 子点での物 性パラメータを
セル中 心の物 性 値を用いて毎 回 計算すること[11] などの方 法によって、グローバルメモリへのアクセスを減らすようにしている。
このようにして作成したプログラムの性能を図５に示す。この性能は、次に述べるマルチGPU 計算の機 能を組み込んだプログラムを用いて、

図 6 3 次元領域分割の概念図。青い部分が袖領域を示す。

単 一 GPUで測 定したものである。比 較のために、TSUBAME-2.0 のホストノードCPU を用いて袖 領 域 処 理やMPI 機 能を含まない差 分 法
プログラムの性 能を測 定した。参 考 的な比 較となるが、これらのプロ

3 次 元 領 域 分 割では隣 接する副 領 域との間で「袖 領 域」
（ ghostグラムに関しては GPU 1 基はホストノードのおよそ 3 倍の性 能を持つ。

zone）を交 換する必 要がある。GPU 計算では、
この袖 領 域のメモリ並びに関して問 題が発 生する。これは、現 状では GPU 間で直 接データ
を交 換することはできないことと関 係している。図２に示すように、
デー
タを交 換するためにはまずホスト計 算 機のメモリにGPUからデータを転 送し、そのデータを別のホストに送るという段 階を経なければならない。
袖 領 域のメモリを割り当てる方 法としては、図７（a）のように内 部 領
域の配 列を延 長する方 法がしばしば採 用される。
しかし、袖 領 域 内 部
ではメモリ並びは不 連 続となるため、袖 領 域のデータをホスト計 算 機のメモリに転 送する際には cudaMemcpyなどのデータ転 送 関 数を
細 切れのデータごとに繰り返し呼び出すことになる。そのため袖 領 域
全 体を転 送するには非 常に長い時 間がかかることになり、この方 法はGPU 計算では実 用にならない。そこで、本 研 究では図７
（b）
のように内
部 領 域とは独 立の連 続なメモリを袖 領 域に割り当てる。こうすると、一度のデータ転 送 関 数 呼び出しで袖 領 域 全 体を GPU からホストに転送することが可能となり、データ転 送 時 間を（a）の場 合よりも大 幅に短図5GPU と CPU での地震波計算プログラムの性能。GPU
の性能は領域サイズ 320 × 320 × 320 の場合について
単一 GPU を用いて測定した。CPU では、袖領域処理を含まない差分法プログラムを Fortran と Open MP を用いて作成し、PGI Fortran コンパイラ（-fastsse オプ
ションを利用）でコンパイルした。領域サイズは 320
× 320 × 320（1 コアの場合）、および 320 × 320 ×
3840（12 コアの場合）とした。

くすることができる。ただしGPU 内 部では、袖 領 域メモリのデータを共有メモリのブロックに転 送するための複 雑なマッピングが追 加される。

袖 領 域 転 送の効 率 化
GPU のグローバルメモリのサイズは大きくない（M2050 の場 合では 3 GB）。そのため大 規 模 計算では、計算領 域を多 数の副 領 域に分割して複 数の GPU に副 領 域を割り当てることが不 可欠となる。本 研 究
では大 規模計算で標 準的に利用される3 次 元領域分割を用いる（図６）。副 領 域 間の通 信には MPI を用いる。

13

図7
袖領域（ghost zone）のメモリ並びの概念図 [1]。（a）内部領域の配列の延長として割り当てた場合。袖領域
内部では数多くの小さなデータが不連続に並んだメモ
リ並びになる。（b）内部領域とは独立の連続なメモリ
領域を割り当てた場合。

計算と通 信のオーバーラップ
食い違い格 子では粒子 速 度 成 分と応力成 分を交 互に時 間 積 分し東北地方太平洋沖地震

4ていく。そのため、それぞれの積 分 演 算を行うGPUカーネルも別々のものになる。さらに我々のプログラムでは、計算に要する全 体の時 間を

ここで 紹 介する計 算 例では、全 世 界 観 測 網の広 帯 域 地 震 波 形

短縮するために、内部領域についての積分と袖領域 通信とをオーバー

データを用いて推 定した震 源モデル [12] を用いる。遠 方の観 測 点のラップさせている。そのため、側 面ブロックのみを処 理するカーネルとデータを用いているため分 解 能に制 限があり、やや長 波 長 側の震 源

内 部ブロックのみを処 理するカーネルの 2 つのカーネルを作 成して利モデルになっている。この震 源モデルは断層面 上におよそ20kmおき用している。計算と通 信のオーバーラップのためには、
GPUとホスト間のに配 置した格 子点（ 23 個× 11 個）に割り当てた点 震 源から構 成されメモリ転送においてCUDA の非同期通信関数（cudaMemcpyAsync）

る。東日本 地 域の 3 次 元 不 均 質 速 度 構 造モデルは、複 数の既 存 構を用いる必 要がある。また、各 MPIプロセス間のデータ転 送には MPI造モデル [13-16] をコンパイルして作 成した。差 分 法 計算のパラメータはの non-blocking 関 数（MPI_Isend 、
MPI_Irec v）
を用いた。

次の通りである：格 子 間 隔 0.15 km 、時 間 間 隔 0.0 05 s 、格 子サイズ

こ の よう に し て 作 成 し た マル チ GPUプログラム に つ い て、

6 4 0 0 × 320 0 × 160 0 、時 間ステップ数 4 4 0 0 0 、周 波 数 帯 域の上 限

TSUBAME-2.0 で全 領 域サイズを変えて計 算を実 行したときの、GPU

0.61 Hz 。この計 算では TSUBAME-2.0 の 10 0 0 GPUを利 用した。こ数に対する性 能を図８に示す。この例では副 領 域サイズを固 定していれは TSUBAME-2.0 の全 GPU の約 24%にあたる。 なお、計 算 時 間るので、グラフは弱スケーリング性 能を示すことになる。図からわかるよは 5768 s、演 算性 能は 33.2 TFlops である。図８に示した例よりも性

うに、80 0 GPU の場 合まで理 想 的なスケーリングに近い（GPU 数に能が低 下したが、これは計算 結 果のファイル出 力に伴うオーバーヘッド

ほぼ比例する）実 効 性 能が得られた。絶 対値としても、80 0 GPU の場が大きいためと考えられる。
この点も今後の改良が必 要な部分である。

合に約 50 TFlops 、
120 0 GPU の場 合に約 61 TFlops という、非 常に高い実 効 性 能を達 成できた
（いずれも単 精 度 性 能）[3] 。

このようにして計 算した地 震 波 動 場を可 視 化したものを図９と図１０に示す [2] 。断 層 各 部からの波 動 場が、地 形や構 造による散 乱もなお、
120 0 GPU のときには理 想 的な場 合よりも少し性 能が低 下し加わって非 常に複 雑に干 渉しながら伝 播していく様 子がわかる。またているように見える。80 0 GPU の場 合の性 能は理 想 的な数 値の約130 秒 後から150 秒 後にかけて、上 向き速 度と下向き速 度の領 域が85%、120 0 GPU の場 合には約 69%である。この計 測では 80 0GPU

顕 著なペアになって福 島 県 海 岸 線 付 近に現れ、震 央 付 近から到 来までは 2 GPU/node 、
120 0 GPU のときは 3 GPU/node の条 件で計する波 動のパタンを掻き消すような様 子になっていることもわかる（図算した。そのため、後 者の場 合にはノード間 通 信が増 加してスケーラビ

１０、黄 色の円 形 枠）
。このペアの領 域をもたらした波 源 域は福 島 県

リティがやや低 下したと考えられる。この点は今 後 改 良を進める必 要

海 岸 線 付 近の深さ約 60km にあるやや大きなすべり領 域と考えられがある。

る。このような領 域は陸 上 観 測 点の強 震 動 記 録 からも推 定されており、福 島 県 海 岸 線 付 近を始め数か所の SMGA（strong-motiongeneration-area）
が推 定されている（例えば [17 ] ）。このような領 域を含めた震 源モデルの全 体像について、研 究を続ける必 要がある。

図8
TSUBAME-2.0 での実効性能。副領域サイズを
320 × 320 × 320 に固定して、全領域サイズを変えて（つまり分割数 =GPU 数を変えて）弱スケーリング性能を測定したもの。 いくつかの点での単位セル数も矢印と数値のペアで示す。

14

大 規 模 並列 GPU 計算による
地 震 波伝 播シミュレーション

図 10 破壊開始点から 130 秒後、150 秒後、
170 秒後、190 秒後のスナップショット。

まとめ5

地 震 波 伝 播シミュレーション、特に東 北 地 方 太平洋 沖 地 震のような巨 大 地 震のシミュレーションを行うために、
マルチGPUを用いた並列 差 分 法 計 算 手 法を開 発した。今 回、TSUBAME-2.0 の 10 0 0 GPUを使うことによって、330 億 個の単 位セルからなる領 域について上 限
周波 数が 0.61 Hz であるようなシミュレーションを実 用 的な時 間 内で
図9
シミュレーションによる地震波伝播のスナップショット。地表および海底の固体側表面での鉛直地動速度を可視化したもの。赤と青はそれぞれ上方向、下方向への動きを示す。上から下にかけて、破壊開始時刻から 10 秒後、30 秒後、50 秒後、70 秒後、90 秒後、110秒後のスナップショットを示す。海水中の音波伝播は計算で考慮さ
れている。

15

実 行することができた。計 算された地 震 波 動 場は、観 測データの（やや長 周期 側での）特 徴を定 性 的に再 現している。このような計算が可能になったことから、
今 後は、日本列島の各 地で観 測された地 震 波 波
形データを用いた震 源モデル解 析を行うための『グリーン関 数 波 形』
をTSUBAME-2.0 によって生 成できる。そのようにして計算したグリーン
関 数 波 形を用いて震 源モデルの改 良と推 定を進める必 要がある。

謝辞e-Science Journal, vol. 2, pp. 9-13 (2010).

構 造モデルデータを提 供してくださった各 機 関の方々に感 謝申し上[11] Takenaka, H., Nakamura, T., Okamoto, T. and Kaneda, Y. “Aげます。また、TSUBAME-2.0 における大 規 模 計 算（グランドチャレン

unified approach implementing land and ocean-bottom

ジ）の機 会をご提 供くださった学 術 国 際 情 報センターの方々に感 謝申

to p o gr ap hie s in th e s t a g g ere d - gr id f ini te - di f feren ce

し上げます。この研 究の一 部に科 学 研 究 費 補 助 金（23310122）
を使

method for seismic wave modeling”, Proc. 9 th SEGJ Int.

用しました。

Symp., CD-ROM Paper No.37 (2009).
[12] O k a m o to, T. , Ta ke n a k a , H . , H a r a , T. , N a k a m u r a , T. a n d

参考文献

Aoki, T. “Rupture Process And Waveform Modeling of The[1]

O k a m o to , T. , Ta k e n a k a , H . , N a k a m u r a , T. a n d A o k i , T.

2 0 11 Tohoku- Oki, Magnitude - 9 Ear thquake”, American“Acce l e r at i n g l a rg e - s c a l e sim u l at i o n o f s e ism i c w ave

Geophysical Union, Fall Meeting, U51B-0038, San Francisco,propagation by multi-GPUs and three-dimensional domain

USA (2011).

decomposition”, Earth, Planets and Space, vol. 6 2 , no. 12 ,pp. 939-942 (2010).
[2] O k a m o to, T. , Ta ke n a k a , H . , N a k a m u r a , T. , a n d Ao k i , T.

[13] Kisimoto, K. “Combined bathymetric and topographic mesh
data: Japan2 50 m.grd”, Geological Survey of Japan, Openfile Report, No. 353 (1999).

“Large-scale simulation of seismic-wave propagation of[14] Fujiwara, H., Kawai, S., Aoi, S., Morikawa, N., Senna, S., Kudo, N.,

the 2 011 Tohoku-Oki M 9 earthquake”, Proceedings of theOoi, M., Hao, K. X.-S., Hayakawa, Y., Toyama, N., Matsuyama,International Symposium on Engineering Lessons Learned

H., Iwamoto, K., Suzuki, H. and Liu, Y. “A study on subsurface

from the 2011 Great East Japan Earthquake, pp. 349-360 (2012).

structure model for deep sedimentary layers of Japan for[3] Okamoto, T., Takenaka, H., Nakamura, T. and Aoki, T. “GPU-strong-motion evaluation”, Technical Note of the National

accelerated simulation of seismic wave propagation”, in

Research Institute for Earth Science and Disaster Prevention,

GPU Solutions to Multi -scale Problems in Science andEngineering, Yuen, D., Wang, J., Johnsson, L., Chi, C.-H., Shi, Y.
(Eds.), 250 pp., Springer, in press.

No.337 (2009).
[15] Baba, T., Ito, A., Kaneda, Y., Hayakawa, T. and Furumura, T. “3-Dseismic wave velocity structures in the Nankai and Japan[4] Graves, R. W. “Simulating seismic wave propagation in 3 D

Trench subduc tion zones derived from marine seismic

elastic media using staggered-grid finite differences”, Bull.

surveys”, Abstr. Japan Geoscience Union Meet., S111- 0 0 6,

Seism. Soc. Am., vol. 86, no. 4, pp. 1091-1106 (1996).

Makuhari, Japan (2006).

[5] Ye e, K . “ N u m e r i c a l s o l u t i o n o f i n i t i a l b o u n d a r y v a l u e[16] Nakamura, T., Okamoto, T., Sugioka, H., Ishihara, Y., Ito, A.,problems involving maxwell’s equations in isotropic media”,

Obana, K., Kodaira, S., Suetsugu, D., Kinoshita, M., Fukao, Y.

IEEE Transactions on Antennas and Propagation”, vol. 14, no.

and Kaneda, Y. “3D FDM simulation for very-low-frequency

3, pp. 302-307 (1966).

earthquakes off Kii Peninsula”, Abstr. Seism. Soc. Japan, P1-[6] Abdelkhalek, R., Calandra, H., Coulaud, O., Roman, J. and Latu, G.

06, Hiroshima, Japan (2010).

“Fast seismic modeling and reverse time migration on a GPU[17] Kurahashi, S. and Irikura, K. “Source model for generatingcluster”, International Conference on High Per formancestrong ground motions during the 2011 off the Pacific coast

Computing & Simulation, pp. 36-43 (2009).

of Tohoku earthquake”, Earth Planets Space, Vol. 63, 571-576[7] Micikevicius, P. “3D finite-difference computation on GPUs

(2011).

using CUDA”, GPGPU-2 : Proc. 2 nd Work shop on General
Purpose Processing on Graphics Processing Units, pp. 79 84, Washington DC, USA (2009).
[8] M i c h éa , D. a n d Ko m at i t s c h , D. “Acce l e r at i n g a t h re e dimensional finite-difference wave propagation code using
GPU graphics cards”, Geophys. J. Int., doi: 10 .1111/j.13 6 5 246X.2010.04616.x (2010).
[9] 岡元 太 郎 , 竹中博士 , 中村 武 史 . “GPUによる地 震 波伝 播シミュ
レーション”, 先進的計算基盤システムシンポジウム（SACSIS2010）
論文 集 , pp. 141-142 (2010).
[10] 下 川 辺 隆 史・青 木 尊 之 , “ 次 世 代 気 象 モデル のフル GPU 計
算‐TSUBAME 2.0 の 39 9 0 GPU で 145 TFlops － ”, TSUBAME

16

● TSUBAME e-Science Journal No.6
2012 年 7 月31 日 東京工業大学 学術国際情報センター発行 ©
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/

17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/をご覧ください。