＊ここに掲載するのは、エジンバラ大学・EPCCのAdrian Jackson博士によるHECToRレポート「Improved Data Distribution Routines for Gyrokinetic Plasma Simulations, Adrian Jackson, EPCC, University of Edinburgh, 27 January 2012」を要約したものです。
概要
GS2は磁気プラズマ内の低周波乱流の研究のために開発された初期値問題シミュレーションコードです。GS2は摂動分布関数に対するジャイロ運動論方程式を、プラズマ内の乱流電磁場に対するマクスウェル方程式と連成して解きます。これは実験において生成されたプラズマの微視的安定性の確認と、不安定性から生じた乱流の様々な性質を計算するのに用いられます。また、宇宙物理学や磁気圏などの自然界で生じるプラズマ中の乱流シミュレーションにも用いられます。
ジャイロ運動論シミュレーションはプラズマ内の荷電粒子種の分布関数の時間発展を解き、セルフコンシステントな電磁場内の荷電粒子の運動を考慮しています。この計算は5次元のデータ空間を用い、その内3次元は粒子の物理的位置、残りの2次元は粒子の速度です。完全な6次元の運動学的プラズマ計算（3次元速度）は、強磁場下の粒子の極めて高速なジャイロ運動のため多くの時間を要し、膨大な計算コストが掛ります。ジャイロ運動論シミュレーションは、この6次元データ空間を磁力線周りの高速なジャイロ運動を平均化して簡略化することにより、問題を6次元から5次元へ変換します（ここで、速度はエネルギーとピッチ角であり、ジャイロ位相角は平均化されます）。粒子の高速なジャイロ運動は完全には計算されませんが、ジャイロ平均として取り込むことにより低周波摂動を低コストで忠実にシミュレートすることが可能です。
GS2は様々な方面で利用されており、例えば線形な微視的安定性シミュレーションでは、成長レートはballooningあるいはflux-tube極限における陰的初期値問題アルゴリズムを用いて波数基底毎に計算されます。任意の波数における最速成長の（あるいは少なくとも減衰の）固有モードの線形および擬似線形性は独立に（すなわち素早く）計算することができます。完全に成長した乱流の非線形ジャイロ運動論的シミュレーションも扱うことが可能です。全てのプラズマ種（電子と様々なイオン種）は、平等にジャイロ運動論的な扱いが可能です。非線形シミュレーションでは、変動スペクトル、異常（乱流）加熱レート、粒子種毎の運動量、エネルギーの異常輸送係数を計算します。しかしながら非線形シミュレーションは計算コストが高く、現実的な時間で計算を行うには並列処理が必要になります。
GS2は完全に並列化され、Fortran90で記述されたオープンソースコードです。並列化はMPIで実装されており、GS2のシミュレーション計算は、分布関数の異なる部分を異なるプロセッサへ割り当てて分割されています。
与えられたGS2入力ファイル（実行すべきシミュレーションと領域分割レイアウトの指定）に対して、プログラムingenはシミュレーションに対する推奨プロセス数（あるいはスィートスポット）のリストと、ユーザが与えたシミュレーションパラメータに関するその他の有用な情報を提供します。こうした推奨値は、入力ファイルで与えられた情報から算出されます。これは可能な限り負荷バランスを良くする様にデータ領域を分割する（つまりプログラム性能を最適化して実行時間を最小化する）のに役立てるためです。推奨プロセス数リストは、シミュレーションの線形部分で用いられるデータレイアウトg_loに基づいています。ingenは非線形部分に適したプロセス数リストも提供します（xxf_lo、yxf_lo）。これはg_loにより推奨されるプロセス数と異なる場合があります。
GS2の性能最適化にはingenが提供する3つのプロセス数リストの内最も有利なものを選ぶ必要があります。しかしながらこれは何時も可能であるわけではなく、特にプロセス数が大きい場合は、ユーザは一般にg_loに最適なプロセス数を選んでいます。
このプロジェクトは、2010年にEPCCが行った前回のプロジェクト（GS2のFFTの更新）をフォローアップするものです。前回のプロジェクトの目的は、FFTW3を組み込み、GS2でのFFT利用を最適化することでした。その中で、FFTに関連するレイアウト間のデータ転送が実行時間の大きな割合をしめることが解りました。
特にルーチンc_redist_22およびc_redist_32がコスト高で、各々8％の実行時間（1024コアでのベンチマーク）を占めていました。これらのルーチンはFFTで必要になり、各々xxf_loからyxf_loへと、g_loからxxf_loへの再配分/再配置を行います。これは他のプロセスとの通信されるデータと受信するデータの処理と、ローカルデータの再配置を要求します。これらルーチンはバッファのパッキング時に間接指標アドレッシングを多用しています。しかしながら、この機能は追加の計算コストが掛り、特に現代的な計算ハードウェアにとっては、この配列変換はメモリーアクセスコストの増加を伴います。
再配分機能
最適化対象の機能は、FFTで必要なデータ変換を実行するルーチンc_redist_22とc_redist_32、およびその逆変換であるc_redist_22_invとc_redist_32_invです。各2DFFT計算では変換あるいは再配置が必要で、これらは2つのフェーズから成ります。最初にc_redist_32はy方向のFFTに必要な転置操作を実行し、中間的なxxf_loデータ構造からyxf_loレイアウトへ変換します。実空間での非線形項計算の後に、変換とFFTを_invルーチンを用いて逆に実行し、k空間での非線形項を得ます。
再配分ルーチンは、ローカルコピーとリモートコピーという2つの機能が有ります。ローカルコピーはデータレイアウト間の同じプロセスに属するデータを再配分し、リモートコピーは異なるプロセス間で通信すべきデータを再配分します（各プロセスへ送るバッファデータの収集、該当プロセスへの送信、および受信データのローカル配列へのアンパッキングが含まれます）。
HECToRフェーズ2bにおける代表的なベンチマークケースの測定により、以下のことが判りました：、
＊ルーチンのローカルコピー部分は、少数コアにおいてそのルーチン性能を支配しています。
＊プロセッサ/コア数が多い場合、例えば1024コアでは、リモートコピーの特にc_redist_32が性能に大きく影響します。
ローカルコピーの改善
結果として、4つのルーチンの内、c_redist_32以外の3ルーチンでは、新しいローカルコピーコードが既存のコードより約40-50％高速化しました。既存コード内の殆どの間接指標アドレッシングを除去しましたが、全てではありません。新しいコードの中で、ループの最初と最後のインデックスを得るために、間接指標アドレッシングがまだ残っています。間接指標アドレッシング（およびtoとfromのデータ構造）を完全に除去するには、新しいコードで必要なインデックスの最小限のセットの計算と保存を分離する必要があります（すなわち、構築したループの先頭でのtoとfromの値に対して）。これは、必要なインデックスの最小限のセットを保存のみする様に、init_x_redistとinit_y_redistを修正すれば可能です（現状はこの修正はしていません）。
更に具体的なメモリーアクセスパターン調査と改善により、c_redist_32は他のローカルコピーコードと同等の性能改善を達成しました。
リモートコピーの改善
既存コードでは、MPIブロッキング通信（MPI_SendおよびMPI_Recv）が、FFTと実空間での再配布におけるプロセス間のデータ転送に用いられています。MPIでのデッドロック回避のために、半分のプロセスを送信後受信の一方向を行い、残りの半分は受診後送信のみを行う様にしています。MPIブロッキング通信は他のプロセスからの受信待ちの不要なオーバーヘッドを生じます。そこでノンブロッキング通信を実装しましたが、大きな改善は見られませんでした。これは、既存の送受信パターン実行は直裁的で、送受信プロセスの組み合わせが極めて洗練されたもので、プロセスの”red-black”タイリング型の選択を行っていることと、ノンブロッキング手法の副作用として、全プロセスが同時に通信を行うためシステムノード上のネットワークに競合を生じさせていることが原因と考えられます。
もう一つの最適化対象領域は、プロセス間の送受信データのパッキングとアンパッキング処理です。リモートコピーコードは既存のローカルコピーコードと同様に間接指標アドレッシングを用いています。これを最適化したところ、ローカルコピーの場合と同様に実行時間が削減されました。
'アンバランス'なデータ分割最適化
先述の様に、GS2はg_lo分割においては”スイートスポット”コア数で用いられるのが一般的です。プログラムingenが与えるスイートスポットは、g_loデータ空間に対する最適なプロセス数ですが、g_loレイアウトとは異なるやり方でシミュレーションデータを分割する他のレイアウトxxf_loやyxf_loには通用しません。具体的に言えば、線形計算に都合が良い並列プロセス数を選ぶことは可能ですが、それは非線形計算を扱うための通信量を大きく増加させます。
各プロセスに対するデータ分割は、単純に全データ空間をプロセス数で割ることで計算されています。割り切れない場合は空のプロセスを生じるます。
また、ユーザが指定する全非線形シミュレーションのレイアウトもまたGS2の並列性能に影響を与えます。xxf_loからyxf_loへのデータ転送には、xおよびyデータ次元からのデータスワップが含まれます。もしxおよびyデータ次元が並列プログラムでプロセス間に分割されなければ、この転送は単純な各プロセスのメモリー内移動になります。しかしながら、もしxおよびyデータ次元がプロセス間で分割される場合、このスワップはプロセス間のデータ送信を含み、ローカルなデータ移動と比較して一般にはコストが増えます。シミュレーションで用いるプロセス数と選択されたレイアウトは、xおよびyデータ次元がプロセス間を跨ぐかどうかに関わらず影響を及ぼします。
そこで最終的に、xxf_loおよびyxf_loからのマッピングにおけるデータ通信の最適化に、アンバランス性を持つ分割機能を実装しました。この修正では、各プロセスに帰属するxxf_loおよびyxf_loのブロックサイズ計算部を置き換えました。均一なブロックサイズからプロセス数に対して異なるブロックサイズへ変更し、xあるいはyのローカルなデータ空間がプロセス数で明確に割り切れない場合に対応させました。
結果として、1536プロセスを用いたyxlesレイアウトの結果から、アンバランス分割を用いることにより、サブルーチンc_redist_22のリモートコピー機能コストが大幅に（ほぼ完全に）削減されました。このことは対応する逆ルーチンにも反映されており、c_redist_32およびその逆ルーチンのコストも大きく削減されました。アンバランス分割を用いる場合にこれはローカルコピーコストの上昇とバランスしており、アンバランス分割による最適化は、1536コア使用時のyxlesレイアウトでのGS2の実行時間を全体として15％削減しました。しかしながらxylesレイアウトを用いた2048プロセスの場合は、アンバランス分割は全実行時間に対して小さな改善しか示しませんでした。これはc_redist_22のリモートコピー実行時間に占める割合が小さいからです。
アンバランス分割機能は、GS2シミュレーションの最適なプロセス数範囲を広げるという利点も持ち合わせます。以前に述べた様に、ユーザは一般的にはingenプログラムが推奨するプロセス数を利用します。これはレイアウトに決められた順序でその変数の因子を個別に考慮しており、それらにマッチしたプロセス数リストを用います。アンバランス分割では、最適なプロセス数リストはより幅広くレイアウト変数の組み合わせの因子を用います。
まとめ
| コア数 | 512 | 1536 |
| 既存コード実行時間[秒] | 4435.20 | 2385.60 |
| 最適化コード実行時間[秒] | 4150.40 | 2074.20 |
| 高速化 | 7% | 17% |
これらの最適化により、以下の表に示す様に代表的なベンチマーク（10000繰り返し計算）において（512コアのxylesレイアウトのデータ分割xxf_loとyxf_loでは、xとyは共にローカルのためアンバランス分割は役に立ちませんが）、関節指標アドレッシングの最適化により7％の性能改善が示されました。1536コアのyxlesレイアウトではアンバランス分割と間接指標アドレッシング最適化は共に有効で、アンバランス分割が約5％の負荷分散を示すにも関わらず、シミュレーションの全実行時間の17％を削減しました。さらに2048コアでのxylesでは、20％弱の実行時間削減を確認しました。
謝辞
このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。