キャッシュマシン向け三重対角化 アルゴリズムの性能予測方式 名古屋大学　計算理工学専攻 山本有作 日立製作所の山本有作です。 「～」について発表いたします。

目次 １. はじめに ２. 三重対角化アルゴリズムとキャッシュ向け最適化 ３. 本研究の目的 ４. 性能予測手法 ５. 実験結果 １.　はじめに ２.　三重対角化アルゴリズムとキャッシュ向け最適化 ３.　本研究の目的 ４.　性能予測手法 ５.　実験結果 ６.　関連研究 ７.　まとめと今後の課題 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

１. はじめに 本研究で対象とする問題 応用分野 標準固有値問題 Ax = λx A ： 実対称またはエルミートの密行列 １.　はじめに 本研究で対象とする問題 標準固有値問題　Ax = λx A ： 実対称またはエルミートの密行列 応用分野 分子計算，統計計算など 分子計算では，10万元以上の超大規模行列の固有値計算も需要あり ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

固有値計算の流れ Q*AQ = T (Q: 直交行列) Tyi =λi yi xi = Qyi 実対称行列 A 計算内容 計算手法 三重対角化 Q*AQ = T (Q: 直交行列) ハウスホルダー法 三重対角行列 T QR法 分割統治法 二分法・逆反復法 MR3アルゴリズム 三重対角行列の 固有値・固有ﾍﾞｸﾄﾙ計算 Tyi =λi yi Tの固有値 {λi }， 　 固有ベクトル {yi } 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 xi = Qyi 逆変換 逆変換 Aの固有ベクトル {xi }

固有値計算の処理時間の内訳 （P4クラスタ，16ノード） Paolo Bientinesi et al.: A Parallel Eigensolver for Dense Symmetric Matrices using Multiple Relatively Robust Representations (2004) より 　　三重対角行列の固有値・固有ベクトル計算に最新のアルゴリズム（MR3）を使った場合，処理時間の大部分を三重対角化が占める。 　　　　　ハウスホルダー法の高速化が必要

２. 三重対角化アルゴリズムとキャッシュ向け最適化 ２.　三重対角化アルゴリズムとキャッシュ向け最適化 ハウスホルダー法による三重対角化 鏡像変換 H = I – a uut による列の消去 Hは対称な直交行列 与えられたベクトルの第１成分以外をゼロにする。 第 k ステップでの処理 ベクトル 左からH を乗算 左からH を乗算 右からH を乗算 k 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

ハウスホルダー法のアルゴリズム [Step 1] k = 1からN –2まで以下の[Step 2] ～ [Step 8]を繰り返す。 　[Step 2]　内積 σ(k) = sqrt(d (k)t d (k)) の計算 　[Step 3]　鏡像変換ベクトル作成： 　　u(k) = (d (k)1 – sgn(d (k)1)σ(k), d (k)2, … , d (k)N – k ) 　[Step 4]　規格化定数の計算：　α(k) = 2 /∥u(k)∥2 　[Step 5]　行列ベクトル積：　p(k) =α(k)C (k)u(k) 　[Step 6]　ベクトルの内積：　β(k) =α(k) u(k)tp(k) /2 　[Step 7]　ベクトルの加算：　q(k) = p(k) –β(k)u(k) 　[Step 8]　行列のrank-2更新：　C (k) = C (k) – u(k)q(k)t – q(k)u(k)t ピボット列 ピボット行 d (k)1 C (k) d (k) k

ハウスホルダー法の特徴と問題点 特徴 問題点 全演算量のほとんどが，level-2 BLAS である行列ベクトル積と行列のrank-2更新で占められる。 全演算量： 約 (4/3)N3 行列ベクトル積の演算量： 約 (2/3)N3 rank-2更新の演算量： 約 (2/3)N3 問題点 level-2 BLASのため，行列データの再利用性なし キャッシュミス，SMPでのメモリ競合の影響大 最近のマイクロプロセッサ上での高性能が期待できない。

キャッシュマシン向けの三重対角化アルゴリズム 原理（Bishof et al., 1993, 1994） 密行列 A をまず半帯幅Lの帯行列 B に変換し，次にこの帯行列を三重対角行列 T に変換する。 三重対角化を２段階で行うことの利点 帯行列への変換は， level-3 BLAS のみを使って実行可能。 帯行列から三重対角行列への変換の演算量はO(N2L)であり， 　　前半部に比べてずっと小さい。 次数 N 半帯幅 L 帯行列化 村田法 約 (4/3)N3 O(N2L) A B T

帯行列化のアルゴリズム ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 ブロックベクトル ブロック鏡像変換によるブロック列の消去 ブロック鏡像変換 H = I – UαUt Hは直交行列 与えられたブロックベクトルの 第１ブロックを上三角行列にし， それ以外をゼロにする。 第kステップでの処理 左からH を乗算 右からH を乗算 左からH を乗算 非ゼロ要素 ゼロにしたい部分 影響を受ける部分

帯行列化のアルゴリズム（続き） 本アルゴリズムの特徴 [Step 1]　K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 [Step 2]　 D(K) を第１ブロックが上三角行列で第２ブロック以下が ゼロ行列であるブロックベクトルに変換するブロック鏡像変換 I – U(K)α(K)U(K)tを求める。 [Step 3]　行列・ブロックベクトル積：　P(K) = C (K)U(K)α(K) [Step 4]　ブロックベクトルの内積：　β(K) = α(K)tU(K)tP(K) / 2 [Step 5]　ブロックベクトルの加算：　Q(K) = P(K) –U(K)β(K) [Step 6]　行列のrank-2L更新：　C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t すべて level-3 BLAS（行列乗算） 本アルゴリズムの特徴 演算が level-3 BLAS 中心のため，キャッシュの有効利用が可能 SMPにおけるメモリ競合の影響を低減可能

多段化による性能向上 原理（Wu et al., 1996） 帯行列化のアルゴリズムにおいて，複数のピボット列・行を溜めておき，後でまとめて更新 M 本溜めることで，rank-2L 更新の部分を rank-2LM 更新にでき，データ再利用性を更に向上可能

Opteron (1.6GHz) 上での三重対角化の性能 L=24, M=4 L=48 Performance (GFLOPS) M=32 Matrix size 本アルゴリズムは他の方法に比べて非常に高性能 N = 3840 のとき，ピークの50%以上の性能を達成

本アルゴリズムの課題 性能を大きく左右するパラメータが２つ存在 各パラメータに関するトレードオフ 帯行列の半帯幅 L 多段化の段数 M L が大　 →　BLAS性能は上がるが，村田法の演算量が増大 M が大　→　BLAS性能は上がるが，帯行列化の演算量が増大 　　対象とするプロセッサおよび行列サイズ N に応じて，L と M を最適化することが必要

３. 本研究の目的 キャッシュマシン向けの三重対角化アルゴリズムに対し，高精度な性能予測モデルを構築する ３.　本研究の目的 キャッシュマシン向けの三重対角化アルゴリズムに対し，高精度な性能予測モデルを構築する 　　　　　　最適な性能パラメータ L，M の値を実行前に推定する 　　　　　　　自動チューニング型ライブラリのための基礎技術 性能予測モデルの他の応用 グリッド上での最適マシン選択 グリッド上でのスケジューリング

４. 性能予測手法 階層的な性能モデリング（Cuenca et al., 2004） ４.　性能予測手法 階層的な性能モデリング（Cuenca et al., 2004） 帯行列化のアルゴリズムは，多種類のlevel-3 BLASへのコールから成る level-3 BLASの各ルーチンの性能を精度良くモデル化できれば，その積み上げによりアルゴリズム全体の性能も精度良くモデル化できるはず

帯行列化で用いるlevel-3 BLASルーチン アルゴリズム 　[Step 1]　K = 1からN /L–1まで以下の[Step 2] ～ [Step 6]を繰り返す。 　[Step 2]　 ブロック鏡像変換 I – U(K)α(K)U(K)t の計算 BlockHouse 　[Step 3]　行列・ブロックベクトル積：　P(K) = C (K)U(K)α(K) DSYMM 　[Step 4]　ブロックベクトルの内積：　β(K) = α(K)tU(K)tP(K) / 2 DGEMM 　[Step 5]　ブロックベクトルの加算：　Q(K) = P(K) –U(K)β(K) DGEMM 　[Step 6]　行列のrank-2L更新：　C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t　　 DSYR2K 必要なBLASルーチン DGEMM　3種 （‘N’,’N’， ‘N’,’T’， ‘T’,’N’） DSYMM DSYR2K BlockHouse （BLASではないがブロック鏡像変換を求めるルーチン）

BLAS性能のモデリング DGEMM（’N’,’N’）の場合 fDGEMM_NNの構成方法 機能：　C := aC + bAB の計算　（A: m×k， B: k×n， C: m×n） 実行時間の予測関数：　fDGEMM_NN(m, n, k) fDGEMM_NNの構成方法 m, n, k の全範囲で f を１つの多項式で近似すると，誤差が大きい n の代表的な値に対し，f を m, k の多項式により近似 多項式としては， m, k の双一次式を用いる 係数 a00n ，a01n ，a10n ，a11n は実測データから最小二乗法で決定 代表値以外の n に対する値は，一次補間により計算 fDGEMM_NN (m, n, k) = fDGEMM_NNn (m, k) 　　　　　　　　　　　 = (a11n m + a10n) k + (a01n m + a00n)

BLAS性能のモデリング（続き） DGEMM（’N’,’T’） および DGEMM（’T’,’N’） DSYMM DSYR2K ‘N’,’N’の場合と同様にして 実行時間の予測関数を構成 DSYMM サイズを決めるパラメータは n, m の2つのみ m の代表的な値に対し，実行時間を n の一次式として近似 DSYR2K DSYMM の場合と同様 BlockHouse サイズを決めるパラメータは，ベクトル長 n とブロック幅 L の2つ level-3 BLAS の場合と異なり，一次式／双一次式近似では実行時間がうまく再現できない (n, L)平面上の格子点で実行時間を測定し，補間により実行時間を予測

帯行列化アルゴリズムの性能モデリング 基本的な考え方 本方式のメリット BLASの各ルーチンと同じ引数を持ち，演算を行う代わりに実行時間の予測値のみを計算して返すルーチンを作成 帯行列化プログラム中のBLASをこれらのルーチンで置き換えることにより，帯行列化の実行時間を予測するプログラムを作成 本方式のメリット 複雑な解析的モデルの構築が不要 予測に必要な時間は O(N/L) 計算プログラム 予測プログラム DO K=1, N/L CALL BlockHouse(m,L,...) CALL DSYMM(m,n,...) CALL DGEMM(m,n,k,...) CALL DSYR2K(m, k,...) END DO DO K=1, N/L T1=T1+BlockHouse_TIME(m,L,...) T1=T1+DSYMM_TIME(m,n,...) T1=T1+DGEMM_TIME(m,n,k,...) T1=T1+DSYR2K_TIME(m, k,...) END DO

５. 数値実験 評価項目 評価環境 評価例題 BLAS性能の予測 帯行列化アルゴリズムの性能の予測 パラメータ L，M の最適化 ５.　数値実験 評価項目 BLAS性能の予測 帯行列化アルゴリズムの性能の予測 パラメータ L，M の最適化 評価環境 Opteron（1.6GHz），g77，GOTO BLAS Alpha 21264A（750MHz），g77，GOTO BLAS Power PC G5（2.0GHz），g77，ATLAS 3.7.0 評価例題 N=480 ～ 9600 の対称密行列の三重対角化

DGEMM(‘N’,’N’), Opteron, n=12 BLAS性能の予測 前述の方式で各 BLAS について実行時間の予測関数を構成し，各プロセッサ上で様々な行列サイズに対して実測データと比較 すべてのサイズパラメータ（m, n, k など）が6を超える領域では，どのプロセッサ／BLASでも予測誤差は5～10%程度と高精度 予測時間／実測時間 予測時間／実測時間 k m k m DGEMM(‘N’,’N’), Opteron, n=12 DGEMM(‘N’,’N’), G5, n=12

帯行列化アルゴリズムの性能の予測 L，M を変えたときの予測時間と実測時間（Opteron, N=1920） 予測誤差は10%以内 時間（s） 時間（s） M M L L 予測時間 実測時間

帯行列化アルゴリズムの性能の予測（続き） L，M を変えたときの予測時間と実測時間（Alpha 21264A, N=3840） 予測誤差は5%以内 モデルは L，M を変えたときの実行時間の変化をよく再現 時間（s） 時間（s） M M L L 予測時間 実測時間

帯行列化アルゴリズムの性能の予測（続き） L，M を変えたときの予測時間と実測時間（PowerPC G5, N=1920） 予測誤差は10%以内 モデルは L，M を変えたときの実行時間の変化をよく再現 時間（s） 時間（s） M M L L 予測時間 実測時間

パラメータ L，M の最適化 N が与えられたとき，（L，M）の各組について，三重対角化の実行時間を次の2つの時間の和として予測 帯行列化の実行時間：　モデルにより予測 村田法の実行時間：　　 実測値 　　　　　　　　　　この予測値を最小とする（L，M）を決定

パラメータ L，M の最適化（続き） 予測した（L，M）による実行結果（G5，N=960～9600） 　　　　　　　　本モデルを用いたパラメータの事前最適化は有効 性能（MFLOPS） 行列サイズ N

６. 関連研究 三重対角化プログラムの自動チューニング（Katagiri et al., 2000） ６.　関連研究 三重対角化プログラムの自動チューニング（Katagiri et al., 2000） 分散メモリ向けの三重対角化プログラムにおいて，ループ展開の段数，通信関数の種類などのパラメータを自動的に最適化する方式 パラメータを変化させてプログラム全体を何回も実行することにより最適値を求めるため，最適化に時間がかかるという問題点がある 階層的な性能モデリング（Cuenca et al., 2004，Dackland et al., 1996） 線形計算プログラムの自然な階層構造を利用して，下位のルーチン（BLASなど）の性能モデルをまず構築し，それを用いて順次上位のルーチンの性能モデルを構築していく方式 本研究のモデルもこの考え方に基づく ただし，従来の適用例は，LU分解やQR分解などの基本的分解，およびヤコビ法などの単純なアルゴリズムに限られている

７. まとめと今後の課題 本研究のまとめ 今後の課題 ７.　まとめと今後の課題 本研究のまとめ キャッシュマシン向けの三重対角化アルゴリズムに対し，階層的なモデリング手法に基づく性能予測モデルを開発した Opteron，Alpha 21264A，Power PC G5上で評価した結果，予測誤差は5～10%であり，L，M を変えたときの実行時間の変化も良く再現できた 本モデルを用いて L，M の最適値を推定し，三重対角化を実行したところ，真の最適値での性能とほぼ等しい性能が得られた 今後の課題 より多くの種類のプロセッサ上での検証 共有メモリ／分散メモリ型並列計算機向けの拡張 性能安定化手法の適用 自動チューニング型ライブラリへの適用

謝辞 日頃から有益な議論をして頂いている自動チューニング研究会のメンバーに感謝いたします