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

＊ここに掲載するのは、HECToRレポート「Adapting QSGW to large multi-core systems」を要約したものです。

概要

準粒子自己無撞着GW近似（Quasi-particle self-consistent GW、QSGW）法は、電子構造の非経験的手法の中で標準的な手法の一つとして注目されつつあります。これは、最適化された単一粒子レベルと多体レベルでの電子構造を同時に扱う手法です。よって、通常のKohn-Sham密度汎関数理論（KS-DFT）のハイエンドの代替手法となります。この方法は、KS-DFTの良く知られた欠点の大部分を乗り越え、LDA + DMFTのような現在普及している方法に固有のあいまいさのないクリーンでパラメータを不要とするやり方で、単粒子理論と多体理論との間のギャップを橋渡しする方法です。支払われる対価は計算コストです。QSGWは、LDAより約100?1000倍のコストが必要です。
このプロジェクトの当初の開発段階では、コードは本質的にシーケンシャルなコードであり、物理的に興味のある系の実行時間は月単位でした。 日常的な計算でQSGW計算を実行可能にするためには、コードを並列化することが重要でした。

ループ構造と並列化手法

QSGW計算において必要となる主たる量は分極率Πです。ループは計算対象の波数ベクトルqと周波数ωから成り、その内側で、波数k、バンドインデックスn（占有）およびn'（非占有）に関する総和のループを持つネストループ構成です。これらをネストしたMPI並列化で処理します。自己エネルギー計算のループも似たループ構造です。
最外側のqループは、計算が互いに独立しているため、グローバルデータを全てのノードにブロードキャストすることを除けばほとんど通信を必要とせず、直截的に実装できます。 出力ファイルは可能な限り、与えられたqベクトルのプロセッサグループの主ノードでローカルに書き込むようにすることが可能でした。
内側のkの総和ループには、リダクション操作を追加しました。 さらに、負荷の高い行列要素計算は、OpenMPで並列化しました（主に基底関数による）。これらの3つのレベルの並列化に関するスケーリングの結果を以下に示します。

実装

· MPI並列

ネストループqとkの並列化においては、N個のMPIプロセスのプール（MPI_COMM_WORLDで記述）をN/Nqプロセスを持つqグループに分割しました。 次に、qベクトルのループをNqグループ上で並列化しました。各qグループ内のkポイント上の並列化は同様に、 qグループのプロセスをさらに分割してkグループのループを並列化しました。各プロセッサでは、各kベクトル毎に結果を総和します。
kグループの実装では、与えられたk点ごとにさらなるレベルのMPI並列化（例えば、基底関数の上で）を導入する柔軟性を持ち合わせています。ですが実際のプロジェクトの現段階では、kグループのサイズは1、すなわちkグループの数Nkは各qグループ内のプロセスの数に等しい状況で行いました。

ここでは、様々なレベルの並列化を処理することを目的としたMPIモジュールを作成しました。このMPIモジュールは、MPIプロセッサグループのデータ構造と、特定のプロセッサグループを複数のサブグループに分割するルーチンを提供可能なものです。このモジュールは非常に一般的に書かれており、ネストされたレベルのループ並列化をサポートします。 現在のプロジェクトでは、2つのレベルが使用されています。

·I/O並列化

既存コードの多くの入出力ファイルは、qポイント毎に別々のファイルとして書き込まれていました。このような場合、I/Oを各qグループの主ノードで行うのが極めて自然な実装です。
既存コードでは、qベクトル毎にスクリーニングされたクーロン相互作用を含むダイレクトアクセス出力ファイルには、全ての周波数ωをレコードとして含んでいました。系のサイズによっては、各ファイルのサイズが約1 GBに達することがありました。ファイル出力を加速するために、各レコード、すなわち各周波数を異なるプロセッサで書き込み、その後、元のフォーマットにファイルを連結する様に修正しました。

·OpenMP並列化

MPI並列化はトップダウンの考え方で行いましたが、OpenMPはボトムアップで実装します。基底関数とブロッホ関数との間の行列要素の計算に使用される行列乗算のような、コードのより低レベルでの単純ループの並列化に適用しました。

プロジェクトの成果

·作業1：並列化

分極関数（q、ωについては直截的、kについてはリダクション実装）：2つのルーチンが並列化されました。
自己エネルギー（qについては直截的、k、、ω'についてはグローバルリダクション実装）：２つのルーチンが並列化されました。
最終的にHECToRの2000コアまでスケールしました。

·作業2：メモリー効率化

基底関数並列化と行列の分散保存：分散メモリモデルを実装し、大型モデル用のΠij行列の処理を実装します。ループ構造は再配列され、分散メモリ、マトリックス乗算と逆行列をハイブリッド並列化します。さらに、数学ライブラリの呼び出しを含まない他の内部ループをOpenMP並列化します。このレベルの並列処理は、第1レベルのスケーリングには影響しません。
ベアおよびスクリーニングされたクーロン相互作用：両方のクーロン相互作用を生成して 分割された形式で保存します。逆行列計算は、系のサイズに応じて異なる並列アルゴリズム（ScaLapack、PLASMA、ダイソン方程式による反復解）を使用します。
自己エネルギーは上述と同じループ構成にします。行列積には並列ライブラリーを用います。
作業時間の制約のため、分散メモリーコードはOpenMP並列のみ実装されました。

スケール性能結果

QSGWで最もコストのかかる、感受率とスクリーニングされたクーロン相互作用、および自己エネルギー計算についてプロファイル分析を行いました。
テストには、16Fe原子のスーパーセルを用いました。既存コードでは12コアで数か月かかる計算でした。この計算では、10x10x10のqポイントを用いています。これは125既約ポイントに相当します。
テストでは5x5x5のkポイントを用いました（10既約ポイント）。

·自己エネルギー計算（exchange部）

外側qループに関するスピードアップは、qベクトルの数（ここでは10）までほぼ線形です。この直線的な挙動は、qポイントの増加に伴い追随することが予想できます。kベクトル上の並列化は、より早く、主にqループに存在しない通信のために劣化します。最適なOMPスレッドは8です。

·スクリーニングされたクーロン相互作用

感受率スクリーニングされたクーロン相互作用を計算するコードは、QSGW中で最も計算集約中的な部分です。4000コア以上まで線形にスケールします。2640コアにおいて1380倍の高速化が示されました。計算はまだ飽和状態には達しておらず、コア数が多いほど性能がさらに向上することが期待されます。
また、qループとkループの両方は、その計算負荷の高さから、独立に50MPタスクまでほぼ完全なスケール性を示すことが確認されました。

·自己エネルギー（correlation部）

これは、周波数ループが追加されることで、exchange部よりもさらに計算負荷が高い部分です。3000コア以上までほぼ線形にスケールしました。
また、qループとkループの両方は、その計算負荷の高さから、独立に30MPタスクまでほぼ完全なスケール性を示すことが確認されました。

謝辞

このプロジェクトは、nAG Ltd.が運営するHECToRの分散計算科学および工学$CSE$サービスの基に実行されました。英国の国立スーパーコンピューティング・サービスである、HECToR：英国リサーチ・カウンシル・ハイエンド計算サービスは、リサーチ・カウンシルを代行するEPSRCが管理しています。そのミッションは英国学術界の科学および工学の研究支援です。HECToRスーパーコンピューターは、UoE HPCx Ltd.およびnAG Ltd.のCSEサポートサービスにより管理運営されています。