京における最適化手法の調査研究Research study of code optimization methods on K computer

hp120031 「京」共用法第 12 条調査研究 Investigation Study Defined by the Public Utilization Act

京における最適化手法の調査研究

Research study of code optimization methods on K computer

山岸孝輝 Takateru Yamagishi 高度情報科学技術研究機構

Research Organization for Information Science and Technology

要旨 本研究の目的は、実アプリケーション、その中でも数値気象モデルの京における最適化手法の 調査研究を行い、京の利用支援業務に資することである。最新の需要に対応すべく、新しい格子 構造の気象モデルを研究対象として選択した。性能評価とその結果解析から、格子点データの再 構築が性能向上への有効な手段と判断した。新しい格子構造の気象モデルであるが故に再構築に は想定以上の作業時間が必要であり、当初予定していた複数の数値スキームの最適化と考察にま で至らなかったが、京の利用支援業務に有用な新しい格子構造の気象モデルを構成する技術の基 礎を把握することが出来た。 キーワード：数値気象モデル、数値予報、数値流体力学、最適化、ステンシル計算 Abstract

To acquire the knowledge for user support on K computer, we tried to study and develop the methods which are applicable to the code optimization of the numerical weather models on K computer. We aimed to meet the latest need in this present research community and chose the numerical weather model with new kind of grid system. We found that reordering of grids would be effective to improve the performance of the model after the survey and analysis of several model executions. Although we could not complete the reordering due to its complexity, this research study will be useful and instructive for our user support on K computer.

Keywords: Numerical Weather Model, Numerical Forecasting, Numerical Fluid Dynamics, Code Optimization, Stencil Calculation

© 2016 Research Organization for Information Science and Technology All rights reserved. Received: 2 December 2015

Accepted: 20 June 2016 Available online: 30 June 2016

1. 研究の背景と目的 本研究の目的は、実アプリケーション、その中でも数値気象モデルの京における最適化手法の 調査研究を行い、京の利用支援に資することである。実施に当たっては利用支援業務において最 新の需要に対応できる体制作りを視野に入れて、従来のモデルに比べて新しい格子構造の気象モ デルを研究の対象として選択した。 2. 計算モデル 数値モデルの選定に当たっては計算科学且つ気象・気候シミュレーションの研究としての観点 に加えて、利用支援業務のための研究であることも考慮した。具体的には、利用支援業務と並行 して行う中での限られた作業時間でも実施出来ること、他の利用者の計算資源を圧迫しないこと に加えて、利用者の今後の支援要請に応えうるような最新の研究動向に沿った数値モデルである ことが必要と考えた。 以上から、準備が容易であることに加え代表的な実験 1 実行あたりの CPU 利用時間が少なく て済む数値モデルを選定の対象とすると共に、本研究に割り当てるノード時間も 360 時間という 他研究と比較して非常に少ないもので申請した。 この半期では CPU 単体での性能に関する調査研究を 1 スレッドまたは 1CPU 内スレッド並列 の実行にて行い、今後の並列数の多い条件への拡張へとつなげることを実施目標とした。その目 標を達成するには、1 スレッドまたは 1CPU 内スレッド並列の実行でも、計算科学並びに気象・ 気候シミュレーションの研究として意味のある実験が利用支援業務に支障をきたさない程度の 実時間内に終了する事に加え、検証のために広く使われている実行テスト設定が存在して先行研 究との比較が容易な数値モデルであることが必要であった。 一般に気象・気候モデルは力学・物理など複数の過程を含む。それら複数の過程を全て数値モ デルに含むと、タイムステップあたりの経過時間が増加することに加えて計算結果が気象・気候 学として意味のあるものとなるまでに長い経過時間を必要とする。また、複数の過程が相互作用 するために数値モデルの非線形性が強くなり、結果のばらつきが大きくなるという問題がある。 よって、本研究では基本的な力学過程のみで構成される数値モデルを用いる事とした。構成要 素を少なくすることで計算時間が短くなり、先行事例が豊富であるため先行研究との比較を通じ て、計算科学からの詳細な考察が可能になる。また、力学過程は気象・気候シミュレーション数 値モデルの核となる部分でもある。 本研究では、以上の条件に加え、最新の研究動向に合うものとして、正 20 面体格子を採用し た ICON 浅水波モデル（参考文献 1）を用いることとした。

3. 並列計算の方法と効果（性能） 用いた並列計算の方法は、コンパイラによるスレッド間の自動並列化である。以下、基本・詳 細プロファイラを用いた性能評価結果並びにコンパイル時におけるメッセージ（コンパイルリス ト）の調査結果をまとめる。 数値モデルの全体を効率よく概観することを目的として、基本プロファイラを用いて実行時間、 演算量、メモリスループット、SIMD 化率らを計測した。後述する詳細プロファイラに比べて、 オーバヘッドが少なく、実行時間に大きな変化を与えにくいことに加えて 1 回の実行でプロファ イルの取得を完了出来るという特徴がある。実際の実行をよく反映するプロファイルを少ない作 業量で取得することができるため、数値モデルの特徴の概観を捉え、その上で詳細プロファイラ の対象とする箇所を決めるのに有効な手段である。 表 1 1 または 8 スレッド実行時の基本性能 Time[s] MFLOPS MFLOPS(理 論性能ピー ク比)[%] メモリスルー プット[GB/s] メモリスループ ット(理論性能 ピーク比)[%] SIMD[%] 1 スレッド 154.0 191.2 1.20 0.32 0.50 21.08 8 スレッド 平均 135.1 217.5 0.17 0.37 0.58 16.19 0 135.1 97.9 0.61 0.37 0.58 8.50 1 135.0 17.1 0.11 0.37 0.58 28.14 2 135.0 17.1 0.11 0.37 0.58 28.15 3 135.0 17.1 0.11 0.37 0.58 28.14 4 135.0 17.1 0.11 0.37 0.58 28.14 5 135.0 17.1 0.11 0.37 0.58 28.13 6 135.0 17.1 0.11 0.37 0.58 28.14 7 135.0 17.1 0.11 0.37 0.58 28.04 FLOPS の実行効率が 1%程度、メモリスループットのピーク比が 0.5%程度であることに加えて SIMD 化率が 21%程度と、ハードウェアの性能を生かし切れていない可能性が高い（表 1）。 8 スレッド実行は、演算、メモリ並びに SIMD 化率が 1 スレッド実行同様に小さい値となって おり、加えてスレッド間のインバランスが存在し、その中でもスレッド 0 の負荷が大きい（表 1）。 スレッド 1 からスレッド 7 はほぼ同じプロファイルであることから、スレッド並列化されず 1 ス レッドで実行されるループが存在していると考えられる。 以上から、1 スレッド実行にてなんらかの性能阻害要因がある事、8 スレッド実行にてスレッ ド間のインバランスがある事がわかった。まずは 1 スレッド実行での性能阻害要因を調べるため に、1 スレッド実行にて最もコストが高いループ（loop1 とする）に対して詳細プロファイラを用 いた計測と解析を行った。

図 1 1 スレッド実行にて最も負荷が高いループ（loop1）の詳細プロファイル 浮動小数点並びに整数のロードキャッシュアクセス待ちが実行時間の大半を占めることがわ かる（図 1）。L1 キャッシュに対してデータを要求した際のミス率（L1dm ミス率）が約 96%と、 L1 キャッシュが殆ど活用されていなかった。ループ長が長く L1 キャッシュにデータが収まりき らないか、L1 キャッシュの特定の領域のみが演算器からアクセスされ、特定の領域のみでデータ の書き換えが発生している可能性がある。L2 キャッシュのミス率は約 2.5%と、倍精度実数を一 様にアクセスした際に、理論的に想定されるキャッシュミス率よりも小さい値であった。L2 キャ ッシュについては有効に利用されていて、最適化を施す必要性は低いことがわかった。 メモリスループットは 0.18GB/s と、ピーク性能の 0.3%程度で、ハードウェアの性能を十分に 活用できていなかった。SIMD 化については、SIMD の対象となる演算命令の内、SIMD 演算が行 われたものが約 78%と改善の余地があることに加え、ロードとストアにおいては命令が全く SIMD 化されていなかった。 loop1 のコード（一部説明加筆）を以下に示す。 DO iv0=1,nvertx：[1:20480] DO j=1,dim_rbf_v：[1:9] temp(:)=temp(:)+rbf_vect_coeff(iv0,j,:)*hin(ivv(iv0,j)) ENDDO ENDDO 同箇所のコンパイルリスト（一部説明加筆）を以下に示す。 <<< Loop-information Start >>> <<< [OPTIMIZATION] <<< SIMD：ループが SIMD 化された <<< SOFTWARE PIPELINING：ループがソフトウェアパイプライニングされた <<< Loop-information End >>> 2 4v DO j=1,dim_rbf_v

<<< [OPTIMIZATION]

<<< FULL UNROLLING：配列式 temp(1:3)に対してアンローリングされた

<<< Loop-information End >>> 2 12v temp(:)=temp(:)+rbf_vect_coeff(iv0,j,:)*hin(ivv(iv0,j)) 2 4v ENDDO SIMD 化並びにソフトウェアパイプライニング化されており、演算は効率的に実行されること が期待される。しかしながら、変数 rbf_vect_coeff は最も内側のインデックス iv0 が 20480 ごとの ストライドアクセスを起こしており、このためにL1キャッシュミスが頻発していると判断した。 4. 研究成果 基本プロファイラ並びに詳細プロファイラを用い、階層的に性能評価を行った。取得したプロ ファイルから、L1 キャッシュミスが多いことがわかった。その結果メモリスループットは非常に 小さい値となった。また SIMD 化された演算の割合も低く、ハードウェアが本来持つ性能を十分 に活用できていなかった。 特に負荷が大きいループのコードとコンパイルリストから、ループの回転に従うデータへのア クセスが、ストライドアクセスとなっている事が、L1 キャッシュミスが頻発している原因と判断 した。 このように、各キャッシュを如何に利用しているかを具体的に追跡する事ができたのは、京に て利用出来る、詳細かつ豊富な種類のプロファイラがあってのことである。 負荷の高いループでのストライドアクセスを解消する方法として、本研究では格子インデック スの並べ替えによるデータ構造の再構築を選択した。 正 20 面体格子を採用した気象モデルは従来のモデルに比べて格子点の並びが特殊であり、格 子点データの構造を把握する作業が困難であった。それぞれの変数の並びを説明する文献等は無 く、またソースコードにも説明が記載されていなかった。そのため、離散化の原理のみを記述し た論文の内容を元に、ソースコードを読み、適宜 write 文の挿入・デバッガによる出力を行い、 いくつもの試行を繰り返しながらひとつひとつ変数の意味を確認していく作業に追われた。変数 それぞれの意味を確認することが困難であることに加えて、対象となる 26 個の変数は独立では 無いために 1 変数だけを並び替えてテスト実行することができなかった。そのため、別途描画ツ ールを用意して変数が想定通りに並び変わっているかを確認する作業が必要になり、半期という 限られた期間内に作業を完了出来なかった。 しかしながら、新しい格子構造の気象モデルについて調査ならびに最適化手法の検討を進めて きたことで、格子構造の気象モデルを構成する技術の基礎を把握することが出来たと考えている。

5. まとめと今後の課題 利用支援業務において最新の需要に対応することを想定し、新しい格子構造の気象モデルを研 究の対象として選択し、実アプリケーション、その中でも数値気象モデルの京における最適化手 法の調査研究を行った。 モデル構造の理解、実行時間に影響を与える要素の検討、性能プロファイルの取得と解析を進 めた。一部の限られたループが性能のボトルネックとなっていて、L1 キャッシュミスが多く発生 していることがわかった。該当部分のコード並びにコンパイルリストから、ストライドアクセス が原因と判断した。 本研究では格子インデックスの並べ替えによるデータ構造の再構築が有効な手段と判断し、京 並びに自所有のデスクトップ PC にて該当作業を利用支援業務と並行して行った。評価を行うに は、全てのデータを並び替える必要があるが、新しい格子構造の気象モデルであったが故に、想 定以上の作業時間が必要であることが判明した。 期間内にてその作業を終了出来ず、データ構造の再構築が実行時の性能に与える影響を評価す るまでには至らなかったが、京の利用支援業務に対して有用な新しい格子構造の気象モデルを構 成する技術の基礎を把握することが出来た。これまでに行った検討を次年度以降の利用支援業務 と並行して継続していく中で技術を高め利用支援業務にフィードバックすると共に、今後の業務 にて扱う他の気象・気候モデルの性能分析・改良の結果を基にメモリスループットや SIMD 化率 の様な代表的な性能指標について事例の蓄積と分析を行い、気象・気候モデルの開発と改良に広 く資するような知見の提供を目指していく予定である。 謝辞 本研究では特定先端大型研究施設の共用の促進に関する法律第 12 条の規定に基づき、特定高 速電子計算機施設のうち研究者等の共用に供する計算資源を、平成 24 年の下半期にて利用させ て頂きました。 参考文献

[1] L. Bonaventura, L. Kornblueh, M. Giorgetta, E. Roeckner, T. Heinze, D. Majewski, P. Ripodas, B. Ritter, T. Tingler, and J. Baudisch. The ICON shallow water model: scienti c documentation and benchmark tests., 2004.