RISTニュース No．
５８（２０１５）

6都市全域の地震等自然災害シミュレションに関する研究*hp130015+7の高度化支援作業結果報告
Repor
to
ft
heRe
sultsofProgram Tuni
ngSupportto(ResearchonSi
mulat
ionso
fEar
t
hquakesandOtherNat
ur
al
Di
s
aster
sint
heEnt
ir
eCity(hp130015))
一般財団法人高度情報科学技術研究機構
志澤 由久、小林寛
RI
ST神戸センター利用支援部では、「京」利用支援の一つとして、東京大学地震研究所の課題「都市全域の地震等自然災害シミュレーションに関する研究（hp
１３００１５）」（HPCI
戦略プログラム分野３の研究開発課題）の高度化支援を実施した。本支援を反映した研究成果は、２０１
４年のゴードン・ベル賞ファイナリスト（最終候補）に選出された。本稿はその高度化支援結果について報告する。
１．はじめに

２．プログラム概要

RI
ST神戸センター利用支援部では、「京」プログラムで採用されているアルゴリズム

利用支援の一つとして、２０１３年度、東京大学の詳細は文献［１］に記述されている。ここ

地震研究所の課題「都市全域の地震等自然災ではプログラムの概要のみを記す。

害シミュレーションに関する研究（hp
１３００１５）」３次元の非線形波動場が計算対象である。

（HPCI
戦略プログラム分野３の研究開発課

変位ベースの四面体二次要素の有限要素法に題）の高度化支援を実施した。本支援を反映より解く。連立一次方程式のソルバーは共役した研究成果［１］は、２０１４年のゴードン・勾配法を使用する。剛性マトリクスと変位ベ

ベル賞ファイナリスト（最終候補）に選出さクトルとの行列ベクトル積は、Element
by-れた。ゴードン・ベル賞はHPC分野で最も権Element法で計算する。ソースはFortr
an７７

威ある賞の一つであり、高い実行性能と計算（約５０００行）で記述されている。ノード間は科学技術への成果に対してACM（Associa-MPIを、ノード内はOpenMPを用いて、ハイtionf
orComputingMachinery）が授与する

ブリッド並列化されている。

賞である。本支援により実行性能の向上が認
３．高度化支援の依頼内容められ、そのことがファイナリスト選出に少なからず貢献できたと考えている。

課題担当者から受けた支援依頼の内容は以以下に高度化支援作業結果の詳細について下のとおりである。解きたい問題規模はフル記述する。まず、プログラムの概要についてノードで１０００億自由度であり、必要な時間ス説明する。次に、課題担当者から受けた支援テップ数は１０万ステップである。現状では、依頼の内容を説明する。それに基づき作業目

昨年度、富士通殿が単体チューニング、及び的と作業方針を設定し、それに沿って実施し６４ノード位までのMPIチューニングを実施た作業内容の詳細を報告する。最後にまとめした。その結果では、１０００ノード位まではスをする。

ケールするが、４０００ノード以降で不穏そうである。そこで、MPI及びI/O等のチューニン−3−

RISTニュース No．
５８（２０１５）

グを行うことで、さらなる高速化と十分なス問題箇所の原因を探るのに使用。

ケーリングの実現を依頼された。
６．使用データ
４．今回の作業の目的および作業方針

プログラムを実行する際の入力データに

以上のような依頼に基づいて、今回の支援

は、課題担当者より提供頂いたデータを使用の目的、及び方針を、「並列性能の向上」とした。使用したデータを表１に示す。上二つ「単体性能の向上」に設定した。

はストロングスケールを調べるためのデータ

並列性能の向上のために、通信、演算、I/Oで、下三つはウィークスケールを調べるための経過時間を処理ブロック毎に計測し、スケのデータである。

ーラビリティを調査した。そして、阻害箇所

ストロングスケールのデータ（領域数：があれば、原因を特定し、改善策を検討する

６４，２５６）は、全体で１億自由度、６４ノードのこととした。

データを基点として、同じ自由度で、分割数単体性能の向上のために、まず性能解析ツ＝ノード数が４倍だけ違う関係にある。

ールを用いてコスト分布を調査した。そして

ウ ィ ー ク ス ケ ー ル の デ ー タ（領 域 数：高コストな箇所を特定し、単体性能向上策を

２
５６，４０９６，１６３８４）は、１億自由度２５６ノード試行することとした。

を基点として、自由度とノード数が同じ倍率を保って、１６倍、６４倍の関係にある。

５．プロファイラについて
表１「京」で提供されている幾つかの性能解析ツ使用データール（以下プロファイラ）を用い、性能分析を行った。各プロファイラの特徴を以下に記述する。
煙基本プロファイラ（fipp）
アプリケーション全体、あるいは指定された区間の基本プロファイラ情報（アプリケーションのハードウェアモニタ情報「京」で実行する際のラージページの設定やコスト情報等）を収集出来る。

は、以後、全て４MBとした。基本プロファイ

一定間隔でサンプリングを行い、その結ラ使用によるオーバーヘッドは無視できるこ果をコスト情報として出力する。と、及び計測したい箇所にタイマーを挿入し比較的少ないオーバーヘッドで容易に基たことによるオーバーヘッドは無視できるこ本プロファイラ情報を取得する事が可能。とを確認した。時間ステップ数は１００とした。

コスト分布を参照することで、負荷の高７．基本プロファイラによるコスト分布調い箇所を特定可能。
煙精密PA（PA：Perfor
manceAnalysis）

査：asis版、２５６ノード指定された区間の詳細なハードウェアモ

asis版（オリジナル版）について基本プロニタ情報（キャッシュミス率や演算数等）ファイラを用いてコスト分布を調査し、高コを収集できる。

ストな箇所を特定する作業を実施した。以下

ハードウェアカウンタ情報取得の都合、にその結果を報告する。ノード数２５６のデー７回プログラムを実行しなければならな

タを使用し、コンパイルオプションは−Kf
ast，い。

ope
nmpとし、１ノード当り８コアを使用した。

−4−RISTニュース No．
５８（２０１５）

図１

基本プロファイラを用いたコスト分布調査結果。一つのプロセスにおける、スレッドゼロでのコスト分布。

図２

基本プロファイラを用いたコスト分布調査結果。図１と同じプロセスにおける８個のスレッド間のコストバランス。コスト２位はスレッド非並列箇所であることが分かった。

図１の円グラフは、一つのプロセスにおけ応している。棒グラフから、コスト２位の赤る、スレッドゼロでのコスト分布である。

い部分はスレッドゼロにのみ存在していて、

コ ス ト １ 位 は１１％ を 占 め て い て、そ れ はスレッド非並列箇所であることが分かった。

nont
et１０kusという名前のサブルーチンであ８．プログラムの構造の概要ることが分かった。ソースコードを見ると、
これはEl
ement
byElement法で行列ベクト

これらの高コスト箇所がプログラムの中でル積を計算しているサブルーチンであることどのような位置にあるのかを把握するためにが分かった。コスト２位は１０％を占めていプログラム構造を調査した。プログラム構造て、それはCG法ソルバーのサブルーチン、図を図３に示す。プログラムは、まずグリッCG４内の部分であることが分かった。これがドデータを読み込み、時間ステップループにどのような部分なのかを、図２の棒グラフで入り、その中で計算を行い終了するという構表した。この棒グラフは、図１と同じプロセ造をしている。時間ステップループの中で、スにおける８個のスレッド間のコストバラン

時間ステップごとの入力データを読み込み、スを示している。図１の円グラフは、図２のCG法のソルバーで解いて、計算結果を時間棒グラフにおける、スレッドゼロの部分に対ステップごとにファイルに出力する。入力は−5−

RISTニュース No．
５８（２０１５）

ループの外で一度だけファイルオープンし、取得したデータについて計算する処理に対応

一つのプロセスは一つのファイルから次々としている。

読んでいく。これに対して出力のほうは、ループ内でその都度ステップ毎の別ファイル名をオープンして書き出しクローズしている。

図５

行列ベクトル積の構造概略。サブルーチ
ンnont
et
１０kus
はコスト１位箇所に対応する。

９．並列性能向上のための改善作業
改善作業結果の報告として、まず、並列性能向上のための改善作業の報告から入ることとする。
図３本プログラムの場合、基本プロファイラにプログラム構造図よる計測ではI/Oのコストは計測されない。
CG法のソルバーは、図４のような計算でI/Oの 時 間 も 把 握 す る た め、asis版 の 時 間構成されている。これらのうち、各領域内でス テ ッ プ ル ー プ 内 に つ い て、タ イ マ ー

ベクトルの内積計算する処理と、ベクトルのMPI_Wtimeを挿入し、通信、演算、及びI/O

スカラー倍や、和・差の計算にスレッド非並の時間を個別に計測することとした。さらに列部分がある。これらの部分の合計がコストその結果からスケーラビリティを調査した。

２位箇所に対応している。

阻害箇所があれば、それらを特定する作業を実施した。スケーラビリティの調査結果を
図６に示す。横軸がノード数であり、６４から２５６ノードがストロングスケール、２５６から１万６０００ノードがウィークスケールデータである。縦軸が経過時間である。青の実線が時間ステップループ合計の経過時間で、点線がその内訳である。

図４

赤の点線の演算部分について見てみると、

CG法のソルバーを構成する処理の種類。
赤い点線で囲った部分の合計はコスト２位箇所に対応する。

６４から２５６ノードで時間が約１/４になっており、ストロングスケーラビリティは良好であるといえる。２５６ノード以降で時間はほぼ一定であり、ウィークスケーラビリティも良好

行列ベクトル積について、構造の概略を図であるといえる。

５に示す。まずスレッドのループがあり、そ対して、ファイル出力がスケーラビリティの内側に各スレッドで担当する四面体要素のを阻害していることがみてとれる。

ループがある。コスト１位箇所のサブルーチ

また、asis版には、もともとMPI_Barrier

ンnont
et１０kusは、一つの四面体要素として−6−

RISTニュース No．
５８（２０１５）が通信ルーチンの前後に挟まれており、１スの書き出しでのサイズは４６％に削減された。

テップ当り約１４００回と数多くコールされてい書き出しをバイナリにすることは、テキストることが分かった。これらのバリアをコール変換の処理が無くなることを意味しており、する時間も計測対象としたところ、バリアもその効果による経過時間の短縮が期待できる。

スケーラビリティを阻害していることが分かった。さらに、バリアは計算ロジックに何

表２ ファイル出力の改善に関する検討結果の影響も無く、削除できることを確認した。
以上の考察から、ファイル出力の改善とバリアの除去を行った。

ファイル出力の改善効果の計測結果
ファイル出力の時間を計測した結果を図７に示す。左のグラフが縦軸を経過時間としたグラフである。ファイル出力の時間は無視できる程度にまで削減された。右のグラフはス
図６スケーラビリティの調査結果ケーラビリティを見るためのグラフである。
縦軸は２５６ノードでの時間を分母として時間９．
１ ファイル出力の改善を規格化したものである。スケーラビリティ

ファイル出力の改善効果の検討結果

はストロング、ウィークともにほぼ理想値へ（ファイルサイズの観点）と改善したことがみてとれる。

まず、ファイル出力の改善作業結果から報告する。表２にファイル出力の改善に関する検討結果、及びasisからtune版への改善点を示す。
asis版では各ランクはステップ毎に別ファイルに出力するという特徴がある。従って、ステップ数が多くなればなるほど出力ファイル数も多くなる。この方法の弱点は、より長図７いステップ数だとステージアウトのファイル数の上限値を超えてしまうということである。

ファイル出力の改善効果の計測結果（経過時間及びスケーラビリティ）

そこで改善方法として、各ランクは時間ステップを通して同一のファイルに出力するよ

次にランク間のばらつきの計測結果を表３うにした（表２の上側の赤い点線部分）。さらに示す。ばらつきの時間を「ランク間の最大に、テキストで書き出していたものを、tune−最小」と定義すると、オープン・クローズ版ではバイナリで書き出すように改めた（表やテキスト変換処理がなくなったことによ

２ の 下 側 の 赤 い 点 線 部 分）。こ れ に よ り、り、特に４０００ノード以上で、ばらつきの時間ファイルサイズの削減が見込まれ、実際一回が大幅に減少した。

−7−

RISTニュース No．
５８（２０１５）

表３

ファイル出力の改善効果の計測結果（出力のばらつき）

図８

バリア除去効果の確認結果

９．
３ バリア除去＋ファイル出力改善による

時間ステップを通して同一のファイルに出総合的な効果の確認

力する改善手法、及びバイナリで書き出すようにする改善手法は、本アプリケーションプバリア除去とファイル出力改善による総合ログラムに限らず、「京」で性能向上させるた的な効果を時間ステップループについて確認

めの必須の改善手法であるといえる。

する作業を実施した。結果を図９に示す。経過時間については、４０９６ノードの時、１．
５倍の９．
２ バリアの除去高速化を達成した。２５６ノードの時間を分母

バリアの有無に関する注意点にして規格化し、スケーラビリティをみると、ストロング、ウィークともに改善した。

バリアの除去結果について報告する前に、まず注意点にふれる。asis版では、バリアは１００ステップで約１４万回コールされている。
asis版では、バリアは通信ルーチンの前後に挟まれている。それ故、計測されたバリアの時間の中には、バリアの純粋な時間に加えて、演算やI/Oのロードインバランスに起因する、同期待ち時間も計上されている。一方、バリア除去したtune版では、同期待ち時間は、allr
educe
や隣接通信の計測時間の中図９に計上されている。バリアを除去することで、少なくともバリアを１００ステップで１４万バリア除去とファイル出力改善による総合的な効果の確認結果

回コールする際に発生するバリアの純粋な時
１０．単体性能向上のための改善作業間の削減が期待できる。

次に、単体性能の改善作業の結果報告に移バリア除去の効果の確認る。ここでは実際に行ったチューニングのう効果を確認するため、バリアとallreduceち、コスト１位：nontet
１０kus
ルーチンの高と隣接通信の時間を合計した値を比較するこ

速化、及びコスト２位：CG４ルーチンのスととした。結果を図８に示す。グラフより、

レッド非並列部分のスレッド並列化について確かに削減効果が得られたことが確認できた。

報告する。

特に４０９６ノードの時２３．７４秒短縮である。
et１０kus
ルーチンの高１０．
１ コスト１位：nont
速化（ElementbyEl
ement
法による行列ベクトル積計算部分）
まず、コスト１位ルーチンの高速化につい−8−

RISTニュース No．
５８（２０１５）て報告する。このサブルーチンは、ひとつのンの高コスト箇所を調査した。各実行文につ四面体二次要素（自由度３０）について、要素いて、コストの高い順に示したのが図１１であ剛性マトリクスを生成して、節点ベクトルる。

（配列要素数３０）に掛ける処理を行う、という一番コストが高いのは、計算結果を格納すものである。計算に用いる１個の四面体デーる配列BDBuのゼロクリアであり、全体のタは小さいのに対して、計算は複雑で多いた

１０％ぐらいを占めていることが分かった。配め、要求B/Fは低い。要素数回コールされる

eと小さいが、他列BDBuのサイズは１２０bytため、高コストである。の積和演算はL１オンキャッシュのため、メモリアクセスが必要なこの部分が相対的に高コ

nont
et
１０k
usas
is
版〜nontet
１０k
usを切り出しストである。

たテストプログラムの性能分析〜

積和演算部分が２％前後で、その後に続いサブルーチンを切り出したテストプログラている。これは、積和演算はL１オンキャッムを作成（６４ノードデータ使用）して、性能シュだが、複雑な計算が多数存在しているた分析と改善策を試行する方法をとった。まめである。

ず、asis版の精密PAの結果を図１０に示す。このサブルーチンは１個の四面体に関する計算である。初回使用時にメモリから取得した計算に必要なデータは、その後L１キャシュ上に全て乗る問題である（メモリスループットの測定値はゼロである）。演算性能ピーク比は１９．
８９％である。

図１１ nont
et１０kusasis版の基本プロファイラによる性能分析結果

nont
et
１０kus
チ ュ ー ニ ン グ：配 列BDBuの ゼ
ロクリア削除
配列BDBuには計算結果を次々と足し込む処理が行われていて、ゼロクリアは不要であることが確認できた。そこで、ゼロクリアを
図１０ nont
et１０kusasis
版 を切り出したテストプログラムの性能分析（６４ノードデータ使 用）。nont
et１０kusの コ ー ル 回 数 は６４と し た。１コ ー ル 当 り の 経 過 時 間 は０．
１２１秒である。

コメントアウトし、BDBuへの一回目の計算
値代入に関して、BDBuの加算処理を削除した（図１２参照）。

基本プロファイラによるnont
et１０k
usasis版の性能分析
これ以上の高速化が可能かどうか検討するために、基本プロファイラでこのサブルーチ−9−

RISTニュース No．
５８（２０１５）

図１２ nontet１０kus
チ ュ ー ニ ン グ 例１：配 列
BDBuのゼロクリア削除

図１４ nont
et１０kus
チューニング例２：重複計算除去（続）

nont
et
１０kus
チューニング： 重複計算除去

nont
et
１０k
usチ ュ ー ニ ン グ 版 〜nont
et
１０k
us

何度も現れる計算パターンについて、一時を切り出したテストプログラムの性能分析〜変数に格納し、それを使いまわすことで演算チューニング後のコードに対する精密PA数を削減した。その例を図１３と図１４に示す。

を 計 測 し た 結 果 を 図１５に 示 す。演 算 数 は

図１
３一行目におけるc（２０）*．
．．という式は８．
１％減少を確認した。演算性能ピーク比は何度も現れるので、左辺の緑色の一時変数に２１％に向上した。経過時間は１５．
８％短縮した。

格納して、図１３の３行目の右辺で使うようにチューニング後の状態で、命令コミット系した。の時間が減少しており、また、各命令コミッ

さらに図１３の３〜４行目にかけての赤及びトの占める割合に、若干の変化が見られる。

緑字で記した長い式も、図１４の上側の四角内

このような現象が、演算数の変化率と経過時で赤字で記した様に同じパターンで何度も現間の変化率が必ずしも一致しない原因の一つれるので、図１３の２行目の左辺のように、変なのかもしれないと推察している。

数名は長いが一時変数に格納して、これらの計算式が現れる箇所で、図１４の下側の四角内の赤字部分の様に置き換えて使いまわした。
このようにすることで演算数を削減した。

図１３ nontet
１０kus
チ ュ ー ニ ン グ 例２：重 複 計算除去
図１５ チューニング前後のコードに対する精密PA計測結果の比較

１０．
２ コスト２位：CG４スレッド非並列部分のスレッド並列化及びその他チューニ
ング（例：内積の計算）
次に、コスト２位のCG４スレッド非並列部分のOpe
nMPスレッド並列化チューニング

−10−

RISTニュース No．
５８（２０１５）について報告する。ここでは、例として、二つの節点ベクトルの内積計算のコードを取り上げる。このコードには、節点内自由度のインデックス計算の除去、及び内積計算をスレッド並列化する際にOMPREDUCTI
ONを用いない手法も適用できるので、合わせて適用した結果を示す。
asis版 の 内 積 計 算 コ ー ド を 図１６に 示 す。

図１７ 内積計算のコード（チューニング版）

asis版のコードでは、まず節点のループがあスレッド並列化、及びREDUCTIONの回避り、その内側に節点当り３つの自由度に関するループがある。節点ベクトルｚとｑを掛け次に、チューニング版のコード（図１７）にて、さらに重みファクター wを掛けて、zqにおける、スレッド並列化及びREDUCTION足し込んで、内積を求めている。その際、節の回避の仕方について説明する。スレッド番

点番号と節点内自由度の番号から配列ｚに対号myidのループと、各スレッドでの担当範するインデックスを図１６の赤い点線で囲まれ囲内での節点iのループ、というように二段た部分の引数の式のように計算している。

階 に す る。外 側 の ス レ ッ ド の ル ー プ をOpenMPでスレッド並列化する。内側の節点ループについて、各スレッドが担当する範囲の下限・上限を明示的に前もって作成し、配列n_istaとn_iendに格納したものを参照して下限・上限を指定する。スレッドmyidでの部分和は、一旦、配列要素zq_dum（myid）

図１６ 内積計算のコード（asis版）に 格 納 す る。配 列zq_dumの 配 列 要 素 を ス
節点内自由度のインデックス計算の除去

レッド番号順に和をとる。このようにするこチューニング版のコードを図１７に示す。

ON指 示 文 を 使 用 せ ず に スと でREDUCTI

チューニング版では、赤い点線で囲まれた部レッド並列化できる。この方法の利点は実行分のように配列ｚを二次元配列に変更し、一形式（逐次／スレッド並列）によらず常に同次元目は節点内自由度のインデックス、二次一な計算結果が得られることである。

元目は節点番号とすることで、余計なイン参考として、スレッド番号順に和をとる別デックス計算が不要となる。さらに、重み法として、REDUCTIONを使用し、かつ、コンファクターｗは節点番号のみに依存すること

パイルオプションKordered_o
mp_reductionに着目して、節点当り三つの自由度のループを用いる方法もあるということに触れておく。

を展開し、まとめて括ってからｗを掛けるこチューニングの効果の確認とで、掛け算の数を減らすことができる。

処理内容と適用した手法との対応は表４のようになる。これらの処理内容の合計に相当する部分（通信含まず）について、チューニング前後で経過時間を比較した結果を図１８に示す。６４ノードで約１．
８倍の高速化を達成した。理想値１/８に届かない理由は、ループ−11−

RISTニュース No．
５８（２０１５）びメモリスループットについて報告する。

内の演算が少なくスレッド並列化のオーバーヘッドが大きいためである。

１１．
１ 最終チューニング版の計測結果の要約
表４

処理内容と適用した手法との対応最終チューニング版の時間ステップループ（I/Oと通信を含む）部分の各種計測結果を図１９に示す。経過時間については、１万６０００ノードの時、３．
２倍の高速化を達成した。スケーラビリティについても、ストロング、ウィークともに改善した。チューニングにより、演算性能、メモリスループットともに向上した。特に、１万６０００ノードの時の演算性能の向上率は２．
９倍、メモリスループットの向上率は３．
２倍となった。
１２．まとめ
並列性能の向上作業として、バリアの除

図１８ 表４の処理内容の合計部分（通信含まず）の経過時間の比較結果

去、及びファイル出力の改善を実施した。
ファイル出力の時間は、ほぼ無視できる程度にまで削減され、スケーラビリティも改善し１１．総合的な改善効果の確認た。また、時間ステップループのスケーラビリティも改善した。

以上の単体性能と並列性能向上策の総合的な効果の確認のための計測を時間ステップル

単体性能の向上作業として、高コストルーープに対して行った。その計測結果として、チンのチューニング、及び各種チューニング

経過時間、スケーラビリティ、演算性能、及を実施した。コスト１位nont
et１０kusについ図１
９ 最終チューニング版の計測結果の要約（時間ステップループ部分の性能）−12−

RISTニュース No．
５８（２０１５）ては、６４ノードの時、１．
１９倍高速化し、演算なお、本高度化申請後のチューニングによ

性能ピーク比は２１．
６８％に向上した。コスト

り、本支援を反映した研究成果［１］では、２位CG４スレッド非並列部分は、６４ノードの２７０億自由度の問題に対して、京コンピュー時、１．
８倍高速化した。

タ全系（８２９４４ノード）で演算性能ピーク比
６．
９３％ （０．
７３６ PFLOPS）の実行性能を出す

並列性能と単体性能の向上策を併せて実施したことで、大幅な改善が達成された。１万に至っている。

６０００ノードでのステップ数１００の時間ステッ
プループ（I/Oと通信を含む）部分について、参考文献asisと最終チューニング版の比較を表５に示［１］Tsuyoshi Ichi
mura, Ko
hei Fujita,す。経過時間は３．
２倍の高速化を達成し、演SeizoTanaka,MuneoHo
ri,Maddege-

算性能ピーク比は２．
２％から６．
４％に２．
９倍向daraLalith,Yoshi
hisaShizawa,and

上し、メモリスループットピーク比は ３．
２％bayashi: Physicsbas
edHiroshi Koから１０．
４％に３．
２倍向上した。

urbaneart
hquakesi
mulationenhanced by 10.
7 BlnDOF x 30 Kti
mestepunstructuredFE nonlinear

表５ ノード数１６３８４での時間ステップループ（I/Oと通信を含む）部分の性能比較
バージョン
asis版

seismicwaves
imulation,Proceedi
ngs

経過時間 演算性能 メモリスルー［sec］ ピーク比 プットピーク比
３
７
２．
４
０

２．
２％

最終チュー
１
１
７．
２
２
６．
４％
ニング版
（３．
２倍） （２．
９倍）
（向上率）

oft
he2014ACM/I
EEE confer
enceon
Super
co
mput
ing (SC 14)
,(ACM

３．
２％

Gor
don BellPrizeFi
nalist),2014,inpress.

１
０．
４％
（３．
２倍）

−13−