＊ここに掲載するのは、エジンバラ大学、EPCCのI. Bethune博士によるHECToRレポート「CP2K - Sparse Linear Algebra on 1000s of cores A dCSE Project, I. Bethune, EPCC, The University of Edinburgh」を要約したものです。
概要
CP2K[1]は、Fortran 95で書かれたオープンソースの原子分子シミュレーションプログラムです。このコードは、Quickstep [2]リニアスケーリングDFTアルゴリズムの実装で最もよく知られていますが、古典的、半経験的、DFT、Hartree -Fock、および最近追加されたM?ller-Plesset二次摂動理論（MP2）の扱いが可能で、拡張性と効率性の高い方法で設計されています。CP2Kは、800,000行以上のコードで構成され、スイスに本拠を置く分散型チームによって開発されています。
CP2Kは、CASTEP（9.21％）とVASP（12.11％）に次いで、システム使用率8.48％（731193 AU）を消費しています（HECToR、HEC 2010年12月21日?2011年10月3日のデータ）。したがって我々は、dCSE資金提供[3] [4]の下で行われた作業は、全体としてのHECToRの科学的成果に大きな影響を与えると考えています。このレポートの3.3節では、HECToRユーザーへのCP2Kに可能な限り最高のパフォーマンスを達成させる推奨方法を示します。
·HECToRについて
英国国立スーパーコンピューティングサービスHECToR [7]（High End Computing Terascale Resource）は、Cray XE6ハードウェアで構成されています。 このプロジェクトでは、システムのフェーズ2bとフェーズ3を使用しました。
フェーズ2b XE6システムは、1856の計算ノードで構成され、各ノードは2.1 GHz、12コア「Magny-cours」AMD Opteronプロセッサと32 GBメインメモリを搭載し、コアあたり1.33 GBを提供します。各プロセッサはCray Geminiルーティングチップで接続され、高帯域幅と低遅延の3Dトーラスネットワークが装備されます。システムのピーク性能は360 TFです。
フェーズ3は、16コアの2.3GHz「Interlagos」AMD Opteronプロセッサに置き換えられ、1ノードあたり合計32コア、1コア当たりメモリ1GBです。 2012年1月に別の10個のキャビネット（928個の計算ノード）が追加され、ピーク性能が820以上に向上します。
·Rosaシステム
Hutter教授とHP2Cプロジェクトのサポートにより、スイスCSCSのRosa（XT5 12コア）、Palu（XE6 24コア）、Todi（XK6 16コア）システムにアクセスしてテストを行いました。
プロジェクトの目的と結果のまとめ
このプロジェクトとは独立にDBCSR[5]ライブラリが開発されました。これらを分離することは困難ですが、以下では全ての作業内容と成果を示します。
コード性能ベンチマークには、最新のHECToRハードウェア（フェーズ3、XE6、ノードあたり32コア）を用い、CP2K 2.1.390（2010年1月10日）とCP2K 2.3.r12105（2011年12月22日）の2つのバージョンのCP2Kを比較します 。ベンチマークデータのH2O-64およびH2O-1024は、それぞれ64および1024水分子の立方体セルの分子動力学実験用データです。これらは現在CP2Kで研究されている典型的でかなり大きなシステムを表しています。 両方ともGTH Basis Set [9]（TZV2P）とPBE [10]、およびPADE交換相関関数を使用します。
OpenMP/MPIハイブリッドモード：計算集約的なpure-MPI実行と同等のパフォーマンスを持ち、通信集約的な実行ではpure-MPI実行より約30％の改善を予想しました。
32コアでのH2O-64ベンチマークは、ほぼ計算集約的です（ランタイムの約15％が通信）。ここでは、ハイブリッドコードがpure-MPIコードよりも低速でした（2スレッド9％、4スレッド42％）。しかし、ランタイムの約40％がDBCSRで費やされているため、コードの他の重要な領域にスレッド性能を上げるべき領域が存在します。
通信集約的な領域（例えば、512コア、39％は通信）では、MPIプロセスあたり2つのスレッドを使用すると、pure-MPIコードに比べて22％高速化します。
より大きなH2O-1024ベンチマークは、実行時間の約70％をDBCSRに費やします。ここでは、ハイブリッドモードとpure-MPIの間で同等のパフォーマンスを達成するという目標は、計算集約的な実行（256コア、12％は通信）で達成されました。 特に、2010年10月版のコードと比較すると、1プロセスあたり2スレッドで性能が18％向上し、4スレッドで35％向上しました。
H2O-64と同様に、コアの総数が増加して通信が支配し始めると、ハイブリッドモード・コードの性能が向上し始めます（4096コアにおいてpure-MPIに比べて29％向上）。
CSCSのJoost VandeVondele博士によりCray XE6 `Monte Rosa 'システムで上において、2?2048個の水分子に同様のベンチマークが行われました。プロセス当たり2OpenMPスレッドを使用した場合、すべての問題サイズに対して最速の最短時間を達成できました。
計算集約的な領域に対しては10%の改善を予測しましたが、上述のように大幅な性能向上を示しました。
通信集約的な領域に対しても10%の改善を予測しましたが、H2O-64データについては2スレッド、512コアで4%の向上がみられましたが、他のケースでは劣化しました。4096コアを用いたH2O-1024データでは、pure-MPI2010年10月版よりも6%性能向上しましたが、プロセス当たり2スレッドでは4%劣化しました。
DBCSRライブラリの最適化
DBCSR (Distributed Block Compressed Sparse Row)は、CP2Kに組み込まれており、スパースブロック構造行列の保存形式と積演算に用いられます。CP2Kで用いられる大きな行列（密度行列、重なり積分行列、コーン・シャム行列など）は空間内のガウシアン基底の局所性により、自然にスパースな性質を持ちます。ブロック構造は各原子が複数の基底関数で表現されるために生じます。よってN原子系の典型的な行列は、N行(列)のブロックを持ち、各ブロックは複数の行と列で構成されます。例えば、最小基底関数(SZV-GTH-MOLOPT)の水素原子は1、より大きな基底関数(DZV-GTH-MOLOPT)の酸素原子は13です。つまり全行列はサイズが異なる小さな密行列を行と列として構成されます。ブロックはCSR保存形式でアドレスされ、個々のブロックは先頭行へのポインターとオフセットでアドレスされます。行列は2Dグリッドを用いてMPI並列で分散されます。
行列積C=A*BはCannonアルゴリズムで実行します。簡単に言えば、行列はP個のプロセスの正方グリッド上に分散され、Pの平方根のステップ数で実行されます。各ステップではローカルデータエリアが掛け算されて結果の行列に加算されます。そして、ノンブロッキングMPIを用いて、AとBそれぞれに対して行方向と列方向にシフトされます。
OpenMPの利点を生かすために、DBCSRは行列のブロックの行を利用するスレッド数で分割します。行をスレッドへ割り当てることで負荷バランスを行い、スレッド当たりの行数はほぼ均一になります（スパース行列なので、各行はほぼ同数の非ゼロブロックを持ちます。）。その後、これら行は、スレッドの持つ行がメモリー上で連続になるように並べ替えられます。
·ノード上のローカルデータ行列の積
行列積演算では最初に、DBCSRは行列の行のループを実行し、各ブロックでは相当する行列Bの列の全てのブロックのループを実行して行列Cへ加算します。各ブロックの行列積はBLASのDGEMMを用います。また、既存の乗算方法はキャッシュ効率が低いため、再帰的乗算スキームを実装しました。
·小型行列積用ライブラリ
特殊な小行列乗算ライブラリ `libsmm 'を組み込みました。 CP2Kではトリックスブロックのサイズは、シミュレーションされる原子種と基底関数に依存します。 通常、これは1x4,5x13,9x22などの特殊なサイズの小さなブロックになります。libsmmは、指定された小さな行列サイズのセットに対して特殊なDGEMMルーチンを作成するためのオートチューニングアプローチを採用しています。
特に小さなブロックサイズ（または1つ以上のディメンションが小さいブロック）では、libsmmがCrayのlibsciのBLASよりも最大10倍優れていることが示されました。
·スレッド最適化
DBCSR行列は行ごとに分解され、各行は特定のOpenMPスレッドによって所有されてます。 既存の負荷分散戦略（行には各行のブロックサイズで重み付けられます）では、各行のスパース性を考慮しないため、負荷の不均衡が発生します。
負荷バランスを改善する方法を調べると、スレッド0は他のスレッドよりも最大20％長く（完全に負荷分散されている人工的な入力データを用いた場合でも）一貫して長くかかることが発見されました。 コードを慎重に調べると、！$ ompのマスタディレクティブを含むすべてのスレッドによって呼び出されるタイミング測定ルーチンが原因であることが判明したため、行列積再帰部分から削除されました。
スレッドに行を割り当てるための基準として、分散された行のブロックの総数を使用して負荷バランスの向上を行いました。 これにより、行列乗算のFLOPの全体的な負荷バランスが良好になります。 しかし、Cannonアルゴリズム乗算の特性から実際には大きな負荷の不均衡が観測されています。
·MPIデータ交換
Craypatプロファイリングから分かることは、乗算ループ中のMPIで費やされた時間の大半はMPI_Waitallであり、非常に大きな負荷分散（一部のプロセスは最速より6倍長い）がありますが、プロセス間で送信されるデータ量は各プロセスによって実行される計算量と同様に、十分にバランスされています。
そこで、MPI_Testany呼び出しを行うことによって、乗算中にマスタースレッドが定期的にMPIをポーリングするスキームを実装しました。
ここで使用したベンチマークは、新しい線形スケーリングDFT [6]の実装コードを使用した6144原子の液体の水の計算です。これは密度行列のSCFを実行し、DBCSRの性能が支配します。 2スレッドでは最大13％、の改善が得られることが6スレッドでは最大20%改善しました。
·行列積におけるデータ形式
プロファイル分析により冗長なコードを削除し、さらにその他の小さなOpenMP最適化の結果、128 MPI x 2 OMPの場合は8％、32 MPI x 8 OMPで43％高速化しました。
謝辞
Ben Slater博士 (University College London, CP2K User), Joost VandeVondele教授 (ETH Zurich, CP2K Developer), Jurg Hutter教授 (University of Zurich, CP2K Developer)の多くのサポートに感謝いたします。
このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | CP2K website - Open Source Molecular Dynamics, http://www.cp2k.org |
| [2] | QUICKSTEP: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach, J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassasing and J. Hutter, Comp. Phys. Comm. 167 (2005) 103-128 |
| [3] | Improving the performance of CP2K on HECToR: A dCSE Project, Iain Bethune, 2009, http://www.hector.ac.uk/cse/distributedcse/reports/cp2k/cp2k_final_report.pdf |
| [4] | Improving the scalability of CP2K on multi-core systems: A dCSE Project, Iain Bethune, 2010, http://www.hector.ac.uk/cse/distributedcse/reports/cp2k02/cp2k02_final_report.pdf |
| [5] | U. Borstnik, J. VandeVondele and J. Hutter, In preparation, 2012. |
| [6] | Linear scaling self-consistent field calculations with millions of atoms in the condensed phase. J. VandeVondele, U. Borstnik and J. Hutter, In preparation, 2012 |
| [7] | HECToR: UK National Supercomputing Service, http://www.hector.ac.uk |
| [8] | Cannon, L. E. Ph.D. thesis, 1969, AAI7010025 |
| [9] | Separable dual-space Gaussian pseudopotentials. Goedecker, S; Teter, M; Hutter, J. PHYSICAL REVIEW B, 54 (3), 1703-1710 (1996). http://dx.doi.org/10.1103/PhysRevB.54.1703 |
| [10] | Generalized gradient approximation made simple. Perdew, JP; Burke, K; Ernzerhof, M. PHYSICAL REVIEW LETTERS, 77 (18), 3865-3868 (1996). http://dx.doi.org/10.1103/PhysRevLett.77.3865 |
| [11] | Efficient Implemention of Sorting on Multi-Core SIMD CPU Architecture, Chhugani et al, Intel Research, http://pcl.intel-research.net/publications/sorting_vldb08.pdf |
| [12] | A Benchmark Parallel Sort for Shared Memory Multiprocessors. R. Francis and I. Mathieson, IEEE Transactions on Computers, 37:1619-1626, 1988. |
| [13] | Optimising the DBCSR GPU implementation, J. Chetty, 2011, http://www.epcc.ed.ac.uk/wp-content/uploads/2011/11/JayChetty.pdf |
| [14] | Million Atom KS-DFT with CP2K, I. Bethune, A. Carter, X. Guo and P. Korosoglou, 2011, http://prace-portal.cscs.ch/uploads/tx_pracetmo/CP2K.pdf |
| [15] | Mixed-mode MPI/OpenMP for the CP2K User Community, I. Bethune, A. Carter, K. Stratford, P. Korosoglou and M-F. Iozzi, 2012, In preparation |