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

＊ここに掲載するのは、STFC デアズベリー研究所のAndrew Sunderl博士らによるHECToRレポートdCSE Project Report:Future Proof Parallelism for Electron-Atom Scattering Codes on the XT4, Andrew Sunderland,Cliff Noble, Martin Plummer, Department of Computational Science and Engineering, STFC Daresbury Laboratory」を要約したものです。

概要

電子-原子および電子-イオン散乱データは、多くの科学技術分野における重要な物理現象の解析に不可欠です。これには、水銀蒸着照明に代わる環境的に安全な代替物の開発、次世代ナノリソグラフィーツールのレーザー生成プラズマ源、大気プロセスの理解、核融合プラズマ中の不純物の診断、天体物理学データの定量的解釈などがあります。さらに、レーザ点火核融合EU（およびSTFC-CLF）HiPERプロジェクト（http://www.hiper-laser.org）や、関連する実験を伴う実験室天体物理学、不透明度測定、レーザ励起中空原子および強磁場中の原子などは、原子散乱データの詳細な計算を要求します。これらアプリケーションの重要性にもかかわらず、利用可能な、多くの複雑な原子およびイオンに対する正確な衝突（理論的または実験的）データが多くありません。これは、電子衝撃励起およびイオン化閾値付近の中間エネルギーでのイオン化の場合に特にそうです。

R行列法[1]は、電子-原子（および分子）散乱問題のシュレーディンガー方程式（およびその相対論的拡張）に対する非経験的変分解法です。配位空間は、標的原子（分子）またはイオンを含む球によって区画され、その外側ではターゲットの波動関数は無視でき、交換相互作用は無視されます。球の内部では、全電子配置間相互作用(CI)計算が、散乱エネルギーとは無関係に系の完全な固有解を生成します。内部領域解の固有値と表面上の振幅から、球の表面上でエネルギー依存性R行列が形成されます。球の外側では、内部領域のCI展開から得られる多重極ポテンシャルによる散乱電子に対して、沢山の連立微分方程式を解かなければなりません。このR行列は、内側および外側領域の解と一致します。高位の「中間エネルギー」状態の励起と電離は、球内の高位励起状態と正のエネルギー状態を表現する平方積分可能な「擬状態」を含めることで可能になります[1]。別の中間エネルギー処理（現在1電子原子の場合）は、球の半径を拡張し、2つの電子がより大きな長距離の擬状態基底に一致するまで、分割されたブロックで解を伝播させます[2]。この大きな半径では、計算を終了するには「標準的な」外側領域のR行列コードが必要です。

ソフトウェアの概要

·PRMATパッケージとPFARMコードの概要

この30年来に渡り、R行列法ベースのプログラム群が作成されて、多くの正確な散乱データが作成されてきました[1]。しかしながら、多くの重要な問題は、シリアル実行向けに設計されたプログラムでは実用的でないため、電子-原子散乱のための並列時間依存シュレディンガー方程式のFortran95コード「PRMAT」がEPSRCにより設計、実装されました[3]。PRMATは英国国立スーパーコンピューティングサービスHECToR [4]とHPCx [5]で提供されるアプリケーションパッケージの1つです。これは、シリアルコードRMATRX2 [6]に基づくRMATRX2/95と、シリアルコードFARM(Flexible Asymptotic R-matrix Package) [7]を用います。RMATRX2/95は内部領域の計算を行います。PFARMは、RMATRX2/95の結果を使用してエネルギー依存R行列を作成し、この行列を球から外向きに伝搬させることによってすべての散乱チャネルで連立微分方程式を解き、解を漸近的境界条件に一致させることで、個々の原子およびエネルギー積分後の温度依存形式の両方を生成します。複雑な原子では、高密度の散乱エネルギーに対して解が要求されます。PRMATパッケージは、Fe、Ni、Sn、中性Oの様々なイオンと電子衝突するデータを計算するために使用されています。軽原子による中間エネルギー散乱の研究にも使用されています[8]。

最近、PRMATパッケージは、PFARMにおける散乱チャネルの数を現在よりも遥かに大きくされて相対論的効果（例えば、開殻d-シェル原子やイオンの詳細な計算に必要とされる）を含むように拡張されました。これは、非常に多数のチャネルが球内の電子連続体の離散化から生じる、中間エネルギー散乱の場合にも当てはまります。 さらに、低エネルギー電子散乱における共鳴構造の複雑さは、正確な有効衝突強度を得るためには、典型的には数万の散乱エネルギー値で決定される衝突断面積を要求します。このプロジェクトの目的は、PFARMコードの並列パフォーマンスを向上させ、メモリをより効率的に使用した、将来的な（Fortran 2003機能も加味した）コードFARM2を開発することです。

PRMATパッケージは、電子原子（イオン）散乱のためのソフトウェアです。2009年初頭には、クイーンズ大学ベルファスト、UCL、オープン大学（Open University）、およびSTFC CSEの4つの共同機関に5年間のソフトウェア開発助成金が与えられました。理論的な電子-原子/分子散乱、ならびに時間依存と時間非依存のレーザー-原子（分子）相互作用の様々な分野の専門知識を集約するグラントです。2009年10月に開始したこのプロジェクト「UK-RAMP」の過程で、最新の機能を搭載した分子R行列コード（UKRmol-inおよびUKRmol-out）の統一セットが構築され、レーザーイオン化および他の相互作用のための2電子系時間依存直接解法コードHELIUMからその機能が導入されました[9]。UK-RAMPもまた、時間非依存R行列内部領域データを初期値として、HELIUMのテクニックをレーザー場内の一般原子に対処するように拡張します。UK-RAMPは、レーザーフリーの時間依存R行列内部領域コードは含みません。目的は、少なくとも原子コードレベルまで分子内部領域を構築、最適化して纏め上げて、同時に一般的な時間依存レーザー相互作用コードを構築することです。

·既存のPFARMコードの概要

PFARMは「外部」配位空間を複数の動径部分に分割して、基底関数展開により各領域においてグリーン関数を解きます。この方法はBaluja-Burke-Morgan (BBM) method [10]に基づいています。
実装においては、BBM法の変形版が外部領域の散乱を定義する連立2次微分方程式の解法として用いられています。連続した大きな動径方向でのR-マトリックスは、有限の動径方向のセクタ内で定義されたグリーン関数を使用して得られます。また、グリーン関数はshfted-Legendre基底を用いて得られます。
PFARMは、標準的な最適化されたMPI並列数値ライブラリを利用することが出来ます。この計算はEXDIGとEXASという2つの独立したステージで行われます。EXDIGでは、各セクションのハミルトニアンが生成され、一度、全てのエネルギーについて分散データ並列で対角化されます。ここではScaLAPACKライブラリのPDSYEVDルーチンを用います。セクターハミルトニアンは全てのプロセスを用いて順番に対角化されます。
次に、EXASの特定のタスクに割り当てられたプロセスのグループで全体が構築されます。生成プロセスグループRMPROD（R行列生成）は、連続する散乱エネルギーの内部領域境界での初期R行列を計算します。伝搬プロセスグループRMPROPは、内部領域境界から外部領域境界に伝播する際のR行列のストリームが通過する可変なパイプラインを形成します。

並列化による利得は次のとおりです。

• 大きなサブ領域行列の負荷の高い並列対角化は、散乱エネルギーとは無関係に計算されます。これは、数百または数千の散乱エネルギーを含む計算に特に適しています。
• サブ領域の伝播計算では、レベル3のシリアルBLAS行列乗算ルーチンを多く使用します。高度に最適化されたBLASライブラリ・ルーチンは、利用するマイクロプロセッサ・アーキテクチャ特性を引き出すように設計されているため、ピーク性能がほぼ達成されます。HPCx [11]の単一ノードでの効率調査からは、PFARMがピーク性能の35％以上で2.1Gflopsを示し、最も高速なアプリケーションコードの1つであることが示されています。
• ハイエンドのほとんどのコンピューティングプラットフォームでは、非同期通信によって通信コストが削減されます。
• この手法は、機能グループ全体で十分なバランスのとれた計算負荷を備えています。機能グループRMPRODとRMASYに割り当てられるプロセスの数は、i）完全に動作可能なパイプラインを維持するのに十分な頻度で初期R行列が生成されること、および ii）衝突強度計算は、ボトルネックを避けるのに十分な速度で処理されることが出来る数であるべきです。

·プロジェクトの概要

2008年、PRMATパッケージ開発者は、英国工学および物理科学研究評議会（EPSRC）[12]に代わり、nAG（Numerical Algorithms Group）LtdとのHECToR分散計算科学工学（dCSE）プロジェクトの契約を結びました。 dCSEプロジェクトの目的は、計算専門家が以下のことを可能にすることです：

• HECToRにコードを移植する。特に新しいコード、または既存のコードでこれまでサポートされていない機能を有効にする。
• HECToRでのコードパフォーマンスとスケーリング性能（理想的には数千コア）を向上させる。
• 長期的なメンテナンス性を改善するためにコードをリファクタリングする。
• o高性能コンピューティング分野におけるアルゴリズム改善を活用する。

PRMATプロジェクトの目標は、PFARMコードの最適化と並列化、Arnoldi log-derivative (ALD) propagator [13]を用いた新たなコードを組み込むことです。ALDコードは、メモリー効率、および全体の性能を大きく改善します。スタンドアローンの新たなコードを開発し、将来的な保守性のために、Fortran2003の新規機能の実験と調査を行います。

ハードウェアの概要

HECToR [4]は、エジンバラ大学に設置されたHPCxコンソーシアムによって運営されている英国の最新の国立スーパーコンピューティングサービスです。HECToR Phase 1 Cray XT4システム（2007 - 夏2009）は、1416の計算ブレードを持ち、それぞれに4つのデュアルコアプロセッサソケットが搭載された合計11,328コアのマシンです。プロセッサはAMD 2.8GHz Opteronで、各デュアルコアソケットは6 GBのメモリを共有し、合計33.2 TBを搭載します。システムの理論ピーク性能は59Tflopsです。HECToR Phase 2a（2009 - 夏2010）XT4は1416の計算ブレードで構成され、それぞれに4つのクアッドコアプロセッサソケットが搭載された合計22,656コアのマシンです。プロセッサーは、AMD 2.3 GHz Opteronで、各クアッドコアソケットは8 GBのメモリを共有し、全体で合計45.3 TBを提供します。システムの理論ピーク性能は208Tflopsで、2009年11月のトップ500リスト[20]で20位に位置しました。

比較のため、STFC Daresbury Laboratory [15]のComputational Science and Engineering Departmentに設置された、160 IBM eServer 575ノードで構成されている前回の英国National Capability ComputingサービスHPCx [5]の結果も示します。各eServerノードには16個の1.5 GHz POWER5プロセッサーが搭載されており、システムには合計2536個のプロセッサーが搭載されています。フレーム当たり合計32GByteのメインメモリは、フレームの16プロセッサ間で共有されています。HPCxシステムのフレームは、IBMのハイパフォーマンス・スイッチを介して接続されています。 最終的な構成での理論ピーク性能は15.4Tflopsです。

XT4上での既存のPFARMの性能

テストはjj-coupling（即ち相対論効果）を含むFeIIIケースで行いました。その伝搬の計算には1181チャネル（R 行列の次元に相当）、細密領域では10677、租領域では205の散乱エネルギーを用いました。PFARMのみPGIコンパイラの-fastオプションでコンパイルし、他はデフォルトです。

EXASの細密領域の伝搬計算は2048コアまで良好にスケールします。この計算コストは全体のコストを支配するため、全体性能も2048コアまでスケールします。コードは4096コアまで高速化はしますが、8129コアではもはや性能向上は示されなくなります。これは散乱エネルギー数により並列効率が制限されてしまうのが原因です。

PFARMの外部領域コードの最適化

初期プロファイルからEXDIGステージが全くスケールしないことが示されました。少ないコアでは全体への影響は小さなものですが、8129コアでは全体の24%を占めるまでになります。事前調査[2,16]から、大規模計算では並列対角化がボトルネックとなることが判明しており、もし最適化された数値計算ライブラリを用いたとしても並列スケール性が制限されることがありえます。

·並列対角化ルーチンの性能調査

ScaLAPACK[17]には、標準的な一般化密対称ハミルトニアン問題に適用可能な固有値問題ルーチンが多くありますが、HPCｘで行われた以前の調査では[16]、分割統治法ベースのPDSYEVDルーチンが、全ての固有値ペアを要求するPFARMの並列対角化に最も適していることが判っています。
ScaLAPACKには最近、Multiple Relatively Robust Representationベースの新たなルーチンPDSYEVRが実装されたため、比較を行いましたが、おおきなコア数ではこれまでのPDSYEVDが良好な性能を示しました。

·EXDIGの最適化

既存のEXDIGは各セクターを順次対角化しています。性能改善のために、サブプロセスグループで各対角化処理を並行して並列計算させるように修正しました。さらに各プロセスの結果を並列出力させるように実装しました。
最適化されたEXDIGコードを使用すると、XT4の8192コアで速度が2倍以上になりました。

·EXASの最適化

グループRMPRODとRMASYは、i）完全に動作可能なパイプラインを維持するのに十分な頻度で初期R行列が生成され、 ii）衝突強度計算がボトルネックを避けるのに十分な速度で処理されるように、パイプライン前後で十分なサイズが用意されることが重要です。
XT4でのプロファイル結果から、RMPRODプールのコア数は288コア、パイプら稲垂のRMASYでは7コアが理想的なバランスでした。この時最良のケースは最悪ケースに比べ28%高速でした。
4096や8192コアではこの配置はほとんど改善しませんでした。プロファイルからは大きなMPI_WAIT時間が示されました。

HECToRフェーズ2aでのPFARMベンチマーク結果

·FeIII、jj-カップリング、部分波：L=3、S=P=0

コア数は1024から16384まで変化させました。新しい―度は既存コードに比べて2.14倍から2.89倍高速化しました。各ステージで適切なコア数を分散させた場合は、既存コードに比べて3.13倍高速化しました。

ALD実装

このプロジェクトではさらに、ALMプロパゲータをPRMATに組み込んだスタンドアローンプログラムを作成しました。ALDは、AlexanderとManolopoulos [13]によって導入された連立Schr?dinger方程式を解くための安定した伝播法です。ALDの行列のサイズは、より多くのセクタを用いますが、BBMアプローチと比較して大幅に低減され、より大きな問題でより良い性能を維持することが期待されます。
Fortran2003標準に従って、非同期I/Oを含めて実装が行われました。

Alexander [13]によって開発されたALD法は、並列対角化ルーチンの効率によって課されるボトルネックを回避する補完的な方法です。基本的な考え方は、各セクタ内のポテンシャルが対角化されて、線形ポテンシャル（linear reference potential）で近似できると仮定することです。LightとWalker（LW）[20]によって提案されたR行列項で定式化されており、固定されたリファレンスポテンシャルを使用します。ALD方式では、リファレンスポテンシャルは、2つの半径方向点における実際のポテンシャルを各セクタと一致させるように選択することができます。これらの2つの半径方向点が2次シフトされたLagendre多項式の根に対応するように選択された場合、この方法は、2次-3次多項式によって近似されたポテンシャルについて正確に一致します[13]。LW法とALD法には2つの利点があります。対角化する必要のある行列はチャネル数と同じ次数であり、各セクタ内の伝搬方程式は行列ではなくベクトルになります（対角行列の行列演算の形でキャストできます）。各計算中に保持しなければならないデータの量は最小限であり、僅か20~30,000の散乱チャネルに対応する行列は少数のプロセッサで対角化することが可能です。セクタ方程式の解は、LW法における三角関数とALD法におけるAiry関数です。 ManolopolousとAlexander [13]は、三角関数よりもはるかに労力が少なく、Airy関数と微係数を計算できることを示しています。したがって、ALDはLWよりも効率的です。小さいセクタがBBMアプローチよりも要求されるますが、大規模並列HPCマシンを使用してより多くのチャネルを処理することができます。さらに、ALDの「潜在的な可能性」として、共鳴領域が共鳴が存在しない領域とは異なるように扱う必要はありません。

文献

[1]	P G Burke, C J Noble and V M Burke, Adv. At. Mol. Opt. Phys. 54 (2007) 237-318.
[2]	NS Scott, MP Scott, PG Burke, T Stitt, V Faro-Maza, C Denis and A Maniopoulou, Computer Physics Communications (CPC) 180 (2009) 2424-2449.
[3]	A G Sunderland, C J Noble, V M Burke and P G Burke, CPC 145 (2002), 311-340.
[4]	HECToR - UK National Supercomputing Service, http://www.hector.ac.uk.
[5]	HPCx - The UK's World-Class Service for World-Class Research, http://www.hpcx.ac.uk
[6]	P G Burke, V M Burke and K M Dunseath, J Phys. B 21 (1994), 5341-5373.
[7]	V M Burke and C J Noble, CPC 85 (1995), 471-500; V M Burke, C J Noble, V Faro-Maza, A Maniopoulou and N S Scott, CPC 180 (2009), 2450-2451.
[8]	BM McLaughlin, MP Scott, AG Sunderland, CJ Noble, VM Burke, CA Ramsbottom, R H G Reid, A Hibbert, K L Bell, P G Burke, Atomic Data Nucl. Data Tables 93, 55-104 (2007); M Lysaght, PhD thesis, University College Dublin, Ireland (2006); MP Scott, CA Ramsbottom, CJ Noble, VM Burke, PG Burke, J Phys B: At. Mol. Opt. Phys. 39, 387-400 (2006); BM McLaughlin, A Hibbert, MP Scott, CJ Noble, VM Burke, PG Burke, J Phys B: At. Mol. Opt. Phys. 38, 2029-2045 (2005); M Plummer, CJ Noble, M Le Dourneuf J Phys B: At. Mol. Opt. Phys. 37, 2979-2996 (2004); see also 'Mathematical and computational methods in R-matrix theory', eds. M Plummer, J D Gorfinkiel and J Tennyson, CCP2, STFC Daresbury Laboratory, UK, 2007 (on line version http://www.ccp2.ac.uk).
[9]	JS Parker, BJS Doherty, KT Taylor, KD Schultz, CI Blaga and LF DiMauro, Phys Rev Lett 96 (2006) 133001.
[10]	K L Baluja, P G Burke and L A Morgan, CPC 27 (1982), 299-307.
[11]	Single Node Performance Analysis of Applications on HPCx, M. Bull, HPCx Technical Report HPCxTR0703 (2007), http://www.hpcx.ac.uk/research/hpc/technical_reports/HPCxTR0703.pdf.
[12]	HECToR - Distributed CSE Support, http://www.hector.ac.uk/cse/distributedcse
[13]	M H Alexander, J Chem Phys 81 (1984) 4510-4516; MH Alexander and D E Manolopoulos, J Chem Phys 86 (1987) 2044-2050.
[14]	TOP 500 Supercomputing Sites,http://www.top500.org.
[15]	STFC's Computational Science and Engineering Department, http://www.cse.scitech.ac.uk/.
[16]	Using Scalable Eigensolvers on HPCx: A Case Study, Ian Bush, Andrew Sunderland, Gavin Pringle, HPCx Technical Report HPCxTR0510 2005, http://www.hpcx.ac.uk/research/hpc/technical_reports/HPCxTR0705.pdf.
[17]	The ScaLAPACK Project, http://www.netlib.org/scalapack/index.html.
[18]	A Parallel Divide and Conquer Algorithm for the Symmetric Eigenvalue problem on distributed memory architectures, F.Tisseur and Jack Dongarra, SIAM J. SCI. COMPUT, Vol.20, No. 6, pp. 2223-2236 (1999).
[19]	PDSYEVR. ScaLAPACK's parallel MRRR algorithm for the symmetric eigenvalue problem}, D.Antonelli, C.Vomel, Lapack working note 168, (2005), http://www.netlib.org/lapack/lawnspdf/lawn168.pdf.
[20]	J C Light and R B Walker, J Chem. Phys. 65 (1976) 4272-4282
[21]	M Gailitis, J Phys. B 9 (1976); C J Noble and R K Nesbet, CPC 33 (1984) 399-411
[22]	N S Scott,F. J?z?quel, C. Denis and J.-M. Chesneaux, CPC 176 (2007) 507-521; http://www.lip6.fr/cadna.
[23]	see, for, example:http://en.wikipedia.org/wiki/Shared_memory
[24]	A G Sunderland, M Ashworth, C J Noble and M Plummer, Proceedings of the Cray User Group (CUG 2009), Atlanta, Georgia, USA, http://www.cug.org/
[25]	http://ccpforge.cse.rl.ac.uk