＊ここに掲載するのは、ユニバーシティカレッジ・ロンドンのLianheng Tong博士によるHECToRレポート「dCSE Project on Improving Scaling in Conquest for HECToR Phase 3-OpenMP-MPI Hybrid Implementation, Lianheng Tong, University College London, 2012/09/27」を要約したものです。
概要
材料科学や生物化学における様々な重要な問題に対して、その物質や分子の特性に関する観測量の理解には大規模系に対する非経験的量子力学研究が必要とされます。例えば、Si 001表面上のGe原子のhut-cluster自己組織化[6]の理解には、少なくとも数層のSi表面と一つのGe hut-clusterを含むシミュレーションが必要です。hut-clusterは20nm正方でSi基盤を含め2万以上の原子に及びます。生体化学分子の研究においては、系は周期性を持たず、小さなシミュレーションセルを多く用意する方法では困難です。処置済みのDNAベースペアは既に数千原子を含んでいます。
Conquest[1]は、このような大規模系の非経験的研究に理想的な線形スケール密度汎関数法(DFT)[5]コードです。HECToRのようなHPCマシン上での実行にはスケール性の良さが本質的です。よって各コア上の仕事量を一定に保ちつつ、コア数を増やすことで、そのスケール性に従って計算対象の系のサイズを大きくすることが可能です。こうすることで、理想的には多くのコアを用いた大規模計算では実行時間の増加を抑えることが出来ます。この性質は「弱スケール性（weak-scaling）」と一般に呼ばれています。
Conquestは、ノード当たり4コアのHECToRフェーズ2a上でほぼ完ぺきな弱スケール性を示しますが、ノード当たり24コアのフェーズ2bではその性能は劣化します。この問題はノードを「不飽和(ノード内のすべてのコアを用いない)」にして実行すると回避できます。例えば、ノード当たり4コアのみを用いるようにすると、弱スケール性が回復します。この事から示唆されることは、この新しい世代のCPUをMPIタスクで飽和させると、ノード間通信チャネルが有限であることから性能劣化が発生することです。全コアが、ノードを跨いで相互に通信する場合は必然的に渋滞します。ノードを非飽和にすれば、MPIプロセス当たりのノード間チャネルを増やすことができます。
HECToRフェーズ3は、2つの16コアCPU、言い換えれば32コアのノードを持ちます。つまりConquestのスケール性にとっては更に状況は悪化しています。これはConquestだけの問題ではなく、新世代のマルチコアCPU上の(純粋な)MPIコードに共通する問題です。さらにスケール性の問題に加えて、pure-MPIコードにはメモリー効率の問題があります。新世代のCPUはノード当たりのコアを増やす傾向を持ち、コア当たりのメモリーはより少なくなります。HECToRフェーズ3のノードは、4つのNUMA領域で構成された全32GBRAMを持ちます。これが意味するのは、飽和ノードで実行した場合は、各MPIプロセスはメモリーは1GBまでしか利用できないということです。フェーズ2aではノード当たり4コアで全8GBRAM、フェーズ2bでは24コアで全32GBでした。つまりもしフェーズ3で全てのコアを用いて実行する場合は、コア当たりのシミュレーションサイズを減らすか、より粗いシミュレーショングリッドを使う必要が生じますが、実際にはこれらは共に満足できる解決策ではありません。前者は、各コアの継続的な作業量が削減されるため、計算の非効率性を招きキャッシュミスが増えます。後者は精度の劣化を招きます。
この問題に対して、コードを変更しない簡単な回避策は、ノード内でより少ないMPIプロセスを使うことです。しかしこれは利用しないコアを多く残すことになり、ノード数のみ多く必要となる「高価」な方法です。こうして、HECToR上の多くのpure-MPIコードはハイブリッド並列へ移行しつつあります。これは、MPIの分散メモリー並列(DMP)とOpenMPのような共有メモリー並列(SMP)をミックスした構造です。このやり方を用いれば、ノード当たりにより少ないMPIプロセス利用が可能で、ノード間通信の競合を減らし、プロセス当たりより多くのメモリーを利用することが出来ます。同時に、ノード内の残りのコアをOpenMPスレッドを用いてMPIプロセスの作業を分担させます。ただしこの場合はマスタースレッドだけがMPI通信を担当します。
このプロジェクトの目的は、フェーズ3上のConquest性能プロファイリングとコード上のホットスポットの特定を行い、OpenMP/MPIハイブリッド並列を実装することです。
プロファイリング
HECToRフェーズ3上のConquestの性能測定用のベンチマークデータは、バルクシリコンの1点計算：セルフコンシステント密度汎関数基底状態エネルギー計算、積分グリッドは0.32Bohr、シングルゼータ基底関数(シリコン原子当たり4基底)を用います。
用いたコンパイラとライブラリは、Cray compiler suite 4.0.46、Cray LibSci 11.0.06です。開発者と研究者が最も関心があるスケール性は弱スケール性です。これは、各コアの作業量を固定した場合に、計算終了までの時間が系のサイズによりどう変化するかを示すものです。系のサイズは、計算に用いるコア数を増やすことで大きくなります。
ノードを飽和されて測定すると、512原子/1ノードでの計算時間と比較して、32768原子/64ノードの計算時間は約2.1倍に増加します。
しかしながらノード当たり4MPIプロセスで計算すれば、同じMPIタスク数(32)の512原子/8ノードの実行時間は855秒から552秒まで縮小し、32768原子の実行時間比は1.875倍まで高速化します。
1飽和ノードと8不飽和ノードでのプロファイリングから、どちらも同じ程度のMPI通信比率が示されています。
にもかかわらず、飽和ノードは系のサイズを増やすと通信時間が増加します。
実際に、32768原子では、飽和ノードではMPI通信時間が全体の53％も費やしてボトルネックになっています。
これらの結果から、OpenMP/MPIハイブリッド並列がフェーズ3では有効であることが示唆されます。
また、プロファイリングから、行列演算カーネルm_kern_minとm_kern_maxがコード内で最も時間のかかる部分であることが判りました。これらは行列積とグリッド上の積分を行うコードです。また、行列積がほとんど(70%)を占めています。よってSMP並列の対象をmultiply_moduleのm_kern_minとm_kern_maxとします。
DGEMMの適用
そこで、m_kern_maxとm_kern_minの行列積演算にマルチスレッドBLASSライブラリのDGEMMを適用しました。
512及び2048個のシリコン原子系を用いて32MPIプロセスでテストしました。512原子系については、シングルゼータ基底(SZ)およびSZ+ポーラリゼーション（SZP） 基底を用いました。テストは、ノードを飽和する様に強スケール性テストを行いました。NUMA領域は8スレッドまでなので、これを超えると通信が発生して性能が劣化します。しかしながら、DGEMM呼び出しがループ再内側だったため、スレッド生成オーバーヘッドが65%を占め、性能は向上しませんでした。一方、シングルスレッド版は既存コードよりも数％高速でした。
行列積カーネルへのOpenMP実装
カーネルには3重ループがあります。最外側は原子、中間はその隣接原子、再内側は両方の隣接原子です。一般的には、スレッド生成とバリア同期のオーバーヘッドを最小にするために最外側ループをOpenMP対象にします。
乗算カーネルのサブルーチン全体を囲むように並列領域を作ることは、まったく問題ありません。殆どの変数は入力データのみなので、並列領域でこれを安全に共有することが可能です。演算結果を加算する場合のみ、競合を避ける注意が必要です。
既存のループにはインデックスのインクリメントが実装されていますが、並列化すると正しく動作しなくなります。よって並列動作前にインデックスを一時的な配列に保存するようにします。
バルクシリコンケースでのスケール性能
OpenMP/MPIハイブリッド並列性能を、バルクシリコンについて非経験的セルフコンシステントDFT基底状態エネルギーの1点計算で行いました。
計算はコア当たり同数の16原子、つまりノード当たり512原子を割り当て、全て飽和ノードを用いました。Pure-MPI計算はノード当たり32MPIプロセスで性能が飽和します。
弱スケール性能において、Pure-MPI並列性性能は小さな系では勝りますが、大きくなるにつれてハイブリッド性能が優位になります。
1万原子以上では、ノード当たり16MPIプロセス、MPIプロセス当たり2スレッドが最良の性能でした。
強スケール性能テストでは、MPIプロセスは固定してスレッド数を増やしました。8スレッドで1割程度高速化されました。
また、OpenMPスケジューリングオプションの中では、STATICが最も良好な性能を示しました。
氷ケースでのスケール性能
Conquestは生体分子系にも用いられます。このような系が面倒な点は二つあります：空間内で原子が不均一に分散され、基底関数に対して異なる相互作用を持つことです。最初の問題は負荷バランスを取ることで解決されますが、このプロジェクトのスコープではありません。第二の問題は、OpenMPスレッドに負荷分散を生じます。
このような系の性能をテストするのに、バルクの氷を用います。これは負荷バランス手法を必要としません。ただし、水素と酸素には異なる基底関数を要し、氷の構造はスパースで、隣接リストに不均一性を導きます。ここでは短いカットオフ（水素に4.8Bohr、酸素に3.9Bohr）を用いて、隣接リストをさらに不均一にしました。
·強スケール性
1536原子系で、32MPIプロセス固定でスレッド数を変化させて強スケール性のテストしました。結果はスレッド数が4を超えると改善率はフラットになりました。この性能劣化は、氷の特性として、隣接リスト長が短いことによるものでした。バルクシリコンと同様に、大きなサイズでは大きく改善します。
·弱スケール性
小さい氷系では純粋なMPI計算が高速で、ノード当りのMPIプロセスが多いほど計算が速くなります。しかし、OpenMP-MPIの実装では、スレッド数が増えるとスケーラビリティが向上し、システムが大規模になると高速化になります。バルクシリコンの場合と同じように、OpenMP/MPIハイブリッド実装は、隣接ノードがよりスパースであるほど弱スケール性を改善します。
謝辞
このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学(CSE)サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。
文献
| [1] | http://www.order-n.org. |
| [2] | D. R. Bowler. Linear-scaling density matrix minimisation and electron number conservation. Technical report, University College London, 2005. |
| [3] | D. R. Bowler, T. Miyazaki, and M. J. Gillan. Parallel sparse matrix multiplication for linear scaling electronic structure calculations. Comput. Phys. Commun., 137(2):255-273, June 2001. |
| [4] | H. Hotelling. Some new methods in matrix calculation. Ann. Math. Stat., 14(1):1-34, 1943. |
| [5] | R. M. Martin. Electronic Structure: Basic Theory and Practical Methods. Cambridge University Press, 2004. |
| [6] | T. Miyazaki, D. R. Bowler, M. J. Gillan, and T. Ohno. The energetics of hut-cluster self-assembly in ge/si(001) from linear-scaling dft calculations. J. Phys. Soc. Jpn., 77:123706, 2008. |
| [7] | A. H. R. Palser and D. E. Manolopoulos. Canonical purification of the density matrix in electronic-structure theory. Phys. Rev. B, 58(19):12704-12711, Nov 1998. |