＊ここに掲載するのは、エジンバラ大学、EPCCのI. Bethune博士によるHECToRレポート「Improving the performance of CP2K on HECToR A dCSE Project, I. Bethune, EPCC, The University of Edinburgh, July 29, 2009」を要約したものです。
概要
ロンドン大学ユニバーシティカレッジのベン・スレーター博士研究グループのメンバーであるMatt Watkins博士の前回の作業に続き、チューリッヒ大学のSlater博士とJoost VandeVondele博士がHECToRでの更なるコード最適化に関するdCSE支援の申請を行い、2008年8月から2009年7月の間、EPCCのアプリケーション・コンサルタントであるIain Bethune博士によって6ヶ月の開発が行われました。
·CP2Kについて
CP2K[1]は、GLPライセンス下でフリーに入手可能な密度汎関数法（DFT)パッケージです。このコードは多くの個人とグループにより開発され、現在（2009年6月時点）は50Kラインを超えるFortran95言語で記述されています。主要な開発者は、チューリッヒ大学、物理化学研究所のJurg Hutter教授グループで構成され、特にプロジェクト全体を通してVandeVondele博士が密接に協力しています。
CP2Kと他のDFTコードの大きな違いは、Quickstepアルゴリズムの実装です。これは、波動関数に原子を中心とするガウシアン関数と、電子密度に対してはグリッド上の平面波を用いるデュアル基底を持ちます。このアルゴリズムの詳細については、Quickstepに関する原論文[2]を参照してください。特にこの方法は、二つの表現間の変換について効率的かつスケーラブルな方法が必要であり、これがプロジェクトの大きな目的となっています。CP2Kは他の多くのDFTコードと同様に、実空間から逆格子空間への変換に3次元FFTを用います。
このプロジェクトでの修正は、レグレッションテストを通過してメインCVSレポジトリに登録され、ユーザにダウンロード可能な状態になっています。
·HECToRについて
プロジェクトは、英国国立スーパーコンピュータサービスのHECToRで行われました。このシステムの詳細はウェッブサイト[3]から参照できます。プロジェクトはHECToRのフェーズ1構成で行われました。これは、Cray XT4ノード：5664個、ノード当りにAMD 2.8GHz dual-core Opteronプロセッサーとコア間共有メモリー：6GB RAMを搭載します。プロジェクトでは、この他に通常のデスクトップマシンやクラスターシステム、インテルやAMD、PowerプロセッサベースのHPCシステムでも行いました。
HECToRへの移植
CP2Kは最低限、Fortran95コンパイラ、BLAS、LAPACK、および並列版にはMPI-2、ScaLAPACKライブラリが必要です。オプションとして、ハートレーフォック交換項（HFX)の計算に必要なFFTライブラリ、libitライブラリ[4]も用いられます。Portland Group (PGI), Pathscale, GNU gfortranを用いてコンパイルしました。
CP2Kは、分子動力学、座標最適化、基準振動モード解析、フリーエネルギー計算、経路積分、モンテカルロ計算など、Quickstep DFT,MM,QM/MM等に拠る様々な力の計算手法を用いた多くの異なる利用方法を持ちます。このプロジェクトでは、多くの利用者に実行され、Quickstepアルゴリズムの主要コンポーネントの検証のために代表される、2つのベンチマーク系を用います。
ベンチマーク"bench 64"は、辺12.42Åの立方体内の64個の水分子の周期境界条件で構成され、TZV2P基底関数(40,000基底）を用いた、温度300KでのNVEアンサンブル計算です。系は0.5fs間隔の50タイムステップを時間進行させます。このシミュレーション実行時間は256コアで数分です。
ベンチマーク"W216"は、1辺34Åのセル内に216個の水分子を含み、さらに大きなmolecularly optimized basis[5]を用います。ですが、この系は周期系ではなく、原子群はシミュレーションセルの中心に配置されます。このため後に議論する興味深い負荷バランス性を生じます。0.5fs間隔の10タイムステップ計算時間は、1024コアで約20分です。
CP2Kには、単体テストのようなコードの特定の機能を分離する一連の"libtests"も含まれています。libtestsの"RS_PW_TRANSFER"（実空間から平面波への変換）や"PW_TRANSFER"（FFT）はプロジェクト中の検証に広く用いました。
実空間から平面波への変換
ここではルーチンrs2pw_transferを最適化します。CP2Kは、2つの異なる基底、波動関数を表現するガウシアン関数の積（スパース行列に保存されます）と、平面波（グリッド上に保存されます）を用います。他のDFTコードと共通して、CP2Kは系の基底状態エネルギーを求めるために自己無撞着場（SCF)計算を行います。各繰り返しにおいては、この2つの表現の間で、以下の変換と逆変換が必要です。
- ガウス関数の積に対応する行列要素は、一組の実空間マルチグリッドに写像される。マルチグリッドの各レベルは完全に分散されているので、各MPIタスクはグリッドの特定の部分を中心とするガウス分布のサブセットを写像するか、複製されます：この場合、タスクは任意のガウスをマッピングできます。これらは、負荷バランス時に重要な性質です。このステップは、「コロケーション」と呼ばれ、逆は「インテグレーション」です。
- その後、実空間マルチグリッドに格納されたデータは、一連の平面波マルチグリッドに転送されます。 平面波グリッドは一般に（並列FFTに合わせてスライスまたはペンシルのいずれかの）実空間グリッドに異なる方法で配布されるため、データを分散された平面波グリッド上の適切な場所へ、MPI_Alltoallvを使用してグローバルに再シャッフルします。 これに加えて「再スワップ」ステップが必要です（実空間から平面波への再分配の前に、またその後の平面波から実空間への逆方向時に）。 これは、ガウシアン関数がグリッドにマップされている場合に、マッピングを実行したMPIタスクのローカルグリッドの境界を超えて拡張する可能性があるために必要となります。 したがって、すべてのプロセスは、可能な限り最大のガウシアンを満たすのに十分広いハロ領域を維持し、マッピングが完了した後、これらのハロー領域を近傍とスワップしてローカルグリッドに合計します。
- データが平面波グリッド上に配置されると、並列3Dフーリエ変換が行われ、実空間から逆格子空間へ変換されます。
ハロ領域スワップ
ハロ領域スワップは、特にMPIタスク数が増えると極めて高コストです。固定グリッドでは、グリッドがプロセスで分割されるほどローカルなグリッド空間は小さくなります。しかしながら、ハロ領域幅は最大のガウシアン関数サイズで決まるため、プロセッサ数が増えればハロ領域はローカル領域に比べてより大きくなります。例えば、125x125x125グリッドは512プロセッサ上であれば、ローカル領域サイズは16x16x16、ハロ領域幅は18です。化火kとしてハロ領域は全グリッドデータの97%に達します。
既存のアルゴリズムでは、ハロは最初にX方向の近傍と交換され、次にY、最後にZ方向の順番です。それぞれにおいて、全ハロが送信されます。ハロスワップが完了した後は、ローカルグリッド内のデータのみが必要であるため、ハロの全幅を各方向に送信する必要はないということが判りました：つまり、ハロになるデータは送信前に破棄することができます。3次元の場合、送信されるデータ量とバッファリングに費やされる時間の両方に関して、実質的な節約が成されます。
また、このままだと、XバッファのパキングとアンパッキングのコストはYやZバッファよりも大きくなります。これはメモリー上のハロデータが、Zバッファでは連続ですが、XやYバッファでは断片化しているためです。これを最適化することで、100%の高速化が達成されました。
これ以外に領域分割の改善を行いました。既存の分割では、グリッドはローカルグリッドと相対的にハロ領域の全サイズを最小化するように決められていました。そこで、重み付けを用いて、Z方向の分断数を最小にするように分割するようにしました。つまり、この時、ハロサイズは最大ですが、スワップコストはより小さくなります。
バッファパッキング処理は相対的に高コストなので、MPI派生タイプの利用も検討しましたが、効果はありませんでした。これは主にハロデータが実空間グリッドへ再統合される際に足し込みを行うため、受信バッファが必要なためです。
ノンブロッキング通信
rs2pw_transferルーチンへの2番目の修正は、MPI_SendRecvをノンブロッキングMPIへ置き換えることです。これにより、バッファパッキング処理を通信とオーバーラップさせることが可能になります。ここで、事前にrecvを呼出し、可能な限り送信を開始することにより、他のバッファのパッキングと重複できるようにしました。
しかしながらこの修正は、送信時間がパッキング時間よりも時間が掛かるため、大きな効果は得られませんでした。
FFT
ここではFFTを最適化します。CP2Kは、FFTW[7]2および3, ESSL from IBM, ACML from AMDとのインターフェイス およびGPUに対応したCUDA実装、またGoedecker et. al.[6]の方法に従ったビルトインFFTライブラリも備わります。これらは全てシリアル実行版です。並列化には、並列3DFFTと、FFT演算の合間に発行するMPI_Alltoallvによるデータのグローバル転置操作が必要です。プロセッサ数が小さい場合は、CP2Kではスラブ領域分割による一回の転置で済みます。大規模プロセッサでは、ペンシル分割を用いて2回の転置操作を行います。これは高コストですが、スケール性はより高いものです。
Alltoallvは、送受信で異なるデータ量を指定することが可能ですが、実際には通常は1列分の違いしかありません。Intel MPI benchmark[8]によれば、HECToR上では同じ量のデータを転送する場合には、AlltoallvよりもAlltoallの方が良い性能を持っています。転置されるデータ量は通常256KBから数MBであり、Alltoallへ変更することで20〜30%の性能向上が期待できます。これを行うには転置の前後でデータに余分なパディングを追加して、同じ量を送信するようにしました。
この修正はlibtestでは43％高速化を示しましたが、全体性能のベンチマークでは高々2％の向上に止まりました。
また、FFTWのプラン変更を試しましたが、大きな性能向上には至りませんでした。
負荷バランスの改善
負荷バランスの改善で特に重要なのは、スパース行列表現のガウシアン基底関数が実空間グリッドへマッピングされる部分です。
実空間から平面波への変換とFFTは、プロセス当たりの作業量がグリッド数に比例しており、均等にスラブ、ペンシル、ブロックへ分散されるため、負荷バランスは良好です。しかしながら、ガウシアン関数を行列からグリッドへマッピングするのはそれほど良いバランスを持ちません。この行列は、ScaLAPACKのブロックサイクリック形式で全プロセッサへ分散されますが、一般に、プロセスの行列ブロックに含まれるガウシアン関数は、実空間グリッドの自身のセクションに配置されない可能性があります。
CP2Kのコストモデルはきわめて正確に機能します。均一な系の場合は各プロセスグリッドセクションの同数の原子が配置され、極めて良好な負荷バランスが観察されます。
しかしながら、界面やクラスターのような不均一系では、あるプロセスでは、そのローカル実空間グリッド内の原子数が少なく、実行するタスクも小さくなる場合があり、その場合には負荷バランスが劣化します。
この問題に対処するために、各プロセスが、各グリッドレベルの実空間グリッドの同じセグメントを所有する代わりに、グリッドブロックのプロセスへのマッピングを各レベルで変えることができるスキームを考えました。1つのグリッドレベルで重い負荷になるブロックを負荷の低いブロックとペアにして、負荷をより均等に分散させる様にします。
W216ケースでベンチマークを行ったところ、128コアで25％、2048コアで18％の実行時間の高速化が示されました。
最終ベンチマーク
全体の実行時間は、bench64ベンチマークでは256コアで30％、W216ベンチマークでは1024コアで300％加速しました。
文献
| [1] | CP2K website, http://cp2k.berlios.de |
| [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. Chassaing and J. Hutter, Comp. Phys. Comm. 167, 103 (2005) |
| [3] | HECToR website, http://www.hector.ac.uk |
| [4] | Libint website, http://www.files.chem.vt.edu/chem-dept/valeev/software/libint/libint.html |
| [5] | Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J.VandeVondele and J. Hutter, J. Chem. Phys. 127, 114105 (2007) |
| [6] | An efficient 3-dim FFT for plane wave electronic structure calculations on massively parrallel machines composed of multiprocessor nodes, S. Goedecker, M. Boulet and T. Deutsch, http://pages.unibas.ch/comphys/comphys/SOFTWARE/FFT/fft.ps |
| [7] | The Design and Implementation of FFTW3, M. Frigo and S. Johnson, Proceedings of the IEEE 93 (2), 216.231 (2005). |
| [8] | Intel® MPI Benchmarks 3.2, http://software.intel.com/en-us/articles/intel-mpi-benchmarks/ |