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

＊ここに掲載するのは、ヨーク大学のPhil Hasnip博士によるHECToRレポート「Bands-parallelism in Castep, A dCSE Project, Phil Hasnip, August 14, 2008」を要約したものです。

概要

このプロジェクトの目的は、密度汎関数計算コードCastep[1]に対して、実行時の利用可能ノード数を増やして実行効率を向上させる新たな並列手法を組み込むことです。この改善により、数千プロセッシングエレメント$PE$まで効率が向上することが期待されます。

この作業は2007年11月から2008年7月まで行われました。主任調査員はK. Refson博士$RAL$で、M.I.J. Probert $York$およびM. Plummer $STFC$の両博士が協力をしました。Castepの主要な計算は、バンドと呼ばれる固有値の低いところから1％程度を大規模固有値問題で解くものです。現状はバンドを構成する変数上で並列化されています。このプロジェクトの目的は、バンドそのものを用いて並列化することです。作業は以下の4つから成ります：

1. １：基本的な並列化：現状の並列手法に加えて、Castepの主要な計算とストレージをバンド上で分割する。
2. ２：逆行列と対角化計算の分散：現在逆行列および対角化処理はシリアル実行されている。そこでこの処理を可能な限り多くのプロセッサを用いて分割して計算する。
3. ３：バンド別の最適化：既知のアルゴリズムを用いて規格化処理の回数を減らし、速度とスケール性を改善する。
4. ４：新たなバンド最適化手法の並列化と最適化：可能な限り高速かつロバストな新しい最適化手法を組み込み、これを並列化する。

目標は少なくとも8倍多くのノードを用いて効率を上げることです。さらにnAGからの要請で、最適なコンパイラとそのフラグ、および最適なライブラリの調査を行いました。

主要な結果：既存の並列スキームに加えてバンド並列手法が組み込まれました。また、大規模逆行列および対角化も並列化され、バンドとは独立な最適化スキームが2種組み込まれました。
CastepのHECToR上での性能は劇的に改善されました。標準ベンチマークal3x3では、Castep4.2に比べて4倍多くのコアを用いてスケールしました。

HECToR上でのCastepの性能調査

· 既存コードの性能

Castepの実効性能は、通常は対角化とFFTで決まります。対角化（および部分空間対角化）は、HECToRに備わるACMLあるいはCray Libsciで提供される標準的なBLASとLAPACKを用いて実行されます。Castepはポータビリティの考慮から内製FFTを持ちますが、FFTWのようなチューニングされたライブラリとは比べ物になりません。
CastepはFortran90で記述され、HECToRはFortran90コンパイラとしてPortland Group $pgf90$と、Pathscale $pathf90$およびGNU gfortranを利用可能です。初期の調査の後に、Pathscaleのpathf90が選択されました。
また、このプロジェクトで用いたCastepは特に明記しない限り、United Kingdom Car-Parinello $UKCP$ consortiumからリリースされた最新のバージョン4.2です。

· ベンチマーク

ベンチマークデータ：Castepの標準ベンチマークデータは並列計算には小さすぎます。最小のベンチマークは酸化アルミニウムal1x1で、実行時間はHECToRの8PE上で6分です。大きなものでは窒化チタンTiNがありますが、これは16PEで1時間です。これはDMアルゴリズムが収束が遅く、SCFサイクル1回当たりの時間はal1x1に比べ2倍強程度であるためです。このため主なベンチマークとしてal3x3を選びましたが、それでもまだ大規模並列には小さなものです。
al3x3データは、270原子（108Al、162O）、88,184 Gベクトル、778バンド（1296電子、非スピン分極、および130伝導帯バンド）、2k点（対称2x2x1MPグリッド）です。ほとんどの計算では、Castepのパラメータファイルal3x3.paramには、最適化パラメータopt_strategy_bias:3を追加しました。

FFT：HECToRでのFFTベンチマークには、FFTWバージョン2とpgf90コンパイラを用いました。CrayのLibSci $10.2.1$にはBLASとLAPACKを用いました。また、サブルーチンwave_recip_to_real_sliceおよびwave_real_to_recip_sliceの測定にはCastep内蔵のトレースモジュールを用いました。このサブルーチンは固有状態の1グループ（波動関数スライス）を扱い、それらを逆格子空間から実空間へ、またはその逆へフーリエ変換します。FFTW 3.1.1はHECToR上で最高速のライブラリです。

数学ライブラリ$BLAS$：Castep内では多くの時間が、倍精度の行列-行列積のサブルーチンZGEMMで費やされます。規格化および部分空間（回転）変換操作は、共にZGEMMを用いた波動関数のユニタリー変換を行いますが、non-local projectorと呼ばれる計算でも多く用いられます。
規格化と対角化サブルーチンは、メモリーコピーやメタデータの更新といった作業も多く含み、これは小さな系では測定時間を不安定にします。このため、ZGEMMの性能測定ではnon-local projectorオーバーラップの時間、特にサブルーチンion_beta_add_multi_recip_allに着目しました。
結果として、少なくともZGEMMに関してはLibSci 10.2.1のBLASがHECToR上で最高速でした。

ノード効率：HECToRのノードは2コア（PE）を搭載するため、Castepの効率のノード内のPE数への依存性を調べました。最高性能は、FFTW3およびlibSciバージョン10.2.0のGoto_BLASでした。最高効率はノード当たり1コアを使用した場合です。2コアを使用してもコストパフォーマンスは減少しますが、CastepのSMP通信機能[1]を用いると、スケール性は劇的に改善し、ノード当たり1コアのケースに近づきます。

·ベースライン

ベンチマークのベースラインとして、Pathscale 3.0コンパイラ、Cray Libsci 10.2.1およびFFTW3ライブラリを用いたコードの性能を選びました。

·解析結果

非ライブラリ部分とnlpot_apply_precon内のZGEMMに大きな時間が掛かっていることが判りました。非ライブラリ部分の時間はパッキングルーチンからの寄与で、これは、波動関数バンドから非局所擬ポテンシャル射影子への射影を含むアンパックされた配列beta_phiを扱うものです。ZGEMMについては、最初の行列が得るミート行列のため、ZHEMMへ置き換えれば負荷は約半分になります。

バンド並列$作業項目1$

この作業は、N_G個のGベクトルとN_k個のk点に加え、N_b個のバンドに渡って波動関数係数、その他の固有値、バンド占有率および部分空間ハミルトニアンといったバンド間オーバーラップを分散することです。

ここでのプログラミング作業の多くは、分散対象音バンドオーバーラップ行列の構築です。この修正はほとんどがwaveモジュールに集約されました。これは波動関数とバンド上の操作の大部分を扱うものですが、当然他のモジュールにも変更が加えられました。

al1x1ベンチマークを16ノード（4バンド並列および4gv並列）でプロファイリングしたところ、サブルーチンwave_rotateが最も高負荷であることが判りました。これは予想通りで、波動関数の変換コストは系のサイズの2乗に比例し、バンド並列においては通信コストが発生するためです。このサブルーチンを最適化し、各ノードがlog2nodeの通信フェーズを行い、第1フェーズは変換されたデータの半分の交換を含み、後続フェーズは前のフェーズのデータの半分を交換するよう、通信の改良を行いました。

この通信の改良によりサブルーチンwave_rotateの性能は改善されますが、保存容量が増加します。実際に第1フェーズでは、新たに変換されたデータの半分の交換が必要となるため、送信バッファと受信バッファは分散されていない波動関数全体を構成します。バンド並列化は、比較的少数のノード（通常<16）でのみ効率的であるため、これまでのところ問題はほぼ実証されていませんが、将来は単一ノード容量の数倍に制限するほうが賢明でしょう。このような変化は無論、opt_strategy_biasパラメータ値に左右される可能性があります。

最終的なal3x3ベンチマーク結果から、このバンド並列は、50%以上の効率を確保するPE数を221から436個へ増加させました（SMP最適化機能を用いない場合）。

対角化と逆行列の並列化

ここで最初に行うのは、波動関数の規格化と部分空間ハミルトニアン両方で必要になるバンドオーバーラップ行列の計算をバンドグループで並列化することです。この行列は対角化が必要ですが、既存コードではシリアル実行されていました。この操作はN_b³に比例し、より大規模なシステムが研究されるにつれて支配的になり始めるでしょう。 このため、プロジェクトの第2段階では、これらの操作をノードに分散します。

Castepは、部分空間対角化にLAPACKサブルーチンZHEEVを用いています。そのテストベッドとして、乱数のエルミート行列を用いてLAPACKの代わりに並列版のScaLAPACKを用いました。
ここで、さらに改良された並列対角化ルーチンPZHEEVR$このRはMultiple Relatively Robust Representations法を意味します$を用いたところ、PZHEEVの性能を凌駕しました。
ScaLAPACKルーチンはblock-cyclic分散並列をベースにしており、これにより、行または列だけでなく、一般的な方法でデータを配信することができます。

新しい並列逆行列および対角化ルーチンにより、Castep全体の顕著なスケール性能向上が示されました。これは大きなコア数を用いた場合に効果があります。Castep4.2に比べて2から4倍高速化されました。標準ベンチマークal3x3は1024コアで50％の効率を示しました。これは原子当たり3コア強に相当し、より大規模計算ではさらに向上すると予測できます。

個別のバンド計算の最適化

大規模計算では、固有状態の規格化計算がボトルネックになります。これにはバンドオーバーラップ行列の計算と逆行列が含まれ、それぞれN_p N_b²およびN_b³に比例します。ここでN_pは平面波基底数、N_bはバンド数です。さらに、バンド並列モードにおいては個々の固有状態が異なるPEにあるため、前者は通信ボトルネックにもなります。
つまり、顕わな規格化を要しない最適化が必要です。

残念ながら、RMM-DIISと新しく開発したバリエーションはロバストかつ高速ではありません。規格化の修正によりSCF１サイクルの時間は減少しますが、収束までに多くのサイクル数が必要になります。RMM-DIIS法は特に、収束点近傍で厳しい数値不安定性に悩まされました。これは、試行固有状態が真の固有状態に近づくにつれて、残差行列がますます特異になるからです。収束法の変化のみを見るためにCastepを固定された密度行列で実行しました。

これらの結果は、収束計算に多く見られる典型的なものでした。つまり、基底状態に近づけるのは比較的容易ですが、必要な精度を得ることが困難です。こうした状況はアルゴリズム固有のもので、アップデート時に正規直交性を適用することで、両方の方法が迅速かつ確実に収束することができました。

文献

[1]	S.J. Clark, M.D. Segall, C.J. Pickard, P.J. Hasnip, M.J. Probert, K. Refson and M.C. Payne, “First principles methods using CASTEP”, Zeit. f¨ur Kryst. 220$5-6$ $2004$ 567-570
[2]	P.J. Hasnip and C.J. Pickard, “Electronic energy minimisation with ultrasoft pseudopotentials”. Comp. Phys. Comm. 174 $2006$ 24-29
[3]	P. Pulay, “Convergence acceleration of iterative sequences. The case of SCF iteration”. Chem. Phys. Lett. 73 $1980$ 393-398