プログラムの高速化・並列化サービスの事例

＊ここに掲載するのは、エジンバラ大学EPCCのSavvas Petrou博士によるHECToRレポート「dCSE Fluidity-ICOM: High Performance Computing Driven Software Development for Next-Generation Modelling of the Worlds Oceans, Xiaohu Guo1), Gerard Gorman2), Michael Lange2), Andrew Sunderland1), Mike Ashworth1), 1)Daresbury Laboratory, 2)Imperial College London」を要約したものです。

概要

Fluidity-ICOM $Imperial College Ocean Model$プロジェクトは、2010年10月から2012年31月まで行われました。Fluidity-ICOMはオープンソースです（https://launchpad.net/fluidity）。
Fluidity-ICOMは、非構造異方性アダプティブメッシュ上に様々な有限要素法と有限体積法を実装したオープンソースの偏微分方程式シミュレータで、様々な外部ソフトウェアパッケージと3種の言語（Fortran, C++, Python）が用いられています。
最深のマルチコアプロセッサを活用することは科学技術ソフトウェアにとって新たな課題を提示しています。そのアークテクチャには、複数コアを持つ複数プロセッサや、異なるレベルのキャッシュの共有、コアのサブセットに対する親和性を持つ複数のメモリーコントローラ、メインメモリーへの不均一なアクセス時間等です。

このため、ノード内ではスレッド並列、プロセス間通信にMPIを用いるハイブリッド並列に多くの関心が集まっています。こうした混合モード並列実装の利点には、第一に、純粋なMPI手法に比べてメモリー容量が少なく済むことが挙げられます。第二に、ノード内で必要なハロ領域を削除することでさらにメモリーが削減されます。例えば、全メッシュハロサイズはパーティション数（プロセス数）に比例して増加します。経験的には、線形四面体メッシュの頂点ハロサイズはパーティション数Pの1.5乗に比例します。最後に、ノードあたり1つのプロセスにI/Oを任せれば、各プロセスに各ファイルI/Oをさせるアプリケーションにとっては、大規模プロセス利用時におけるファイルシステム上のメタデータ操作を大きく削減することが可能です。つまり、ハイブリッドOpenMP/MPI並列手法は、計算ノードあたりのメモリー要求と、総ディスク出力データ量、およびメタデータの総操作量を削減します。

最新のマルチコアアーキテクチャマシンにおいて、ハイブリッドOepnMP/MPIは純粋な分散メモリー並列を超える数値アルゴリズムの最適化を提供します。例えば、代数的マルチグリッド法のスケール性は、領域境界を跨ったコースニングが困難なことからサブ領域数が増加すると悪化します。メッシュアダプティビティ手法のスケール性についても、領域境界を跨ったアダプテーションが出来ないため悪影響を受けます。

異なるシステム間の可搬性はアプリケーションパッケージに重要な性質です。指示行ベースの手法は可搬性を維持しながら並列処理を表現することが可能な方法です。OpenMPは、組み込みシステムやアクセラレータを含む如何なるシステムでも同じコードを用いることができます。これがFluidity-ICOMへOpenMP並列を実装する理由です。
しかしながら、そのループ並列の簡潔性にもかかわらず、スケール性を持つOpenMPの実装は自明な事ではありません。その競合条件や性能の落とし穴に注意が必要です。

このプロジェクトは、最初に有限要素アセンブリステージのOpenMP/MPI並列化、次に大規模コア数における線形プレコンディショナ/ソルバーのHYPERライブラリ活用による最適化、最後にスレッド化PETScライブラリを含むベンチマークを行います。

有限要素アセンブリステージのOpenMP/MPI並列化

事前の性能解析[2]から、シミュレーションのコストは、スパース行列アセンブリ（30-40%）とスパース線形系の求解の二つに占められていることが判っています。Fluidity-ICOMはPETScインターフェイスを通して、HYPERライブラリのハイブリッドスパース線形系ソルバー/プリコンディショナを用いることができます。これは純粋なMPI実装とも比較可能です。よって、残るOpenMP化が必要なのは、有限要素行列アセンブリカーネルです。これが高コストなのは、最深ループがサイズの4次で大きくなる深いループネスト、結合運動量や圧力や自由表面や各移流量といった構築すべき行列の多さ、間接参照アドレス（これは有限差分に比べてよく知られる有限要素法の欠点です）、キャッシュ再利用（非構造メッシュにとって特に重要な問題）といった事情によるものです。行列アセンブリコストは高次の増加を示し、不連続ガラーキン（DG）法が用いられます。

シミュレーションでは、例えば速度や圧力の連続、非連続有限要素形式や、Navier-Stokes方程式に対するトレーサー場やストークス流れなど、アセンブリ対象の行列が多く存在します。これらは個別にOpenMPで並列化すべきです。グローバル行列は、メッシュ（あるいは領域分割が適用されている場合はサブ領域）の全要素ループ内で、要素からの寄与を加算していくことで計算されます。スパース行列はPETScのCSR$Compressed sparse row$コンテナ（これには例えば流速ベクトルやDG場に用いるblock-CSRを含みます）に保存されます。要素の寄与は、CSRフォーマットで保存されているスパース行列へ加算されます。

·Greedyグラフ彩色法

アセンブリループは、要素からローカル行列を構築して、それをスレッドセーフにグローバル行列へ追加します。これはよく知られたグラフ彩色法[3]で実行できます。これは、最初にグラフ作成で実装されます。ここで、グラフのノードはメッシュ要素に対応し、グラフのエッジはグラフ構築における要素間の関係を定義します。ここで各色はそれぞれ、平行してグローバル行列へ加算される要素の独立なセットを定義します。こうしてデータの競合が避けられ、OpenMPのクリティカルセクションによる効率的な並列化が可能になります。

行列アセンブリの彩色法を用いた並列化では、メインアセンブリループの外側に色ループを作成します。要素に関するメインアセンブリループは並列指示行でOpenMP並列化されてループは指定されたサイズのループに分割されます。このループ内で、要素はその色が同じ場合にのみ行列へアセンブリされます。

·行列アセンブリカーネルの性能改善

ベンチマークはHECToR Cray XE6-Magny Cours $phase2b$およびCray XE6-interlagos $phase3$を用いました。テストケースは風成傾圧循環です。DG法による200M流速自由度、10M頂点メッシュまで用いました。基本とした実行はメッシュアダプティビティを用いない4時間ステップです。ここでは行列アセンブリと線形ソルバーを主に検討します。運動方程式と設定の詳細は文献[2]を参照してください。
運動量方程式カーネルは不連続ガラーキン（DG）法を用い、連続ガラーキン（CG）法は上述の足り方で並列化されます。

·ローカルアセンブリ

PETScでは、行列へ要素が追加される時スタッシュ$stash$が用いられます。並列行列フォーマットにおいてこれはある一つの重要な機能を提供します。それは、異なるプロセスにおいてローカル行列の一部として格納される要素を1つのプロセスに追加することができるものです。アセンブリの最中に、スタッシュされた値は正しいプロセッサーへ移動されます。これをノンローカルアセンブリを呼びますが、アセンブリループ内でスレッドセーフに関する問題を引き起こします。

ローカルとノンローカルアセンブリの性能差は、少数コア数では目立ちません。しかしながらコア数を大きくして768コアの場合には、ローカルアセンブリコードは40%高速です。
また、Criticalのような集団同期指示行をループ内で用いないようにすれば性能は50%以上向上する場合があります。

·メモリー帯域の最適化

ccNUMAノードの性能の鍵の一つはメモリー帯域です。メモリーバンド幅の最適化には以下の手法を用います：

1. 初期化にファーストタッチを用いると、CPUに直接接続されたページフォルトを発生させるメモリバンクにおいてページフォルトが生じることが保証されます
2. 個々のスレッドが計算中に同じコアにバインドされていることを保証するThread pinningを用いる。

Thread pinningは全ベンチマークを通して使用しました。ファーストタッチは12スレッドまでは効果がありましたが、それでも24スレッドでの急激な性能劣化は残ります。この問題は、メモリー確保が大きな寄与を占めています。

Thread-Caching malloc$TCMalloc$（http://goog-perftools.sourceforge.net/doc/tcmalloc.html）は、メモリ確保と解放にロック機構を用いない同期手法を提供します。これにより純粋なMPI版よりも純粋なOpenMP版のほうが優れた性能を示しました。

大規模コア利用時の線形プリコンディショナ/ソルバーに対するHYPERライブラリ利用の最適化

Fluidity-ICOMはスパース線形系の求解にPETScライブラリを用います。以前の解析から、実行時間の大半がスパース繰返し法線形ソルバーに占められていることが判っています。ここでは主に、スレッド化されたプレコンディショナとしてHYPERライブラリのBoomerAMGを調査します。

BoomerAMGはセットアップと求解という2つのフェーズから成ります。セットアップでの計算カーネルは、粗グリッドの選択と、内挿演算子の生成、各粗グリッド上のグリッド行列演算子の表現を行います。求解フェーズの計算カーネルは、行列ベクトル積（MatVec）とスムーシング演算子です。
MatVecや内積などの基本的な行列/ベクトル演算は、ループレベルでOpenMP並列化されています。セットアップでは粗グリッド演算子（3つの行列積）のみがスレッド化されています。コースニングと内挿はともにOpenMPを用いた並列化は用いません。求解フェーズは完全にスレッド化されています[9]。
使用したHYPER-2.9.1aはPETScのサイト（http://ftp.mcs.anl.gov/pub/petsc/externalpackages/hypre-2.9.1a.tar.gz）からダウンロード可能です。HYPRE BoomerAMGのマニュアルは[10]を参照してください。

PETScは、線形代数に必要な高レベルの部品群を実装したライブラリ系列で、Index Setクラス、ベクトルクラス、行列クラスの他に、クリロフサブ空間法およびプレコンディショナ、非線形ソルバー、時間発展演算子といった異なるクラスから構成されます。ベクトルおよび行列クラスは最低レベルの抽象化であり、様々な機能のコアブロックです。

PETScのOpenMP版ではベクトル、行列クラスがスレッド化されています。クリロフサブ空間法とプレコンディショナはKSPクラスとPCクラスに実装されています。これらはベクトル/行列クラスを用いるため明示的なスレッド化はされていません。SORや不完全LU分解などのプレコンディショナは、その複雑なデータ依存性のためスレッド化されていません。スレッド化には並列効率を改善するようなアルゴリズムでの再設計が必要です[11]。

lock exchange問題をテストケースとして、BoomerAMGプレコンディショナとFluidity自身の持つマルチグリッドプレコンディショナの性能を圧力ポアソンソルバーで比較しました。結果としてBoomerAMGはマルチグリッドに比べ良好な性能を示しました。ただし長さスケールに大きな多様性が存在する場合は、Fluidityのマルチスケール法のみが選択肢です。

最終ベンチマークテスト

lock exchange問題をテストケースとしました。これは古典的なCFDテスト問題で、最初にlockはタンク内の異なる密度の2つの流体を分離しており、lockを取り去ると2つの重力による流れがタンクに沿って伝搬するという問題です。3.4M頂点メッシュ数と82M流速自由度数を用いました。ベンチマークはHECToR Interlargosの2ノードから1024ノードまで行いました。
行列アセンブリは32Kコアまでスケールしました。256ノードでは、純粋なMPI版高速化率が99.3倍に対して、ハイブリッド版の高速化率は107.1倍でした。これはローカルアセンブリにより、コードが本質的にローカルプロセス内に閉じているためです。

プレコンディショナにHYPER BoomerAMG、ソルバーに共役勾配法を用いた圧力ポアソンソルバーの性能は、純粋なMPI版が64ノードから性能劣化し、ハイブリッド版がその性能を上回ることを示しています。

ここではファイル出力をせずにメッシュファイルなどの入力のみを計測しました。結果として大きなI/O効率の向上が示されました。例えばノードあたり4プロセスのみを用いると、ノード上のメタデータ操作量が削減されます。逆にこうしないと性能に悪影響が出ます。さらに、全ハロメッシュサイズはパーティション数（つまりプロセス数）と共に増加します。よって、ハイブリッドOpenMP/MPI手法は、計算ノード上の要求メモリー量、ディスクへの書き込み総量、メタデータ操作量を削減します。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学$CSE$サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。

文献

[1]	Xiaohu Guo, G. Gorman, M Ashworth, A. Sunderland, Developing hybrid OpenMP/MPI parallelism for Fluidity-ICOM - next generation geophysical fluid modelling technology, Cray User Group 2012: Greengineering the Future $CUG2012$, Stuttgart, Germany, 29th April-3rd May 2012
[2]	Xiaohu Guo, G. Gorman, M Ashworth, S. Kramer, M. Piggott, A. Sunderland, High performance computing driven software development for next-generation modelling of the Worlds oceans, Cray User Group 2010: Simulation Comes of Age $CUG2010$, Edingburgh, UK, 24th-27th May 2010
[3]	Welsh, D. J. A.; Powell, M. B., An upper bound for the chromatic number of a graph and its application to timetabling problems, The Computer Journal, 10$1$:8586, 1967 doi:10.1093/comjnl/10.1.85
[4]	P. Berger, P. Brouaye, J.C. Syre, A mesh coloring method for efficient MIMD processing in finite element problems, in: Proceedings of the International Conference on Parallel Processing, ICPP'82, August 24-27, 1982, Bellaire, Michigan, USA, IEEE Computer Society, 1982, pp. 41-46.
[5]	T.J.R. Hughes, R.M. Ferencz, J.O. Hallquist, Large-scale vectorized implicit calculations in solid mechanics on a Cray X-MP/48 utilizing EBE preconditioned conjugate gradients, Comput. Methods Appl. Mech. Engrg. 61$2$, $1987a$, 215-248.
[6]	C. Farhat, L. Crivelli, A general approach to nonlinear finite-element computations on shared-memory multiprocessors, Comput. Methods Appl. Mech. Engry. 72$2$, $1989$, 153-171.
[7]	D. Komatitsch, D. Michaa, G. Erlebacher, Porting a high-order finite-element earthquake modeling application to NVIDIA graphics cards using CUDA, J. Parallel Distrib. Comput. 69, $2009$, 451-460.
[8]	C. Cecka, A.J. Lew, and E. Darve, Assembly of Finite Element Methods on Graphics Processors, Int. J. Numer. Meth. Engng 2000, 1-6.
[9]	A.H. Baker, M. Schulz and U. M. Yang, On the Performance of an Algebraic Multigrid Solver on Multicore Clusters, in VECPAR 2010, J.M.L.M. Palma et al., eds., vol. 6449 of Lecture Notes in Computer Science, Springer-Verlag $2011$, pp. 102-115
[10]	Fluidity Manual. Applied Modelling and Computation Group, Department of Earth Science and Engineering, South Kensington Campus, Imperial College London, London, SW7 2AZ, UK, version 4.1 edn. $May 2012$, available at https://launchpadlibrarian.net/99636503/fluidity-manual-4.1.9.pdf
[11]	M. Weiland, L. Mitchell, G. Gorman, S. Kramer, M. Parsons, and J. Southern, Mixed-mode implementation of PETSc for scalable linear algebra on multi-core processors, In Proceedings of CoRR. 2012.
[12]	Michael Lange, Gerard Gorman, Michele Weiland, Lawrence Mitchell, James Southern, Achieving efficient strong scaling with PETSc using Hybrid MPI/OpenMP optimisation, submitted to ISC'13, 2013