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

目次

1   イントロダクション

2    シリアル実行コードの改善

3    並列化の方法

4    並列性能

5    新しい通信ライブラリーの実装

6    結び

1   イントロダクション

本プロジェクトは、非圧縮性ナビエストークス方程式ソルバーのIncompact3Dの追加の開発として、その姉妹コード：圧縮性流体シミュレーション向けCompact3DおよびQuasiCompact3Dを扱います。これらのアプリケーションは、インペリアル・カレッジ・ロンドンの乱流、混流、流体制御(Turbulence,Mixing and Flow Control)グループおよび共同研究者らにより、最先端の乱流研究に用いられています。

本プロジェクトは、Incompact3Dに関する2番目のプロジェクトとして行われます。最初のプロジェクトでは、Incompact3Dのスケーラビリティは、2DECOMP&FFTライブラリー開発を通して大きく改善しました。2DECOMP&FFTは並列分散FFTインターフェイスを持つ2Dペンシル領域分割の通信フレームワークです[3]。その結果、Incompact3Dは現在、HECToR上の実研究計算において、4,000-16,000コアを用いて通常利用されています。

本プロジェクトは2つの作業項目から成ります。最初に、現状シリアル実行版である2つのコードを、2DECOMP&FFTフレームワークで並列化します。これにより大規模圧縮性流体シミュレーションが可能になります。次に、新しい通信ルーチンを導入して、空間的陰的スキームフレームワークに陽的スキームを埋め込みます。これにより粒子-乱流相互作用などのより幅広いレンジの研究が可能になります。本プロジェクトは6人月予算、作業期間は1年とし、2011年3月に開始されました。

2    シリアル実行コードの改善

Incompact3Dの姉妹コードは長い歴史を持ち、前回のプロジェクトで現代的なスタイルに書き直されたIncompact3Dとは違い、これらのコードはFortran77形式のため、並列化作業前に以下の項目について書き直されました：

1. 全てのソースを固定形式から自由形式へ変換しました。
2. 全てのcommonブロックをFrotran moduleへ置き換えました。
3. 主要なデータ構造全てを静的なものからallocatable配列へ変換しました。
4. 重要な変数を含むほとんどの配列を、重要な機能単位で共有される一つのモジュールで定義するようにしました。これにより、サブルーチン間で渡されるパラメーターの数が激減し、コードの可読性と保守性が向上しました。
5. nAGコンパイラーの整形機能により、プログラムのコーディングスタイルの整備とベストプラクティスのために、ソースコードの形式を整えました。
6. Incompact3Dにおいて新たに再作成されたサブルーチンを、これら姉妹コードへ適用しました。非圧縮性と圧縮性ナビエストークス方程式は基本的には別のものですが、空間微分の計算といった幾つかの基本的なアルゴリズムは同じです。こうした作業はこのようなルーチンの検証にも役立ちます。Compact3Dの微分計算においては、幾つかのバグが発見されました。これらは同等機能のIncompact3Dのルーチンとのクロスチェックで見つかりました。本プロジェクトの終了時には、これらコード間で約2500行のソースコードが共有されました。
7. 3D配列をスワップするような幾つかの旧式のコードがあり、これらは古いベクトルマシンでの性能最適化(ベクトル長を長くするなど)に用いられていましたが、将来のコード開発をシンプルにするために削除しました。
8. この他にも不要なコードが見つかりました：応力テンソルを保存するための配列(移流項と拡散項の計算時に用いられます)が、他の場所でのテンポラリーの作業領域として用いられていました。これは、メモリー負荷低減のためですが、コードの解読と保守を困難にさせます。並列実装においては、メモリーの制約は、通常は重要な要素ではありません。この問題に対処するために、コード開発者が新しい変数セットを導入することになりましたが、本プロジェクト期間内には実装されませんでした。

元のシリアル実行版のコードを、バージョン管理のためのSVNサーバーへ登録しました。これは、以前のIncompact3Dプロジェクトで用いた同じサーバーです。改修毎に標準テストケースの検証が行われました。全ての変更が完了した後に、新しいコードは元のコードと同じ結果を生成しました。
こうしたコードの最新化作業が、高品質の並列化作業には必要です。

3    並列化の方法

並列化は、以前のプロジェクトで作成された2DECOMP&FFTライブラリーを用いて行います。
2DECOMP&FFTライブラリーは柔軟に設計されているため、Compact3DおよびQuasiCompact3Dの並列化には幾つかのやり方が存在します。関係者へのコンサルティングにより以下の方法が採用されました：

1. 画面出力は、全てMPIランク0で処理する。
2. 全ての2D/3D変数は、2DECOMP&FFTで分散させる。
3. 1D配列はグローバル配列として置く。すなわち、1プロセスに対してローカルでない情報が使われない場合でも、各MPIプロセスはそれらのコピーを持つ。これは、メモリ使用量を多少増加させるが、コーディングを簡素化する。
4. 2D/3D分散データは常にグローバル座標を用いる。これは不便さを生じさせずに、様々な通信を簡素化する。よって、幾つかの2DECOMP&FFTの通信ルーチンやI/Oルーチンは、ローカルおよびグローバルデータ両方で定義されたデータを操作するように更新されました。
5. 分散された全ての配列は、最初はXペンシル形式で定義され、元のコードの変数と1対1に対応しています。これらは必要な場合にのみ転置されます。
6. 速度、密度、圧力、温度と言った主要な変数に対し、2D領域分割で定義された各ペンシル方向に対応して、3つのコピーが保存されます。これはコーディングを簡素化しますが、余分な通信オーバーヘッドを導きます。これは後で議論されます。
7. グローバルな最大/最小値計算やグローバルメッシュライン上の平均値計算などのグローバル操作に対しては、MPIのリダクション操作を用います。これらは他のアプリケーションでも同じく必要なため、2DECOMP&FFTライブラリーへの実装が次のリリースで計画されています。
8. 転置操作ベースの並列コード作成においては、元のシリアル実行版アルゴリズムは、通信数の削減のため、多くの場合書き直されます(例えば、計算順序の並べ替え、中間結果保存のためのテンポラリー作業領域の確保、計算の複製化など)。コード内のアルゴリズムを良く知ることから、コード開発者が上述の実装を行いました。よって、並列アルゴリズムのあるセクションと元のコードの間で1対1に対応するように、並列化が行われました。これは、スケーラビリティーよりも絶対的な性能に影響を与えるでしょう。

この手続きは、将来のコード開発に役立つように考えられました。開発者が新たな機能を追加する際は、シリアル実行コードを作成するように新しいアルゴリズムを考えることが出来ます。上述のやり方に従えば、新しいアルゴリズムの並列化は容易になります。

MPIコードの検証は、通常は全ての作業を終えた後に行いますが、本プロジェクトでは、途中のコードでも単一のMPIランクで実行して試験するやり方を定式化しました。これはシリアル実行コードと同じ結果を出力するべきものです。このやり方はコーディングのバグ発見に役立ちました。

並列化作業を完了後、コードはデバッグのためクラスターへ持ち込まれました。デバッグは、分散配列内の中間データを、シリアル実行によるデータと比較することで行いました。検証は、セクション毎に全ての結果がマシン精度で一致するまで行われました。さらに、領域に関する積分値(エントロピーなど)の検証を行いました。これらは領域分割手法に関わらず一定でなければなりません。最終的に、並列化コードは、シリアル実行コードと数値的に一致する結果を示しました。

4    並列性能

新しいコードはHECToR上へ移植され、検証とベンチマークが行われました。HECToRでのベンチマークには、ユーザーが提供したテストケースを用いました。これは75×101×75の3Dメッシュ点によるシミュレーションです。このサイズは実際の用途におけるものよりも2ケタ小さいですが、時間の制約から大規模なベンチマークは実行していません。

図1：Compact3Dの強スケーリング性能

図1は、Compact3Dの強スケーリング性能を示しています。まず、シリアル版に対する性能向上は、64コアを用いて12倍であり、十分ではありません。しかしながらこの性能は、ベースとなる2DECOMPライブラリーの振る舞いに合わせて大きく振れます。ALLTOALL型の通信は、小さなケースでは大きなオーバーヘッドを生じますが、大きなケースでは良好にスケールすることが可能です。これは、小さなテストケースのメッシュ数を各次元で倍にした、人工的なベンチマークテスト結果から支持されました。32から128コアの間で並列効率がほぼ80%に達しました。

5    新しい通信ライブラリーの実装

Incompact3Dは、空間的に陰的スキームで境界から境界へのデータアクセスが要求されることから、転置ベースの手法で並列化されました。しかしながら、ステンシルベースの陽的スキームが望ましい状況は多くあります。例えば、

1. Large-eddyシミュレーション(LES)においては、小さなスケールの乱流は解かれた流れ場を用いてモデル化されます。サブグリッドスケールのモデル構築には、多くの場合空間微分が必要です。このような微分計算に陰的スキームを用いるのは高価で不必要です(通常はモデル化された量に対して数値精度を求めることはありません)。
2. 粒子-乱流相互作用のようなアプリケーションで用いられる、粒子追跡(Particle tracking)は、固定した領域分割と陽的スキームで容易に実装できます（注意すべきことは、粒子追跡に計算領域を分割するのは大規模問題に対しては適していません。粒子追跡は流れ場の計算よりも計算量が大きく、粒子の分布が一様でない場合は、領域分割は大きな負荷分散を生じます。もしそうであればタスクベースの分割を考えるべきです。）。

陽的スキームのサポートのために、第2の通信ルーチンセットが導入されました。これは、ペンシル(2D領域分割で定義される)周囲にハロ(halo)セルを構築し、1対1通信コードを作成して、隣接ペンシルでデータ交換を可能にするものです。

目標は、Xペンシルでのみ粒子追跡を行うことです。ここでは空間補間に4次の差分スキームを利用します(5点ステンシルを用いるので、通信コード構築時に2層のハロセルが必要です)。ここで、この陽的スキームの利便性のために、ハロセル数に関わらずサポート可能なように、より一般的な通信ルーチンを設計し、2DECOMP&FFTへ実装しました。

関連する仕事として、P3DFFTパッケージ[2](これも2Dペンシル領域分割を使用しています)は、以下のような保存の仕方でハロセルをサポートしています：

|000000000000000000000000000|11111|22222|33333|44444|55555|66666|
              ^                ^     ^     ^     ^     ^     ^
         Pencil data           W     E     S     N     B     T

メモリーの最初のセクションは、ペンシル内のデータを含みます。ハロセルのデータはこれに追加される形で、west/east(第1次元)、south/north(第2次元)、bottom/top(第3次元)と続きます。この保存形式の利点は効率性です。同じメモリー参照が、グローバル転置ルーチンで用いられています。そこではハロセルは単純に無視されます。不利な点は、ハロセルのデータを、それを利用する前により便利なijk順の形式へマッピングする追加のコードが必要となる事です。

現状ではよりユーザーフレンドリーな設計が施されています。2DECOMP&FFTへ、以下の形式のライブラリールーチン一つが導入されています：

call update_halo(var, var_halo, level)

ここで、最初の引数varは入力配列で、領域分割で定義された通常のペンシル分散データです。サブルーチン呼出し後、2番目の引数var_haloが出力で、元のデータにMPIで隣接プロセスから集められたハロデータが追加されています。3番目の引数levelは、オーバーラップするレイヤーが幾つ必要かを定義します。var_haloは、呼出しルーチンで、3Dのallocatable配列あるいはポインターとして定義されていなくてはなりません。ライブラリーは、そのメモリー空間を計算し、アロケーションして、サブルーチンから戻った時点では、呼出しルーチンで通常のijk順の形式で用いることが可能です。更に、以下の一般的な形式で用いる事も可能です：

call update_halo(var, var_halo, level, opt_decomp, opt_global)

2つのオプショナル・パラメーターopt_decompおよびopt_globalは、任意のグローバルデータサイズとグローバル座標を可能にするものです。
転置ベースのアルゴリズムにおいては、境界条件は、その次元がローカルメモリー内に全て存在するとして、シリアル実行コードと同様に適用できます。このライブラリーコードは、周期境界条件をサポートします。よって領域の境界となるペンシルは、反対側の境界のペンシルと情報を交換します。
ライブラリーは内部的には、標準的なMPIノンブロッキング・1対1通信ルーチンを用いて、隣接ペンシル間通信を行います。ハロセル・サポートの詳細は、2DECOMP&FFTのドキュメントを参照してください[1]。

新しい通信ルーチンは、現在Incompact3Dにおいて、数種の重要なアルゴリズム実装において用いられています。例えば、埋め込み境界法を用いて計算領域に固体をモデリングする場合、Incompact3Dは、流体-固体界面の正しい境界条件を規定するために3次元補間を実行する必要があります。こうした補間は、陽的スキームを用いて、新しい通信ルーチンを用いる事により可能になります。その粒子追跡アルゴリズムは、メッシュ点から任意の空間位置への速度場の補間を含み、ハロセル・サポートも要求します。

6    結び

本プロジェクトを通して、Incompact3Dの姉妹コードが並列化されました。第2の通信フレームワークがIncompact3Dの能力拡張のために導入されました。当初のプロジェクト目的は全て達成されました。

セクション3で議論したように、新しい並列コードの性能には、依然として改善余地があります。それは、アルゴリズムに通信ライブラリーを有効に活用させることです。これは、ソフトウェア工学よりもアルゴリズム設計に多く依存します。

前回のプロジェクトと同様に、再利用可能なソフトウェアコンポーネントが開発されて2DECOMP&FFTライブラリーに導入されました。2DECOMP&FFTの新しいバージョン1.4は、2011年10月にリリースされました。

謝辞

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

文献

[1]	2DECOMP&FFT website. http://www.2decomp.org.
[2]	P3DFFT website. http://code.google.com/p/p3dfft/.
[3]	N. Li and S. Laizet. 2DECOMP&FFT - a highly scalable 2D decomposition library and FFT interface. In Cray User Group 2010 conference, May 2010.