＊ここに掲載するのは、エジンバラ大学、EPCCのI. Bethune博士によるHECToRレポート「Improving the scalability of CP2K on multi-core systems A dCSE Project, I. Bethune, EPCC, The University of Edinburgh, September 22, 2010」を要約したものです。
概要
このプロジェクトは、CP2KへOpenMP並列を実装するHECToRのdCSEプロジェクトです。前回のプロジェクト[1]から引き続き、エジンバラ大学EPCCのアプリケーションコンサルタントIain Bethune博士により、2009年9月から2010年9月まで6人月工数で行われました。UCLのBen Slater博士、およびチューリッヒ大学のCP2K開発者Joost Vandevondele博士からの支援に対し感謝いたします。
·CP2Kについて
CP2K[2]は、フリーの原子分子シミュレーションコードです。分子からバルク物質まで幅広い物質を対象に古典的なポテンシャルを含めた手法と、密度汎関数法(DFT)およびハートレーフォック法がベースです。さらにDFTにはQuickstep線形スケール法[3]を搭載しています。詳細は前回レポートを参照してください[1]。
·HECToRについて
HECToR(High End Computinging Terascale Resource)[4]は、Cray XT,X2で構成された英国国立スーパーコンピューティングサービスです。このプロジェクトでは、フェーズ2a(XT4)とフェーズ2b(XT6)を用いました。XT4は、5664計算ノード、各ノードは2.3 GHz quad-core`Barcelona' AMD Opteronプロセッサーと8 GBメインメモリーから構成されます。ノードはSeaStar2+ルータチップによるCray高速ネットワークで接続され、全システムはピーク性能208Tfopsです。
XT6は、1856計算ノード、各ノード2つのは2.1 GHz 12-core`Magny-cours' AMD Opteronプロセッサと32 GBメインメモリーで構成されています。ノードはSeaStar2+ルータチップによるCray高速ネットワークで接続されていますが、2010年後半にはGeminiネットワークに置き換えられます。ピーク性能は360Tflopsです。
·Rosaシステムについて
`Monte Rosa'[6]は、スイスCSCSに設置されたCray XT5 です。1844計算ノード、各ノード2つのは2.4 GHz hexa-core`Istanbul' AMD Opteronプロセッサと16 GBメインメモリーで構成されています。このシステムは、フェーズ2b設置前のプロジェクト半ばで利用しました。
·共同作業について
この作業と並行して、チューリッヒ大学Urban Borstnik博士が開発した新しいスパース行列ライブラリDBCSR(Distributed Block Compressed Sparse Row)がCP2Kに導入されました。このライブラリは当初から、スケーラビリティが高く、行列のブロック構造を利用して密度汎関数理論の計算に特に適するように設計されています。また、ノード間はMPI、ノード内のSMPにはOpenMPを利用します。スパース行列演算は、プロジェクトで用いたデータでは多くの場合全体の30％以上のコストが掛かかるため、このライブラリの導入により大きな性能向上が望めますが、まだ開発途上です。そのため、最終ベンチマークにのみ結果を示します。
OpenMP実装
効率的にMPI並列のCP2KをOpenMP/MPIハイブリッド並列へ更新させるために、プロファイルに現れた」全体コストに占める割合が高いルーチンを集中的に並列化しました。これには、前回のプロジェクトで報告された、FFT、実空間から平面波への変換、コロケート、積分ルーチンが含まれます。さらに、密行列演算はBLASSライブラリ置き換え(Cray LibSci)を行います。
ハイブリッドコードはpure-MPIコードに比べ、複数の理由からより良いスケール性を示すことが期待されます。まず、同数のコアを用いた場合、コアをOpenMPスレッドに適用する事でMPIプロセス数を削減できます。これは、例えばFFT内のMPIAlltoallv集団通信のようなプロセッサ数に関してスケール性が悪い場合に効果があります。第二に、HPCシステムはコア数を増加させる傾向にあり、MPIのみを用いてノードを飽和させた場合、ネットワーク上の競合が増加します。特にフェーズ2bでは顕著なため、プロジェクトではそのアーキテクチャに密接に対応したスレッド化を行い、良好な性能を出しました。
勿論、通信時間の削減効果もありますが、MPIプロセスとは独立に、同量の分散データを同時に処理する様に効率的な実装をしなくてはなりません。これはいつも直截的に可能なわけではありません。OpenMPとMPIを組み合わせる場合には、MPI実装側からの要請である種のスレッドセーフ性が必要になります。このプロジェクトの実装はシンプルに、OpenMP並列領域の常に外側でMPI実装をするようにしました。これはMPI THREAD FUNNELEDモデルに相当し、MPIの呼び出しはマスタースレッドによってのみ行われることが保証されるものです。
·FFT
CP2Kでは、3D-FFTがQuickstepアルゴリズムの重要な要素です。ここでは、電子密度やポテンシャルの実空間表現をフーリエ空間（平面波）表現への変換とその逆を意味します。その計算は、一般下記のようなに5つの段階に分かれます（2次元分散平面はグリッドの場合。1次元の場合はより少ない）。詳細な議論はJagode 2006[7]などを参照してください。
１．1D-FFT実行（例えばZ方向）
２．MPI_AlltoallvでデータをY方向へ変換
３．1D-FFT実行（Y方向）
４．MPI_AlltoallvでデータをX方向へ変換
５．1D-FFT実行（X方向）
この処理にはOpenMP化が可能な部分が2つあります。一つは1DFFT自身もう一つはMPI_Alltoallvを用いたバッファのパッキング/アンパッキングです。
様々な並列パターンの中で、長さNの1DFFTをシングルスレッドにして、複数FFTのブロックM個を並列化する方法が最も効率的でした。シングルスレッドに比べて4スレッドで2.9倍に高速化しました。
また、データ転置に用いるMPIバッファパッキング処理も並列化しました。125x125x125グリッドに対し、512コアを用いたPure-MPI版に比べて2048コアによる4スレッド版は2.8倍に高速化しました。
·実空間から平面波への変換
これは、実空間グリッドから平面波グリッド間の変換に必要なハロスワップとalltoall通信から成り、平面波グリッドを逆格子空間へ変換するFFTの準備のための処理です。このルーチンで最も時間が掛かるのは通信で、OpenMPで並列化可能なループが複数あります。MPIバッファは並列化されており、各スレッドはバッファの一部をパッキングします。
この処理はOpenMPのworkshare指示行で容易に並列化可能で、その場合には各スレッドは特定の配列のサブセットを複製します。gfortran4.5.0の実装は、配列操作のシリアライズを単一スレッドで行っています。性能面の考慮から、ここでは、スレッドのIDと配列の境界を元に範囲を定義して主導で並列化しました。また、その他の並列化対象ループとして、プロセス間の送受信データ量の計算処理を並列化しました。
最後に、alltoallvがスパース（少数のプロセスペア間通信しかない）のため、ノンブロッキング1対1MPI通信へ置換えました。これにより5～10％の性能向上が得られました。さらに、全ての変更を組み込んだコードに対して、125x125x125グリッドの場合、6スレッドを用いて1152コアまでスケールしました。
·コロケーションおよび積分
3つ目の並列化は、コロケーションと積分です。これらはGaussians（原子中心のガウス基底関数の積）のリストを取り、これをフーリエ変換の前に通常のグリッドにマッピングします。逆もまた同様です。このステップがQuickstep dual-basisアルゴリズムの重要な部分を構成する方法については、[1]を参照してください。
積分は24スレッドで20倍まで高速化しました。コロケーションでは、高価なグリッド・リダクション・ステップのため、前述のリモート・メモリ・アクセスおよび同期のコストのために、最大6コア（シングル・プロセッサー）までしか高速化できません。
この他、相関関数計算サブルーチンについてもOpenMP並列化を行いました。 今回はPBE機能[9]のみを対象としました。関数計算の大部分は、実空間密度グリッド上の点に対して単一のループとして実行されます。 各点で関数計算が実行され、その結果がグリッド上の対応する点に書き込まれます。 このループは各反復が完全に独立しているため、並列化可能です。非常に優れたOpenMP並列効率（6スレッドで効率93％、ノードで24コアすべてを使用して74％）が示されました。
ベンチマーク結果
このセクションでは、3つのCrayシステム、HECToRフェーズ2a、Rosa、HECToRフェーズ2bでのベンチマーク結果を示します。
3つのシステムを比較すると、ノード内のコア数（HECToR 2aでは4、Rosaでは12、HECToR 2bでは24）が増加するにつれて、コードのスケーラビリティが低下します。Pure-MPIコードはそれぞれ256、144、および144コアで最高性能が達成されています。ただし、Seastar 2+ネットワークインターフェイス（3つのシステムすべてで使用）は、間もなくHECToRフェーズ2bの新しいGeminiインターコネクトに置き換えられます。これにより、メッセージスループットが向上し、特に小メッセージのパフォーマンスが向上します。
第2に、pure-MPIコードのパフォーマンスが約1000コアで40?70％向上したことが示されました。
第3に、スレッドを使用することが実際にコードのスケーラビリティを向上っています。 使用するスレッドの適切な数は2?6（単一プロセッサ内のコア数）です。 ハイブリッド並列モードでOpenMPを使用する場合、HECToRフェーズ2aおよびRosaではコード全体のピーク性能が約30％向上しました。
また、HECToRフェーズ2bでは60％増加しました。これは、SeaStarを通過するメッセージ数が大きく減少したためです。
ノードあたり1つのMPIタスクのみを使用し、12個のOpenMPスレッドを使用すると、最高性能はpure-MPIコードの2.5倍になります。タスクあたり2スレッドを使用した場合は、576コアで20％の高速化が示されました。
文献
| [1] | Improving the performance of CP2K on HECToR: A dCSE Project, Iain Bethune, 2009, http://www.hector.ac.uk/cse/distributedcse/reports/cp2k/cp2k _nal report.pdf |
| [2] | CP2K website, http://cp2k.berlios.de |
| [3] | 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 |
| [4] | HECToR: UK National Supercomputing Service, http://www.hector.ac.uk |
| [5] | TOP 500 List - June 2010, http://www.top500.org/list/2010/06/100 |
| [6] | CSCS Swiss National Computing Centre, http://www.cscs.ch |
| [7] | Fourier Transforms for the BlueGene/L Communications Network, H. Jagode, 2006, http://www2.epcc.ed.ac.uk/msc/dissertations/dissertations-0506/hjagode.pdf |
| [8] | The Design and Implementation of FFTW3, M. Frigo and S. Johnson, Proceedings of the IEEE 93 (2), 216.23 1 (2005). |
| [9] | Generalized gradient approximation made simple. J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett., 77 (1996) 3865-68 |