非対称行列向けマルチシフト QR 法の性能予測方式

非対称行列向けマルチシフトQR法の性能予測方式非対称行列向けマルチシフトQR法の性能予測方式名古屋大学　計算理工学専攻山本有作 2006年6月12日 HPC研究会

目次１.　はじめに２.　マルチシフトQR法３.　本研究の目的４.　性能予測手法５.　実験結果６.　関連研究７.　まとめと今後の課題

１.　はじめに • 本研究で対象とする問題 • 標準固有値問題　Ax = lx • A： n×n 非対称密行列 • 応用分野 • MHD • 化学工学 • 量子化学 • 流体力学 • Cf. Bai and Demmel: A test matrix collection for non-Hermitian eigenvalue problems.

高い信頼性 非対称行列の固有値計算の流れ計算内容計算手法密行列 A QTAQ = H (Q: 直交行列) ハウスホルダー法ヘッセンベルグ化ヘッセンベルグ行列H 固有値の計算 QR法分割統治法 A の固有値 {li} 固有ベクトルの計算 Hui =li ui Hの固有ﾍﾞｸﾄﾙ {ui } 逆変換逆変換 vi = Qui Aの固有ベクトル {vi }

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

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

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

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

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

マルチシフトQR法におけるバルジ追跡 m = 6 の場合 • 方式１：　大きなバルジ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行離れていれば，干渉なしに計算が可能 • 密に詰めることで，データ参照の局所性を向上 • 得られる行列は，無限精度演算では方式１と同じ 0 0

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

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

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

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

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

基本演算ルーチンの性能モデリング • DGEMM（’N’,’N’）の場合 • 機能：　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) = (a11nm + a10n) k +(a01nm + a00n)

基本ルーチンの性能モデリング（続き） • DGEMM（’N’,’T’） • ‘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 • 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に規格化） • モデルは mを変えたときの相対実行時間の変化を定性的に再現 • 絶対時間の誤差は PowerPC G5 の場合と同程度相対実行時間相対実行時間予測結果実測結果

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

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

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

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

非対称行列向けマルチシフト QR 法の性能予測方式