非対称行列向けマルチシフトQR法の 性能予測方式 名古屋大学　計算理工学専攻 山本有作 2006年6月12日 HPC研究会 日立製作所の山本有作です。 「～」について発表いたします。

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

１. はじめに 本研究で対象とする問題 応用分野 標準固有値問題 Ax = lx A ： n×n 非対称密行列 MHD 化学工学 量子化学 １.　はじめに 本研究で対象とする問題 標準固有値問題　Ax = lx A ： n×n 非対称密行列 応用分野 MHD 化学工学 量子化学 流体力学 Cf. Bai and Demmel: A test matrix collection for non-Hermitian eigenvalue problems. ここでは特に，Aが実対称またはエルミートの密行列の場合を考える。

非対称行列の固有値計算の流れ 密行列 A 計算内容 計算手法 QTAQ = H (Q: 直交行列) ハウスホルダー法 固有値の計算 QR法 ヘッセンベルグ化 ヘッセンベルグ行列 H 固有値の計算 高い信頼性 QR法 分割統治法 A の固有値 {li } 固有ベクトルの計算 Hui =li ui 密行列Aをまず三重対角行列Tに相似変換してからTの固有値・固有ベクトルを求めるのが最も一般的な計算法。 三重対角化には，後に述べるハウスホルダー法を使う場合がほとんど。 三重対角行列の固有値・固有ベクトルの計算には，色々なアルゴリズムがある。 Hの固有ﾍﾞｸﾄﾙ {ui } 逆変換 vi = Q ui 逆変換 Aの固有ベクトル {vi }

各部分の実行時間 全固有値を求める場合の演算量 ヘッセンベルグ 化： (10/3) n3 QR法： 10n3 （経験値） Execution time (min) Origin2000 1PU （R12000，400MHz）上での実行時間 K. Braman et al: “The Multishift QR Algorithm I” より抜粋 Hessenberg 化の時間は推定値

QR法の特徴と高性能計算 QR法の特徴 高性能計算の条件 計算の逐次性 （bulge-chasing 型演算） 計算の並列性 低いデータ再利用性 高いデータ再利用性 　　（キャッシュの有効利用） オリジナルのQR法のままでは高性能化が困難 高性能計算の条件を満たす新しいアルゴリズムが必要 マルチシフトQR法

ダブルシフトQR法 原理 陰的ダブルシフトQR法 Ak の右下の 2×2 行列の固有値 s1，s2 をシフトとして用い，QR法の2ステップを一度に実行 (Ak – s1 I)(Ak – s2 I) = Qk Rk Ak+2 = Qk–1Ak Qk 陰的ダブルシフトQR法 (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホルダー変換H0 を求める Ak’ = H0tAk H0 　（バルジ導入） 直交行列による相似変換を繰り返すことにより， Ak’ を再びヘッセンベルグ行列に変形する（バルジ追跡） これにより得られる行列を Ak+2 とする Ak H0 AkH0 バルジ Ak+2

陰的ダブルシフトQR法の演算パターン バルジ追跡における演算 演算の特徴 3×3のハウスホルダー変換を左右からかけることにより，bulge を1つ右下に動かす 演算の特徴 並列粒度は O(n) 並列アルゴリズム：　R. Suda et al.（1999），G. Henry et al. (2002) QR法の1ステップで，各行列要素は3回のみ更新 　　データ再利用性が低く，キャッシュの有効利用が困難 左から Hl を乗算 右から Hl を乗算 0にしたい要素 更新される要素

２. マルチシフトQR法 原理 （Bai & Demmel, 1989） シフト数増加の効果 Ak の右下の m×m 行列の固有値 s1，s2 , … , sm をシフトとして用い，QR法の m ステップを一度に実行 (Ak – smI) ・・・ (Ak – s2I)(Ak – s1I) = Qk Rk Ak+m = Qk–1Ak Qk シフト数増加の効果 並列粒度は O(m) 倍に増加 行列の各要素に対する更新回数も O(m) 倍に増加 データ再利用性の向上により，キャッシュの有効利用が可能 　　（WY representation を用いた BLAS3化）

マルチシフトQR法におけるバルジ追跡 方式１： 大きなバルジ1個を追跡 （Bai & Demmel, 1989） シフト s1, s2, …, sm を含む (m+1)×(m+1) のバルジを導入 (m+1)×(m+1) のハウスホルダー変換を用いてこれを追跡 方式２： 小さなバルジ複数個を追跡 （Braman et al., 2003） シフトを2つずつ組にし， 3×3 の小さなバルジ m/2 個を導入 ダブルシフトQR法と同様にして，これらを追跡 バルジ同士は3行離れていれば，干渉なしに計算が可能 密に詰めることで，データ参照の局所性を向上 得られる行列は，無限精度演算では方式１と同じ m = 6 の場合

K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 方式１と方式２の収束性比較 シフトの数と総演算量との関係 DHSEQR： 方式１（LAPACK） TTQR： 方式２ K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用 有限精度演算では，方式２が有利 （Cf. Watkins による理論的解析） 以下では方式２について考える

レベル3 BLAS の利用 更新処理の分割 r の決め方 m/2 個のバルジをそれぞれ r 行追跡する際， まず対角ブロックのみを更新 　 まず対角ブロックのみを更新 非対角ブロックは，更新に使ったハウスホルダー変換を1個の行列に蓄積し，後でまとめて更新 　　非対角ブロックの更新で レベル3 BLAS が使用可能 r の決め方 演算量の点から，r ～ 3m とするのが最適 このとき，演算量は方式１の約2倍 　　（非ゼロ構造を利用した場合） バルジ（3×3） 最初に更新 まとめてBLAS3で更新 蓄積されたハウスホルダー変換の非ゼロ構造

３. 本研究の目的 シフト数 m の最適化 r （1回のバルジ追跡行数）の最適化 m，r に対する自動チューニングの必要性 ３.　本研究の目的 シフト数 m の最適化 m を大きくすると，BLAS3 部分の性能は向上 しかし，シフトの計算量は増加（O(m3)） 最適な m の値は計算機と問題サイズに依存 r （1回のバルジ追跡行数）の最適化 演算量最小化のためには r ～ 3m が最適 BLAS3 の実行時間は必ずしも演算量に比例しないため，実行時間の最適値はこれとは異なる 　　　　m，r に対する自動チューニングの必要性 n m の最適値 1000～1999 60 2000～2499 116 2500～3999 150 4000～ 180 実験的に求めた m の最適値 （Origin2000上，Braman et al. より引用） 本発表では，性能予測モデルに基づくシフト数 m の自動最適化について報告する

４. 性能予測手法 階層的な性能モデリング（Cuenca et al., 2004） 反復回数の推定 ４.　性能予測手法 階層的な性能モデリング（Cuenca et al., 2004） マルチシフトQR法のアルゴリズムは，レベル3 BLAS，シフト計算のためのQR法など，数種類の基本演算ルーチンから構成される 各ルーチンの性能を精度良くモデル化できれば，その積み上げによりアルゴリズム全体の性能も精度良くモデル化できるはず 反復回数の推定 QR法は反復法のため，実行時間予測には反復回数の推定が必要 固有値１個当たり，平均4反復で収束すると仮定

マルチシフトQR法を構成する基本演算ルーチン 8種の基本演算ルーチン HQR：　シフト計算のためのQR法 BCHASE1：　m/2個のバルジの導入 BCHASE2：　対角ブロック内でのバルジ追跡 BCHASE3：　バルジの追い出し DGEMM(‘N’, ‘N’)：　非対角行ブロックの更新 DGEMM(‘N’, ‘T’)：　非対角列ブロックの更新 COPY1：　コピー (1) COPY2：　コピー (2) バルジ（3×3） 最初に更新 まとめてBLAS3で更新

基本演算ルーチンの性能モデリング 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)

基本ルーチンの性能モデリング（続き） DGEMM（’N’,’T’） BCHASE2，COPY1，COPY2 ‘N’,’N’の場合と同様にして 実行時間の予測関数を構成 BCHASE2，COPY1，COPY2 サイズを決めるパラメータは m, r の2つ r の代表的な値に対し，実行時間を m の多項式として近似 BCHASE1，BCHASE3，HQR サイズを決めるパラメータは m のみ 実行時間を m の多項式（3次式）として近似

マルチシフトQR法全体の性能モデリング 基本的な考え方 本方式のメリット 基本演算の各ルーチンと同じ引数を持ち，演算を行う代わりに実行時間の予測値のみを計算して返すルーチンを作成 マルチシフトQR法中の基本演算ルーチンをこれらのルーチンで置き換えることにより，帯行列化の実行時間を予測するプログラムを作成 本方式のメリット 複雑な解析的モデルの構築が不要 予測に必要な時間は O(N2/m2) 計算プログラム 予測プログラム DO K=1, N/L CALL HQR(m,L,...) CALL BCHASE1(m,n,...) CALL BCHASE2(m,n,k,...) CALL DGEMM(m,n,k,...) CALL BCHASE3(m, k,...) END DO DO K=1, N/L T1=T1+HQR_TIME(m,L,...) T1=T1+BCHASE1_TIME(m,n,...) T1=T1+BCHASE2_TIME(m,n,k,...) T1=T1+DGEMM_TIME(m,n,k,...) T1=T1+BCHASE3_TIME(m, k,...) END DO

５. 実験結果 評価環境 評価例題 Power PC G5（2.0GHz） 2way，IBM XL Fortran，GOTO BLAS ５.　実験結果 評価環境 Power PC G5（2.0GHz） 2way，IBM XL Fortran，GOTO BLAS Opteron（1.6GHz） 4way，g77，GOTO BLAS 評価例題 N=1000 ～ 8000 の行列の固有値計算 m = 30，60，90，120 の４通りのシフト数 入力行列は， [0,1] の乱数行列をハウスホルダー法でヘッセンベルグ化して使用

PowerPC G5上での予測結果（1CPU） m を変えたときの相対実行時間（最小値を1に規格化） モデルは m を変えたときの相対実行時間の変化を定性的に再現 　　　→　シフト数の自動最適化に使える可能性あり ただし，絶対時間では20～40%の誤差 相対実行時間 相対実行時間 予測結果 実測結果

PowerPC G5上での予測結果（1CPU，続き） m を変えたときの相対実行時間（実測時間の最小値を1に規格化） 予測時間は，N が小さいとき過大評価，大きいとき過小評価 相対実行時間 相対実行時間 予測結果 実測結果

PowerPC G5上での予測結果（2CPU） m を変えたときの相対実行時間（最小値を1に規格化） モデルは m を変えたときの相対実行時間の変化を定性的に再現 絶対時間の誤差は1CPUの場合と同程度 相対実行時間 相対実行時間 予測結果 実測結果

Opteron 上での予測結果（4CPU） m を変えたときの相対実行時間（最小値を1に規格化） 絶対時間の誤差は PowerPC G5 の場合と同程度 相対実行時間 相対実行時間 予測結果 実測結果

検討すべき課題 誤差の原因の究明 N による誤差の系統的な変化 その他の系統誤差 基本演算ルーチンごとの誤差の調査など，より詳細な解析が必要 平均反復回数の違いが影響している可能性あり 基本演算ルーチンのモデルの誤差 その他の系統誤差 基本演算ルーチンごとの誤差の調査など，より詳細な解析が必要

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

７. まとめと今後の課題 まとめ 今後の課題 マルチシフトQR法に対し，階層的なモデリング手法を用いて実行時間を予測するモデルを開発した ７.　まとめと今後の課題 まとめ マルチシフトQR法に対し，階層的なモデリング手法を用いて実行時間を予測するモデルを開発した 1000元から8000元の行列による評価では，モデルはシフト数を変えたときの相対実行時間の変化を定性的に再現できた しかし，絶対時間では誤差が大きく，原因の究明が必要 今後の課題 誤差の原因解明とモデルの精密化 基本演算ルーチンの性能モデリングの自動化 自動チューニング型ライブラリへの展開

共有メモリ向けの並列化手法 Level-3 BLAS 内部での並列化 より効率的な並列化 非対角ブロックの更新において，次の対角ブロックの更新に必要な部分のみを先に更新 　　並列化困難な対角ブロックの更新を，非対角ブロックの更新と並列に実行可能 次の対角ブロックの 更新で必要 Bulge（3×3） 最初に更新 まとめてBLAS3で更新