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

＊ここに掲載するのは、STFC デアズベリー研究所のC J Noble, A G Sunderland, M Plummer両博士によるHECToRレポートdCSE project report, Combined-Multicore Parallelism for the UK electron-atom scattering Inner Region R-matrix codes on HECToR, C J Noble, A G Sunderland and M Plummer, Department of Computational Science and Engineering, STFC Daresbury Laboratory」を要約したものです。

概要

電子-原子および電子-イオン散乱データは、多くの科学技術分野における重要な物理現象の解析に不可欠です。これには、水銀蒸着照明に代わる環境的に安全な代替物の開発、次世代ナノリソグラフィーツールのレーザー生成プラズマ源、大気プロセスの理解、核融合プラズマ中の不純物の診断、天体物理学データの定量的解釈などがあります。さらに、レーザ点火核融合EU（およびSTFC-CLF）HiPERプロジェクト（http://www.hiper-laser.org）や、関連する実験を伴う実験室天体物理学、不透明度測定、レーザ励起中空原子および強磁場中の原子などは、原子散乱データの詳細な計算を要求します。これらアプリケーションの重要性にもかかわらず、利用可能な、多くの複雑な原子およびイオンに対する正確な衝突（理論的または実験的）データが多くありません。また、現在の（グローバルな）関心事は、多光子イオン化、高強度超高速プロセス、電界支援再結合および散乱、高調波生成およびコヒーレント制御のような、レーザ場における電子 - 原子（イオン）相互作用です。英国では、このような原子・分子のレーザー場プロセスを扱う非経験的コードは、例えば、コラボレーティブコンピューティングプロジェクトCCPQ（以前のCCP2）[1]によってサポートされているUK-RAMPコンソーシアムによって開発されています。

PRMAT [2]は、電子 - 原子および電子 - イオン散乱のための多数電子シュレーディンガー方程式の変分解法[3]に基づく非経験的R行列法に基づくプログラムの集合です。相対論的拡張が開発されており、より正確な散乱データを生成することが可能になっています。このパッケージは、例えば、天体物理的用途（星間媒質、惑星大気など）のFeおよびNiおよび中性Oの様々なイオンの電子散乱データ計算や、プラズマモデリングや核融合炉の不純物（例えばSn、Coのイオン、進行中のWなど）などの他の用途にも使用されます。 R行列計算では、配位空間は原子または分子の「ターゲット」を中心とした球によって2つの領域に分割されます。球の内部では、全電子対配置間相互作用(CI)計算を実行して、各散乱対称性に対する有限体積内の完全な（エネルギー非依存性の）ハミルトニアンを構築し対角化して、境界上のエネルギー依存R行列の構築を準備します。超並列外側領域コードPFARMは、以前のdCSEプロジェクト[4]の主題でした。このプロジェクトは、以下に示すように、3つまたは4つのスタンドアロンモジュールを含む内部領域コード、特に非相対論的内部領域コードに焦点を当てています。このコードはHECToRにおいて、UK-RAMP、特に'Atoms for Astrophysics'プロジェクトで使用されています。

RADモジュールは、計算に関与する全ての動径束縛状態および連続軌道について、1電子および2電子の動径クーロンおよび交換ハミルトニアン積分と多重極積分を計算します。系全体がN + 1電子（N電子ターゲット）から成る場合、 N個の結合軌道（「擬似軌道」を含む）と1個の連続軌道の反対称積と、N + 1個の結合軌道の反対称積を含みます。この時の連続軌道は球教会でゼロではありません。ANGモジュールは、様々なオープンおよびクローズチャネルのCIターゲット状態と連続軌道の組み合わせの角運動量カップリングと「N + 1」積カップリング計算を行います。これは最も複雑なモジュールです。HAMモジュールは、N電子ターゲット状態とN + 1電子散乱ハミルトニアン行列とN電子多重極、 必要に応じて、N + 1電子電気双極子マトリックス要素を作るための、各非相対論的「LSπ」（軌道、スピン角運動量、パリティ）散乱対称性に必要な多様な動径積分および角運動量カップリングを収集します。その後ハミルトニアンは対角化され、固有値と固有ベクトルの境界上の振幅および電気双極子を出力します。これはPFARM（またはシリアル版FARM [5]）コードまたは相対論的（ 'Jπ'）計算の出力を修正および再合成するコードのいずれかに出力されるか、 またはUK-RAMPマルチフォトンコードに適用されます。

プロジェクトの開始時に、パッケージは2つまたは3つの有効なパッケージに分割されました。散乱計算のために、固定された球半径境界条件での入力束縛軌道と直交する動径連続関数を数値格子上に生成しました（ 従来の方法[2,3]）。RAD、ANG、HAMはシリアルコードでした（ANGにはいくつかのOpenMPディレクティブが含まれていました）。RADコードはH van der Hart（Queen's University Belfast）によりOpenMP並列化されていました。

ユーザコミュニティは、CCPForge [8]で管理されている単一のパッケージとしてパッケージをまとめて戻すことを望んでいたため、最適化と拡張はアプリケーション全体に適用されました。RADの並列処理を拡張し、HAMを並列化してより野心的な計算（複雑な原子と固定された拡張境界条件ベースが問題となった拡張の両方）を可能にすることも望まれていました。最も重要なのはANGコードの並列化で、これは球面テンソル代数を用いる非常に洗練された手順でしたが、複雑な計算では大きなボトルネックになっていました[9]。

プロジェクトの目標と結果

18人月で、HECToRの XE6 32マルチコアの複数ノードを用いて作業を行いました。これは、OpenMP/MPIハイブリッドコーディングとIPC共有メモリセグメント利用の両方を探索して、pure-MPI構造でマルチコアを活用することでした。内部領域コード間とよりポータブルな最終出力の並列I/Oも意味します。ANGコードは数値計算やHPC技術と同様に理論原子物理学（および他電子角運動量カップリング）を必要とする最も数学的に複雑なコードのため、C J Noble博士が担当し、Forran2003構造を用いて高水準コードと理論の関係を簡素化しながら、コードを高性能かつ将来的に保証します（以前のdCSEプロジェクト[4]への貢献のように）。A G Sunderland博士は、RAD開発を担当してHAMの開発を協力して行いました。全計画は12人月として計画されました。RADのノード内並列はOpenMPで実装され、ANGとHAMは 共有メモリセグメントを含むMPIを介して並列化されました。

·RADの目標

1.	すべてのアプリケーション用に単一のRADコードを生成するために、B-スプラインと「通常の」連続関数生成を単一のインタフェースで結合します。
2.	OpenMPマルチコアの並列性をバインドされた軌道に拡張し、一般的なスキームを「通常の」連続関数に拡張します。
3.	ｌ個の変数（生成と積分の組み合わせ）と、関連するメモリ分布を持つ個別の境界またはオープン（連続）セクションにわたって、より高いレベルのタスク管理並列処理とMPI通信を導入する。 'xstream'並列I/Oライブラリを使用して、並列データの追跡と最終データの書き込みを行います。

目標1と２は完了しました。OpenMPコードはより効率的になりました。これまでのすべてのテストケースでは、従来のRADコードの実行時間が大幅に短縮され、Bスプラインコードの軌道と格子点が増えましたが、積分計算は改善されました。
目標3は、外側ループで作業を並列化する修正を完了しました。

·ANGの目標

1.	スケーラビリティとパフォーマンスを向上させるためのOpenMPコードの変更、「threadprivate」とバインディング要件の簡素化を行い、'xstream'並列I/Oライブラリを組み込んで書き込まれたデータを追跡します。
2.	SMP並列処理を他のコードセクションに組み込みます。System Vモジュールを使用して並列処理を実装します。
3.	対称性と「オープン」および「クローズ」セクションの外側ループに、タスク管理MPIコミュニケータを導入します。

プロジェクト開始時のANGコードには、並列処が可能な2つの場所にOpenMP指示行が挿入されていました。IBM Power-series HPCxでは正常に動作していましたが、HECToRでは競合が発生していました。ここではMPI並列化の作業を優先することとしました。
新しいANGコードと2つのユーティリティパッケージshmとpfilehandによるMPI / shmアプローチは大きな成功を収めました。shmパッケージは、各共有メモリセグメントをポインタ経由でアクセスするFortran 2003オブジェクトとなるように開発されました。ANGでは、shmは、並列に計算された「表面係数」のテーブルを格納するために使用されます。

'pfilehand'は、非同期MPI-IO呼び出しを使用するダブルバッファ並列I/Oモジュールです。これは、 バルクデータを検索するための最小限の情報用に予約されたインデックスレコードを持ち、後続の並列読取りが、並列書込みを実行するために使用されるコミュニケータとは完全に独立していることを可能にします。これはFortran 2003オブジェクトとして再度アクセスされ、MPI並列で計算された最終データを書き込むためにメインのANGコードで使用されます。この並列処理は、さまざまなANGルーチンすべてに拡張されていますが、束縛状態間のカップリングは依然として最も時間が掛かります。pfilehandモジュールは、R行列コード全体で使用されるダイレクトアクセスファイル用の旧シリアルファイルハンドルルーチンから開発されましたが、一般用とにも使いやすいものです。これはANGでは、表面係数から必要なカップリングを計算するために、カップリングインデックス上の並列化された外部ループと組み合わせて利用されます。新しい並列ANGコードは、これまでのシリアルコードから最大256コアまでのスケーリングがテストされました。Shmとpfilehandモジュールはテスト実行のスタンドアロンパッケージとして利用できます。これらはCCPforge [8]で公開を目指しています。

·HAMの目標

1.	MPIサブタスクコミュニケータは、LSπ散乱対称性および行列要素の分類（開殻-開殻、開殻-閉殻、閉殻-閉殻、これはハミルトニアン要素とダイポール要素の両方に）をループする。ハミルトニアンの出力データファイルを同義に並べて、PDGによって並列にピックアップする。異なる対称性の組み合わせに対して双極子データを分配する。
2.	マトリックス要素の構築にOpenMP/System V並列性を導入する。
3.	（ハミルトニアン対角化後の固有ベクトルの係数から）ダイポール行列を再計算するためのルーチンと、フィールドアトムの他の固有関数データを適合させ、並列PDGコードのオプションとして組み込む。

行列要素作成ルーチンに並列ループが導入し、pfilehand がANGデータ（およびRADデータをオプションとして）をピックアップするために使用され、必要に応じてshmを実装しています。前のプロジェクトで開発された 'xstream'ルーチンが最終データのハンドリングに導入され、 '従来の'バイナリ出力またはPFARMおよびPDG優先XDR出力が自動的に選択されるように、最終データの処理に導入されました。以前のバイナリ出力またはPFARMおよびPDG優先XDR出力を自動的に選択することができます。HAMコードは、現在、散乱と多光子アプリケーションの両方に対する独立した並列コードです。

RADの詳細

OpenMP並列は2電子交換積分の後半と電子間potential積分に適用しました。ここで、不要な再計算を回避するためにループの入れ替えを行いました。また、個別にバッファへの書き込みをまとめて出力する様に再構成しました。さらに内部ループを展開して、その上位ループへOpenMP並列を適用しました。
角運動量毎に180個のB-スプライン関数を用いた酸素テストケースを実施しました。
ループの並べ替えは、実質的なシリアル最適化を意味し、少ないコアでは高速化しました。8以上のコアの場合、ファイル操作がスケール性を劣化させています。

新しく実装された並列レベル開発は、（i）実行開始時に全てのMPIタスクに入力構成を配布すること、（ii）一時ファイルをメモリベースのストレージに置き換えること、および(iii)基底関数と動径積分の両方の計算の角運動量lの値による外側ループをMPIタスク間へ分散すること、が含まれています。これにより、複数ノードの利用だけでなくノード内に最適なMPIタスク/OpenMPスレッドの組み合わせを決定することで、全体的なパフォーマンスを向上させることが可能です。
OpenMPの並列化を補完するMPIタスクの導入により、1ノード内の全体的なパフォーマンスが大幅に向上し、2ノードを使用して性能は妥当なスケール性を示しました。
データサイズを倍にすると、4ノードの実行が2ノードの実行よりも高速になりましたが、負荷分散の影響が出ました。タスクへの割り当てがラウンドロビン方式を取ることが原因です。
最終的に、負荷バランスのために2つのループを融合したところ、性能とスケーラビリティの向上が示されました。HECToRでの最適なスレッド数は8です。

ANGの詳細

R行列は、WignerとEisenbudにより原子核へのプロトン/中性子衝突を扱うために開発されました[3]。このアイデアは、この核衝突を、入射核が標的核の外部にありクーロン力が衝突粒子と入射核子との間の相互作用を媒介する場合と、内側で短距離力の核力で相互作用する場合に応じて、これらの核衝突を2つの別個の衝突問題として扱うことでした。ターゲット核の境界でこれらの2つの解の逆対数微分を一致させると、断面積、共鳴パラメータおよび他の物理的観測値を計算するために利用できる実行列とR行列が生成されます。
1970年代、P G Burkeは電子原子衝突を扱うためにこの形式を拡張しました。原子系では、入射電子がターゲット原子の電荷雲の内側にあるかイオン交換効果が重要である場合はPauli排他原理が必要となり、結果として衝突過程を記述する複雑な積分 - 微分方程式が生成されます。しかしながら、入射電子がターゲットの電荷雲の外側にある場合は、より簡単な非交換散乱問題を解くことになります。このR行列形式は、第一原理で衝突問題を扱うことを要する計算を単純化します。

クイーンズユニバーシティ・ベルファストで開発された最初の（そして広く使われている）パッケージRMATRIX-I [13]の後に、RMATRX-IIは、より複雑な電子配置へ対応するために、PG Burke、VM Burke（DL）およびKM Dunseath[9]によって原子系R行列理論の実装の効率を改善する目的で開発されました。より大きくより複雑な原子システムが考慮されるにつれて、初期コードの限界が明らかになりました。特に、ハミルトニアン行列要素を導出するのに必要な角度方向積分の計算は、より深くカップリングされたシェルの寄与を考慮せねばならず、計算上のボトルネックとなりました。 多くのシェルが衝突プロセスに関与するにつれて、配置の数は急激に増加します。Burke-Burke-Dunseath[9]は、表面衝突係数と呼ばれる一連の角運動量変換係数が導出される、原子衝突を処理するための改良されたアルゴリズムを発表しました。これは、ターゲット占有シェル内の一対の電子のカップリングツリーの頂部への再結合変換を表し、そこでは、相互作用に関与する1電子および2電子積分を介して相互作用することができる。これらの係数は、（電子の反対称性を考慮して）fpc係数、Racah係数および運動学的因子を含む。これらの係数の表が与えられると、角度積分は迅速かつ効率的に評価される。このプログラムは主にV M Burkeによって書かれました。
ANGでは、電子構成は2つのセットに分けられます：ユニークセットと同等セットです。後者の組のメンバーは、ユニークセットのメンバーと同じ角運動量特性を有し、それらが関与する動径波動関数のみが異なります。こうしてユニークなセットのみを詳細に扱う必要があります。ほとんどの内部ファイルは、filehandサブパッケージを使用して読み書きされます。固定レコード長の未フォーマットFortranファイルは、filehandインタフェース経由でダブルバッファリングを使用して処理されます。 プログラムの次元は、1994年に前処理によって決定されました。
その後RMATRXIIは、（1）標準Fortran 90、Fortran 95への変換と 動的メモリ割り当ての導入とモジュール構造への再編成、（2）ANGの主要部分に疎行列技術と派生データ型が導入され、（3）マイナーなアルゴリズム的およびプログラム的な改良が実施されました。また（4）ANGプログラムは、OpenMPを使用してIBMアーキテクチャ上で並列化されましたが、4?5スレッドを超えるとスケール性は急速に悪化していました。

ANGの計算には、与えられた電子配置の集合からのすべての可能なカップリングと、正規直交LSπカップリング基底の定義に必要な詳細な角運動量代数（または、JKπ結合基底：Kについては[3]を参照）が含まれます。RMATRXII / PRMAT内部領域形式、特にANG形式の数学的基礎は[9]に詳述されています。ANGは2段階で進行します。まず、必要な表面係数を特定し、次に適切な表面係数テーブルが完成するまで計算します。次に、テーブルを用いて、必要なすべての角運動量カップリング（直接および交換の1電子および2電子ハミルトニアンおよび多極マトリックス要素における連続-連続、連続-束縛、束縛-束縛）を計算します。新しいコードでは、両方のステージが並列化されています。表面係数計算が分散され、係数が共有メモリセグメントに書き込まれ、共有メモリセグメントがSMPノード間で更新されます。第2段階は、各ケース（c-c、b-c、b-b等）のカップリングインデックスを介して外側ループ上で並列化され、出力データは並列に書き込まれます。

RADテストで使用された酸素ケース、およびシリアルANGがHECToR上で実行できないFe +の大規模ケースについてテストをしました。両方の場合において、実行時間の大部分は、「束縛-束縛」角運動量カップリング計算ルーチン「angbb」が支配しました。
酸素テストは32コアまで非常に良くスケールし、MPIオーバーヘッドが顕著になるにつれスケーリングが低下し始めます。鉄ケースは64コアまで非常に良くスケールし、128コアで顕著なMPIオーバヘッドが示され、256コアが限界でした。

HAMの詳細

HAMモジュールは、一組の散乱対称性に対応するハミルトニアン行列要素を構築します。これらの行列要素は、PDGへの入力としてディスク上のファイルに書き込まれるか、直接対角化されます。その他、対象となるハミルトニアンの構築と対角化と、漸近ポテンシャル係数（入射電子とターゲットが離れている領域での散乱を規定するポテンシャルの多重極ポテンシャル係数）、多光子コード内の光イオン化断面積計算に用いる双極子係数を計算します。これらは全てCPU負荷が高いものではありません。しかしながらハミルトニアン要素を構築する作業はデータ集約的であるため、これに対処するためにpfhコードを組み込みました。
'ij'番目のハミルトニアン行列要素は、ParAngコードで計算された係数をRADで掛け算した1電子と2電子動径積分の線形結合です。 多くのR行列計算では、動径方向要素の数は多くなく、これらの要素はメモリに読み込まれて必要に応じて選択されます。特定の中間エネルギー計算ではこの数は十分に大きくなるため、より精巧な分割スキームが必要になります。新しいコードParHamは、pfhを用いて動径要素を読み取ることが可能です。角度係数および関連するインデックスはpfhのファイルから並列に読み取られます。ParHamは多くのテストで検証されました。

謝辞

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

文献

[1]	The UK-RAMP consortium is a 5-year EPSRC Software Development Call Grant shared between Queen's University Belfast, UCL, the Open University and STFC Daresbury Laboratory, commencing 1 October 2009 (EP/G055416/1, EP/G055475/1, EP/G055556/1, EP/G055599/1); CCPQ details at http://www.ccpq.ac.uk
[2]	A G Sunderland, C J Noble, V M Burke and P G Burke, Comput. Phys. Commun. 145 (2002) 311-340
[3]	P G Burke, 'R-matrix Theory of Atomic Collisions: Application to Atomic, Molecular and Optical Processes', Springer, 2011
[4]	'Future-proof parallelism for electron scattering codes,' A G Sunderland, C J Noble and M Plummer, http://www.hector.ac.uk/cse/distributedcse/reports/prmat
[5]	V M Burke and C J Noble, Comput. Phys. Commun. 85 (1995) 471-500; V M Burke, C J Noble, V Faro-Maza, A Maniopoulou and N S Scott, Comput. Phys. Commun. 180 (2009) 2450-2451
[6]	M A Lysaght, P G Burke and H W van der Hart, Phys. Rev. Lett. 101 (2008) 253001; Phys. Rev. A 79 (2009) 053411; A C Brown, S Hutchinson, M A Lysaght and H W van der Hart, Phys. Rev. Lett. 108 (2012) 063006
[7]	L R Moore, M A Lysaght, L A A Nikolopoulos, J S Parker, H W van der Hart and K T Taylor, J. Mod. Optics 58 (2011) 1132-1140; L R Moore, M A Lysaght, J S Parker, H W van der Hart and K T Taylor, Phys. Rev. A 84 (2011) 061404
[8]	http://ccpforge.cse.rl.ac.uk
[9]	P G Burke, V M Burke and K M Dunseath, J. Phys. B: At. Mol. Opt. Phys. 27 (1994) 5341
[10]	http://www.prace-project.eu
[11]	http://software.intel.com/en-us/articles/intel-parallel-inspector/
[12]	I J Bush, Technical Report HPCxTR0701 www.hpcx.ac.uk/research/hpc/technical_reports/HPCxTR0701.pdf , A Maniopoulou and C Armstrong, nAG HECToR Case Study CSE Report (2009) http://www.hector.ac.uk/cse/reports/castep_m.pdf
[13]	K A Berrington, W B Eissner and P H Norrington, Comput. Phys. Commun 92 (1995) 290-420 (this article subsumes earlier publications of the package)
[14]	http://www.pgroup.com/lit/articles/insider/v3n2a2.htm