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

＊ここに掲載するのは、エジンバラ大学,EPCCのLawrence Mitchell博士によるHECToRレポート「Performance optimisations for CARP: a dCSE project, Lawrence Mitchell, EPCC, The University of Edinburgh, July 27, 2010」を要約したものです。

概要

不整脈研究パッケージ（CARP,The Cardiac Arrhythmia Research Package）は心臓の大規模シミュレーション研究で広く利用されているソフトウェアです。その目的は不整脈や不随意細動といった病状の処置の詳細な個別治療法を提供することです。in-silicoモデリングを臨床作業に組み込むこの領域は実用可能で、心臓再同期療法（不整脈）や臨床試験（副作用確認）および除細動器の開発などが含まれます。

シミュレーションには心臓のモデル、通常はMRIスキャン情報の離散化から得られたものが入力として必要です$Bishop et al., 2009$。CARPは様々なアーキテクチャをターゲットに開発されてきました。コンパイル時のフラグにより次の並列化の選択（同時も）が可能です：シリアル実行、OpenMP並列、GPU利用、MPI並列。このプロジェクトではMPI並列のみを対象とします。今回施す最適化は他のコンパイル時フラグにも適用されます。

主要な目的と背景

プロジェクト開始前の性能プロファイルから、大規模HPC環境でCARPを十分効率的に実行するために必要と判断された2つの領域が特定されました。その並列領域分割スキームは大きな負荷分散のため性能を劣化させ、出力ルーチンは大規模シミュレーションにおいて大きなボトルネックになっています。この作業の目的は、第1に領域分割の改良によりCARPの強スケール性を改善すること、第2に出力ルーチンを最適化することです

CARPは、全非構造3次元有限要素メッシュ上で心臓バイドメイン方程式$cardiac bidomain equations, Plonsey, 1988$を数値的に解きます。微分演算子はメッシュ上で離散化され、シミュレーションのカーネルループは、Ax=bの形式の大規模スパース線形システムを解きます。バイドメイン方程式には演算子分割法を適用し、3つの分離された方程式へ変形し、これらをleap-frog法で解いて解を求めます。この3つの方程式は、楕円型線形PDE、放物型線形PDE、および非線形ODEから構成されます。楕円型PDEの求解が計算時間を支配します$Linge et al., 2009$。

並列化

心臓は有限個の小セクションに分割されて入力メッシュとなります。各セクションは要素、要素の角は頂点として扱います。各頂点は固定された空間点を与えられ、心臓の形状を形成します。こうした要素は些か大きめですが、心臓を構成する細胞の近似と考えることが出来ます。シミュレーションの目的は、これら要素の各々上の電位ポテンシャルを解くことです。このポテンシャルは心筋の反発を引き起こし、心拍を生じます。
CARPが採用する並列手法は、この有限要素メッシュを領域分割することです。バイドメイン方程式はこれらの領域上で解かれ、領域間で通信するデータは必要に応じて異なるMPIプロセスで計算されます。線形システムの観点では、これは各プロセスローカルに行と列のサブセットを保存することに対応します。このブロック対角分割においては、行列要素は通信を必要とします。良好な並列負荷バランスのためには入力メッシュの領域分割性が決め手となります。

既存の分割方法

既存のCARPは頂点ベースの分割手法を採用しています。各プrセスは連続した頂点群をローカルとして捉え、要素もローカルかどうかを単純なタイブレークルールで決めます。この方法は頂点数をプロセス間に極めて平坦に割当てます。

残念ながら各プロセスがローカルと考える要素の分布は極めてアンバランスです。さらに、この分割手法はインターフェイスサイズの最小化を考慮せず、並列システムで大きな通信コストを生じます。

このアンバランスな要素配分は、行列のある部分は別の部分より密度が高いことを意味します。つまり分割内で大きなインターフェイスを持つ領域は、ブロック対角外の多くのエントリを含むことを意味します。この領域分割を用いる行列ベクトル積は大量の通信が必要です。こうした事から負荷分散とそれによる並列性能劣化を生じます。

新しいメッシュ分割

最適な負荷バランスのために、頂点ベースではなく要素ベースの分割を採用しました。目的はメッシュを同一サイズで、インターフェイスを最小限にして分割することです。これはグラフ分割問題です。そこでメッシュの双対グラフ$dual graph$を構築しました。ここで、メッシュ内の各要素をグラフ内のノードとして、要素間を共有する面をリンクとします。こうしてグラフの最小カットのk-wayバイセクション分割を行います。

最適なk-wayバイセクション分割はNP困難問題ですが$Garey and Johnson, 1979$、ここでは並列分割ライブラリParMetis $Karypis et al., 2003; Schloegel et al., 2002$を呼出すこととしました。

頂点の再番号付与

良好なメッシュ分割の後には、2つやるべきことが有ります。最初に、新しい分割に応じた要素データの再配分を行います。このため、各プロセスにおいて、他の全てのプロセスに配るべき要素を格納したバッファをパックします。このバッファはMPI_Alltoallvで交換され、その後バッファを新しい要素データ構造へアンパックします。

次に、新しい要素分割に応じた頂点を再度番号付けする必要があります。これはプロセス内ローカルな頂点番号を連続にするためです。プロセス数の増加に伴いインターフェイス要素が重要になります。この時表面/体積比は増加し、インターフェイスの頂点の所属を決めるやり方が重要になります。

ここでシリアル版と並列版の再番号付けアルゴリズムを実装しました。並列版は150百万個の要素ケースでは、プロセス数に依らず20秒で完了ました。

方程式を解く

前述の方程式を並列ライブラリPETSc$Balay et al., 1997, 2008, 2009$を用いて解きます。特に楕円型および放物型PDEについては前処理付き共役勾配法を用います。ODEは陽解法で解きます。放物型PDEもオプションで陽解法を使うことが出来ます。

楕円型PDEのプリコンディショナ

CARPでは、16プロセッサシステムでのhyper AMGプレコンディショナ$Falgout and Yang, 2002$のパラメータセットが用意されていました。そこで大規模計算において最適な性能を引き出すパラメータセットを再調査しました。PETScは多くのプレコンディショナの実行時オプションを持っています。更に他のマルチグリッド法、Trilinosパッケージやブロックヤコビ・プレコンディショナも調査しました。CARPにおいては、Trilinosマルチグリッド・プレコンディショナはhyper AMGと同等のスケール性を示しますが、速度は3倍遅いです。ブロックヤコビ・プレコンディショナはhyper AMGよりスケール性が高いですが、速度は20倍遅いです（大規模系ではこの差は1.5倍程度まで縮まります）。

hyper AMGが最速です。以下のオプションを設定しています。
-ksp_type cg
-pc_type hypre
-pc_hypre_type boomeramg
-pc_hypre_boomeramg_max_iter 1
-pc_hypre_boomeramg_coarsen_type HMIS
-pc_hypre_boomeramg_interp_type ext+i
-pc_hypre_boomeramg_P_max 4
-pc_hypre_boomeramg_strong_threshold 0

AMGプレコンディショナは、粗視化演算子が粗視化レベルで密行列を生じる可能性が大きいという問題を孕んでいます。この場合は疎行列よりも多くの通信が発生し、行列ベクトル積のスケーラビリティを阻害します。ここで選んだオプションは、粗視化レベルでも行列がスパース性を保っていることを仮定したものです。

放物型PDEのプリコンディショナ

放物型問題ではブロックヤコビ・プレコンディショナを用い、ブロック行列の近似逆行列計算には不完全コレスキー分解を用います。ここで、プレコンディション行列の計算中に通信を発生させないようにするために、ブロック対角行列の外に在る非零要素は全てゼロにしました。

マルチグリッド法を放物型PDEに適用して性能を調査しましたが、繰り返し数は少ないですが、個々の繰返しの計算コストは極めて高くなり、結果としてブロックヤコビを選択しました。

ODEソルバーの性能

CARPではODEの計算は完全にローカルに実行されます。これは各グリッドにおけるイオンレベルの計算に用いられます。次の時間ステップでの媒体の導電性に影響を与えます。結局、各プロセスが同数のODEを扱うならば、良好な負荷バランスを達成できます。更に、通信をさせないようにするために、PDEの場合と同じローカルグリッドをODEにも設定します。こうしてODEの頂点分割は自然にPDEと同じ物を利用可能になります。この方法で良好な負荷バランスとデータ局所性を得ました。

ODEのデータは分散されたベクトルに格納されます（イオンレベルと各頂点です）。既存コードはこの格納にPETScルーチンを用いています（VecSetValues、VecAssemblyBegin、VecAssemblyEndの一連の呼出し）。この一般性を確保するPETScルーチンによる全ベクトル組立には少なくとも一回のMPI_Allreduce呼出しが含まれます。ここでは分散配置されたベクトルへの格納は全てローカルであることが判っているため、このグローバル通信は不要です。よってこのローカルな値のセットをループで置き換え、16,384コアまでスケールすることを確認しました。

出力の改善

CARPの出力ルーチンはシリアル実行で記述されています。その殆どは、グリッドポイント上で定義された電位ポテンシャルが記述されたベクトルを出力ファイルに書き出す処理です。出力が要請された時点で、ランク0はそのローカル部分を書出し、続いて順次全ての他のプロセスからデータを受信し書出します。その際全てのプロセスは出力ルーチンが終了するまで待ち合わせをします。

この方法ではプロセス数が増えるに従って書き出し時間は増加します。ディスクへの書き出し時間は一定ですが、プロセス数が増えるとデータ収集に掛かる時間は増加します。このシリアル実行ボトルネックの除去のために、並列出力ルーチンを実装しました。

ポスト処理

問題はメッシュを再分割したため結果ベクトルと入力メッシュのグリッドの順序が異なることです。
出力データを基の順序に戻すためには、例えばPETScのVecScatterを利用して順序を入れ替える方法が有りますが、グローバル通信が発生してスケール性を劣化させます。そこで全てをMPI-IOで行う案を避けて、計算を阻害しない非同期出力を採用しました。

ノンブロッキング出力

計算を行うプロセスはディスク書込みをしないように出力ルーチンを書き換え、少数のプロセスが出力するようにしました。このプロセス数はユーザ指定可能です。出力ルーチンが呼び出されると、計算プロセスで定義された分散ベクトルは出力プロセスへスキャッターされます。この時点で計算は完了しており、シミュレーションは続行可能です。

出力プロセスでは、データ書込み前に基の順序へ変換されます。書込みにはシリアル実行とMPI-IO出力の選択を可能にしました。
この実装は真の非同期出力ではありません。特に出力ノードは全てのデータ受信をプリポスト処理しますが、メッセージは並列に処理されません。

結果として、大きな系で20倍（8129コア上の出力性能が、以前には約25MB/sだったが、yaku460MB/sへ）高速化しました。

ベンチマーク

以下の3種のデータ/式を用いました。
・TBunnyC1：ラビット心臓3Dメッシュ、862,515頂点、5,082,272要素、完全バイドメイン方程式
・OxfordHeart2：上記の詳細版、6,901,583頂点、40,992,163要素、完全バイドメイン方程式
・HumanHeart3：ヒト心臓3Dメッシュ、26,190,199頂点、152,891,134要素、モノドメイン方程式（放物型PDEおよびODEのみ）

TBunnyC1データにおいては、スケール性能は、既存コードは64プロセスまでしかスケールしなかったのに対し、新しいコードは256プロセスに達しました。実質的に2.5倍の高速化です。楕円型および放物型PDEソルバーのスケール領域は4倍に拡大しました。

OxfordHeartデータはTBunnyCメッシュよりも大きなデータです。TBunnyCと同様に楕円型および放物型ソルバはODEと比較して性能が良く有りませんが、メッシュ分割を選び、非同期出力ルーチンを加えれば大きなプロセス数では約4倍高速化し、良好な効率を示します。小さなプロセス数でも約3倍高速化しました。全シミュレーションは既存コードに比べ、約740％高速化しました。

HumanHeartデータは、我々の知る限りこのモデルは現状研究されている最大のモデルよりも2倍大きなものです$Niederer, 2010$。
CARPは楕円型ソルバーに2つの手法を用います。一つは共役勾配法を用いた陰的クリロフ法、もう一つは陽解法で行列ベクトル積が要求されます。陽解法は時間ステップ当たりの計算と通信量が少ないものです。
クリロフ法は内積計算のために集団通信が発生します。一方、解の安定性確保のためには、時間積分ステップを減らさなければなりません。各ステップは高コストでは有りませんが、通常は5倍ステップを実行しなくてはなりません。この時間ステップの追加にも関わらず、陽的な解法は陰的なものに比べ高速で、大きなコア数でスケール性が良いことが分かりました。
前の例ではスケール性は1024コアまででしたが、このケースは16394コアまでスケールします。実時間に換算すれば、既存コードは最良のケースで、1秒の実時間に対して72分のシミュレーション時間が必要です。新しいコードはこれが4分になり、1800％の加速です。数週間かかっていたものが数日で完了できます。

結論

プロジェクトの当初の目的は、並列分割と出力ルーチンの改良によりCARPの強スケール性を改善することでした。大規模系での性能は当初の目的を実現しました。ヒト心臓シミュレーションは5分以内で完了することが可能になりました。

この作業の結果はCARPの開発リポジトリへ組み込まれ、現在は世界のユーザが利用可能です。この成果はCUG2010で発表され、PNASへ論文投稿準備中です。その他vph20105でも発表を予定しています（本レポート作成時）。

謝辞

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

文献

[1]	M. Ainsworth. Personal communication, 2010.
[2]	S.Balay,W.D.Gropp,L.C.McInnes,andB.F.Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163-202. Birkhauser Press, 1997.
[3]	S. Balay, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.0.0, Argonne National Laboratory, 2008.
[4]	S. Balay, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang. PETSc Web page, 2009. http://www.mcs.anl.gov/petsc.
[5]	M.J.Bishop,G.Plank,R.A.Burton,J.E.Schneider,D.J.Gavaghan,V.Grau, and P. Kohl. Development of an anatomically detailed mri-derived rabbit ventricular model and assessment of its impact on simulations of electrophysiological function. Am J Physiol Heart Circ Physiol., 298$2$: H699-718, 2009.
[6]	R. Bordas, B. Carpentieri, G. Fotia, F. Maggio, R. Nobes, J. Pitt-Francis, and J. Southern. Simulation of cardiac electrophysiology on next-generation high-performance computers. Philosophical Transactions of the Royal Society A, 367:1951-1969, 2009. doi: 10.1098/rsta.2008.0298.
[7]	C. Chevalier and F. Pellegrini. PT-Scotch: a tool for efficient parallel graph ordering. ParallelComputing,34:318-331,2008. doi: 10.1016/j.parco.2007. 12.001.
[8]	M. Deo, P. Boyle, G. Plank, and E. Vigmond. Role of Purkinje system in cardiac arrhythmias. In Engineering in Medicine and Biology Society, 2008. EMBS 2008. 30th Annual International Conference of the IEEE, pages 149-152, 20-25 2008. doi: 10.1109/IEMBS.2008.4649112.
[9]	T. Edwards and K. Roy. Using I/O servers to improve application performance on Cray XT™ technology. In Cray Users Group Meeting $CUG$ 2010, Edinburgh, Scotland, 2010.
[10]	R.D.Falgout. An introduction to algebraic multigrid. Computingin Science and Engineering, 8:24-33, 2006.
[11]	R. D. Falgout and U. M. Yang. hypre: a library of high performance preconditioners. In P. M. A. Sloot, C. J. K. Tan., J. J. Dongarra, and A. G. Hoekstra, editors, Computational Science - ICCS 2002, volume 2331 of Lecture notes in computer science. Springer, 2002.
[12]	M. R. Garey and D. S. Johnson. Computersand Intractability: A Guideto the Theory of NP-Completeness. W. H. Freeman, 1979.
[13]	G. Karypis, K. Schloegel, and V. Kumar. ParMETIS: parallel graph partitioning and sparse matrix reordering library. University of Minnesota, version 3.1 edition, 2003. URL http://glaros.dtc.umn.edu/gkhome/ metis/parmetis/overview.
[14]	S. Linge, J. Sundnes, M. Hanslien, G. T. Lines, and A. Tveito. Numerical solution of the bidomain equations. Philosophical Transactions of the Royal Society A, 367:1931-1950, 2009. doi: 10.1098/rsta.2008.0306.
[15]	S. Niederer. Personal communication, 2010.
[16]	R. Plonsey. Bioelectric sources arising in excitable fibers $Alza lecture$. Annals of Biomedical Engineering, 16:519-546, 1988.
[17]	K. Schloegel, G. Karypis, and V. Kumar. Parallel statis and dynamic multiconstraint graph partitioning. Concurrency and computation: practice and experience, 14:219-240, 2002. doi: 10.1002/cpe.605.
[18]	K. St¨uben. Algebraic multigrid $AMG$: an introduction with applications. Technical Report 70, GMD - Forschungszentrum Informationstechnik GmbH, 1999.