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

＊ここに掲載するのは、エジンバラ大学EPCCのMark Filipiak博士によるHECToRレポート「Mesh reordering in Fluidity using Hilbert space-filling curves, Mark Filipiak, EPCC, University of Edinburgh, March 2013」を要約したものです。

概要

Fluidityは、Applied Modelling and Computation Group $AMCG$ at Imperial Collegeで開発されたオープンソースの汎用マルチスケールCFDモデルです。Fluidityの開発目的は海洋モデリングです。Fluidityは、全球、海盆、地域、プロセススケールを同時に解析することが可能で、他のモデルでは不可能だった海洋研究を可能にします。またこれはシミュレーション中に計算自由度の数や配置を動的に最適化可能なアダプティブ法モデルです。長期的な目標は大循環シミュレーションおよび、特に海洋境界流、対流プルーム、北西ヨーロッパ棚の潮汐前線の1km解像度までの結合ダイナミクスの解析シミュレーションを行うことです。

現在Fluidityは、デスクトップワークステーションからHECToRのようなHPCプラットフォームまで様々な計算機において実行できます。Fluidityは現状では、解析対象の問題サイズは特にHECToRのようなNUMA$Non Uniform Memory Access$マルチコアプロセッサシステムにおいて制限が生じます。本プロジェクトの目的は、有限要素アセンブリとソルバーにおいてOpenMP/MPIハイブリッド並列の性能をより引き出すメモリーアクセス手法を実装し、Fluidityのスケーラビリティを改善することです。

プロジェクトの目的はFluidityのトポロジカルメッシュと計算メッシュの並び方を改善することです。これは、Fluidityの分割ソフトウェアfldecompに新たな機能を追加し、Fluidityのコアにメッシュ再順序化コードを組み込むことで行い、NUMAの性能をベンチマークしました。

最終的にヒルベルト空間充填曲線法による再順序付けが実装され、トポロジカルメッシュに適用されました。不連続ガラーキン要素やより高次の要素が実行中に生成されますが、これらの再順序付けの品質はトポロジカルメッシュの並びの品質に強く依存します。よって再順序付けの直接的な利点を見るために、一次区間線形の連続ガラーキン離散化を用いた後方ステップ流れのテストケースを検証しました。

再順序付け手法

NUMAアーキテクチャでローカルメモリーにメモリーページが確実に確保されるように、Fluidityにはファーストタッチ原則が実装されています。スレッド・プロセス親和性はユーザ環境で指定されます。しかしながらこれは、構築ステップ中に要素を構成する物理ノードのデータがメモリ内に互いに近接して配置されていない場合、またはソルバーステップ中にノードの物理的近傍がメモリ内に互いに近接して配置されていない場合には、メモリーアクセスの局所性を得られません。

一般に、メッシュの構築手法はデータ局所性よりも正しいメッシュを構築するように設計されています。メッシュ生成においてメッシュを再順序付けしたとしても、続くMPI用の領域分割で局所性は劣化します。Fluidityのデータ構造でのメッシュ再順序付けがデータ局所性を改善するのは、要素アセンブリ中に各要素データがメモリー内で近接しており、ソルバーではノード間メモリー距離が小さく、行列のバンド幅が狭いためです。

ノード並びの改善には多くの手法が用いられています。例えば、計算メッシュではnested dissection、四分木/八分木法、空間充填曲線など、行列では、逆Cuthill-McKee法[1]等です。これら手法はトップダウン法かボトムアップ法のどちらかです。典型的なトップダウン法はnested dissection（入れ子分解）で、メッシュは指定された最小サイズまで再帰的に分割され、最終的な二分木は深さ優先順に生成されて最小サイズのメッシュ順序を得ます。各最小サイズメッシュ内のノードの任意の順序付けにより、すべてのノードの順序付けが得られます。典型的なボトムアップ法はヒルベルト空間充填曲線法で、メッシュで覆われた物理領域に対して（有限）ヒルベルト曲線が計算されます。この時ノードのヒルベルト座標はノードの順序を示します。

Dissectionや曲線による階層構造の生成は、順序付けにおいて望ましいノードの局所性を与えます。 ここで、順番的に近いノードはメッシュ内または物理領域内内でも近い場所にあります。しかしながら、メッシュまたは物理的な領域では近くても、順序でははるかに離れているスケールにおける境界ノードが存在し、これは避けることができません。どちらの方法も物理空間に階層構造を持っているため、すべてのスケールで局所性を持ちます。局所性特性は、物理的局所性がメモリの局所性に対応することを保証するために必要なものです。 階層特性は、実際のUMA領域数に関わらず一般的に適用できることを意味します。 階層構造はまた、キャッシュサイズまたはその数に関わらず、結果として得られる順序付けが良好なキャッシュ局所性を有することを意味します。

Nested dissectionはメッシュベースのため、海洋の薄い球形の殻のような、如何なるトポロジーにも適用できる利点を持ちます。
一方、再帰処理を停止するサブメッシュサイズの選択は、プロセッサのキャッシュサイズに依存しますが、通常はUMA領域のサイズよりもはるかに小さいためNUMAとは無関係になります。さらにバイセクション法よりも遅いです。

一方、ヒルベルト空間充填曲線法は、マシン精度がヒルベルト曲線の高次近似を生成するのに十分な大きさであると仮定すると、すべてのスケールで局所性が達成されるのでその順序付けはキャッシュやNUMA構成とは無関係となることや、高速であるという利点があります。
一方、順序付けは物理座標に基づいているため、海洋の球形シェルといった、アスペクト比が極端に高い物理領域や高次元の空間に埋め込まれた物理領域では結果が悪くなる可能性があります。近年の自己回避歩行（self-avoiding walk）によるメッシュ生成の研究[3,4]は、メッシュベースの空間充填曲線と類似していますが、所望の階層構造を持つことは明らかではありません。

再順序付けライブラリ

検討されたライブラリは、Boost graph library、Leda、Scotch、Zoltan、Chaco、Jostle、ParMetisです。Fluidityには、ヒルベルト空間充填曲線法と、ParMetisとScotchを用いた再帰的dissection法を用いるライブラリZoltanがあらかじめ実装されています。文献[5]は、これが非構造計算の特に良い選択であることを示しています。

実装

Fluidityの中心となる計算は、運動量（速度）、圧力、輸送方程式の構築と求解です。速度、圧力およびトレーサ場は、基本メッシュと同じ要素上に定義されますが、異なる順序（例えば、2次要素）および異なる連続性を有することができます。これら場のメッシュは、基本メッシュから派生し、同じ要素および要素の接続性を持ちますが、ノード数とノード接続数が異なり、別々に順序付けする必要があります。インペリアルカレッジ（Imperial College、London）では、このような導出されたメッシュの順序付けスキームを開発していますが、開始する前にも基本メッシュに対する良い順序付けも必要です。

Fluidityのデータ構造は階層構造で、現在のシミュレーション状態を保持する状態変数を持ちます。 この状態変数は複数の場を保持しており、各場はそれが用いる離散化に関する情報、つまりトポロジカルメッシュに関する情報を持ちます。さらに領域分割は、押し出しメッシュや周期メッシュでも機能する必要があります。押し出しメッシュは海洋シミュレーションで一般的に使用されるもので、まず2次元サーフェスメッシュが作成され、次に決められた深度へメッシュは下方向に押し出されてz層メッシュまたはσ層メッシュのいずれかが作成されます。この場合は2Dメッシュだけが領域分割されますが、全3Dメッシュも正しく順序付けられる必要があります。周期メッシュには、1つまたは複数の領域エッジ上の要素間の結合が含まれるため、要素の順序付けには考慮が必要です。

トポロジカルメッシュの番号を変更する必要があるケースは3つあり、$i$初期メッシュ入力時、$ii$メッシュのアダプテーション後、および$iii$負荷バランス中に、基本メッシュの負荷がバランスされた後、かつ場がプロセッサ間で再分散される前です。この場合には、すべてのMPIプロセスに渡るユニバーサル番号から、1つのMPIプロセスにローカルなグローバル番号へのマッピングも置き換える必要があります。
各プロセス個別にメッシュは最順序付けられます。これはNUMAとキャッシュのパフォーマンスを改善するのには十分ですが、プロセス間で共有されるハロ領域の最順序付けはできないことを意味します。Fluidityは、ハロの順序を固定することで通信と計算をオーバーラップさせます。このコードはFluidityの開発ブランチとして利用可能です$https://code.launchpad.net/~fluidity-core/fluidity/hilbert$。

結果

50以上のテストを行いましたが、構築計算の性能は、順序付けしていないメッシュの場合に比べて95%から105%となりました。流れの計算時間は、順序付けしていないメッシュの場合に比べて90%から120%となりました。これらは実際の研究に用いられるデータよりも小さいため、すでに良い順序となっています。
隣接ノードの配列インデックスの差（メモリー空間上の距離に相当）の計測結果から、最順序化したほうがノード間メモリー距離が小さくなっていることが観測されました。

·大規模ベンチマーク

Fluidityを用いた代表的な大規模データは、後方ステップ流れのDNSです[5]。これは全場に基本計算メッシュを用い（P1-P1シミュレーション）、32-128コアを用いて計算しました。メッシュは20時間ステップ毎にアダプテーションされます。

·MPI

問題の明確化のため計算途中のI/Oは止めて実行しました。最順序化の性能向上はMPIプロセス数に依存しました。2プロセスでの性能向上は、経過時間の全体で5%、圧力計算で10%、流速計算で7%でした。この計算では流速の構築とメッシュアダプテーションにほとんどの時間がかかっています。プロセス数を増やすと全体性能は向上せず、16プロセスでは5%劣化しました。

·OpenMP

ハイブリッド性能の確認のため、データサイズを大きくして、HECToRの1ノード上で1MPIプロセスを用いてOpenMP並列実行を行いました。ここではハロは存在しません（大規模データではその割合は非常に小さく、その状況のアナログとなります）。
シミュレーションサイズを約3万ノードから1Mノード（メモリーは約4Gを使用）まで増加させて、シミュレーションには20時間ステップと1アダプテーションのみを実行しました。
HECToRノードは32コアが搭載され、4つのUMA領域に分散されています。HyperTransportリンクはUMA1と2および0と3の間が最もバンド幅が狭いため、スレッドの配置はNUMAの効果を生かすように選択しました。

1スレッドの場合は全体として20%性能向上が見られましたが、2スレッドの性能は劣化しています。４、8スレッドではほぼ横ばいです。これは、再順序化による性能向上よりも、コードの非スレッドセクションのUMA間の移動に時間がかかるためです。

謝辞

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

文献

[1]	dCSE website, “Distributed CSE Support.” Available at : http://www.hector.ac.uk/cse/distributedcse/, accessed 10 March 2010.
[2]	T. Sloan and M. Mewissen, “SPRINTing with HECToR dCSE Proposal.”
[3]	G. Grimes, S. Moodie, J. Beattie, M. Craigon, P. Dickinson, T. Forster, A. Livingston, M. Mewissen, K. Robertson, A. Ross, G. Sing, and P. Ghazal, “GPXMacrophage Expression Atlas: A database for expression profiles of macrophages challenged with a variety of pro-inflammatory, anti-inflammatory, benign and pathogen insults,” BMC Genomics, vol.6, p.178, 2005.
[4]	A. Brazma, H. Parkinson, U. Sarkans, M. Shojatalab, J. Vilo, N. Abeygunawardena, E. Holloway, M. Kapushesky, P. Kemmeren, G. G. Lara, A. Oezcimen, P. Rocca-Serra, and S.-A. Sansone, “ArrayExpress? a public repository for microarray gene expression data at the EBI,” Nucl. Acids Res., vol.31, pp.68-71, 2003.
[5]	The R package, “The R Project for Statistical Computing.” Available at : http://www.r-project.org/, accessed 10 March 2010.
[6]	BBSRC web site, “Biotechnology and Biological Sciences Research Council.” Available at : http://www.bbsrc.ac.uk/.
[7]	EPCC website, “Edinburgh Parallel Computing Centre.” Available at : http://www.epcc.ed.ac.uk/, accessed 10 March 2010.
[8]	DPM website, “Division of Pathway Medicine.” Available at : http://www.ed.ac.uk/schools-departments/pathway-medicine, accessed 10 March 2010.
[9]	SPRINT’s website, “SPRINT: A new parallel framework for R.” Available at : http://www.r-sprint.org/, accessed 10 March 2010.
[10]	Jon Hill, Matthew Hambley, Thorsten Forster, Muriel Mewissen, Terence M Sloan, Florian Scharinger, Arthur Trew and Peter Ghazal, “SPRINT: A new parallel framework for R,” BMC Bioinformatics, December 2008.
[11]	The MPI parallel library, “The Message Passing Interface $MPI$ standard.” Available at : http://www.mcs.anl.gov/research/projects/mpi/, accessed 10 March 2010.
[12]	HECToR’s user support, “HECToR : Hardware overview.” Available at : http://www.hector.ac.uk/service/hardware/, accessed 10 March 2010.
[13]	J. Hill, “Parallel statistics with R,” EPCC News, vol.62, p.7, 2008
[14]	M. Dublin, “New software tool lets researchers user R across a computer cluster.” Available at : http://www.genomeweb.com/new-software-tool-lets-researchers-use-r-across-compute-cluster, accessed 29 March 2010, 2009.
[15]	Savvas Petrou, “D1.1 Work Package 1 deliverable. Performance analysis of parallel correlation function pcor.”
[16]	Savvas Petrou, “D2.1 Work Package 2 deliverable. Performance analysis of parallel permutation testing function pmaxT.”
[17]	NWS, “NetWorkSpaces for R.” Available at : http://nws-r.sourceforge.net/, accessed 3 March 2010.
[18]	Rmpi, “Rmpi for R.” Available at : http://www.stats.uwo.ca/faculty/yu/Rmpi/, accessed 3 March 2010.
[19]	R/parallel, “R/parallel framework.” Available at : http://www.rparallel.org/, accessed 3 March 2010.
[20]	papply, “papply - R parallel apply function.” Available at : http://math.acadiau.ca/ACMMaC/software/papply.html, accessed 3 March 2010.
[21]	SNOW, “Simple Network of Workstations for R.” Available at : http://www.stat.uiowa.edu/~luke/R/cluster/cluster.html, accessed 3 March 2010.
[22]	Wellcome Trust web site, “The Wellcome Trust.” Available at : http://www.wellcome.ac.uk/.
[23]	ECDF, “The Edinburgh Compute and Data Facility.” Available at : http://www.ecdf.ed.ac.uk/, accessed 10 March 2010.
[24]	The FF R package, “ff: memory-efficient storage of large data on disk and fast access functions.” Available at : http://cran.r-project.org/web/packages/ff/index.html, accessed 10 March 2010.
[25]	Y. Ge, S. Dudoit and T. P. Speed, “Resampling-based multiple testing for microarray data hypothesis,” TEST, vol.12, pp.1-77, June 2003.
[26]	The multtest R package, “multtest: Resampling-based multiple hypothesis testing.” Available at : http://cran.r-project.org/web/packages/multtest/index.html, accessed 10 March 2010.
[27]	Westfall,P.H. and Young, S.S., “Resampling-Based Multiple Testing: Examples and Methods for P-Value Adjustment.,” p.?340, 1993. Wiley, New York.
[28]	R User Meeting, “SPRINT Project.” Available at : http://www.nesc.ac.uk/action/esi/contribution.cfm?Title=1044, accessed 30 March 2010.
[29]	Data-Intensive Research Workshop, “SPRINT Project.” Available at : http://www.nesc.ac.uk/esi/events/1047/, accessed 30 March 2010.
[30]	HPDC 2010 Workshop, “Emerging Computational Methods for the Life Sciences $ECML$.” Available at : http://salsahpc.indiana.edu/ECMLS2010/ Accessed 25 April 2010.
[31]	ACM International Symposium, “High Performance Distributed Computing $HPDC$.” Available at : http://hpdc2010.eecs.northwestern.edu/ Accessed 25 April 2010.
[32]	useR! 2010, “The R User Conference 2010.” Available at : http://user2010.org/ Accessed 25 April 2010.
[33]	CRAN, “The Comprehensive R Archive Network.” Available at : http://cran.r-project.org/, accessed 30 March 2010.