非対称行列向けマルチシフト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 } 逆変換 Aの固有ベクトル {vi } vi = Q ui 逆変換各部分の実行時間 • 全固有値を求める場合の演算量 – ヘッセンベルグ化： (10/3) – QR法： 10n3 （経験値） n3 Execution time (min) 45 QR 40 Hessenberg reduction 35 30 25 – 演算量（時間）の大部分をQR法が占める。 – QR法の高速化が必要 20 15 10 5 O Origin2000 1PU （R12000，400MHz）上での実行時間 K. Braman et al: “The Multishift QR Algorithm I” より抜粋 Hessenberg 化の時間は推定値 LM 10 00 TU B1 00 0 RD B1 25 TO 0 LS 20 0 PD 0 E2 9 M HD 61 32 00 B 0 QR法の特徴と高性能計算高性能計算の条件 QR法の特徴 •計算の逐次性（bulge-chasing 型演算） •計算の並列性 •低いデータ再利用性 •高いデータ再利用性（キャッシュの有効利用）オリジナルの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 H0 AkH0 0 0 バルジ 0 • (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホルダー変換H0 を求める • Ak’ = H0tAk H0 （バルジ導入） 0 • 直交行列による相似変換を繰り返すことにより， Ak’ を再びヘッセンベルグ行列に変形する（バルジ追跡） • これにより得られる行列を Ak+2 とする Ak+2 0 陰的ダブルシフトQR法の演算パターン • バルジ追跡における演算 – 3×3のハウスホルダー変換を左右からかけることにより，bulge を1つ右下に動かす左から Hl を乗算 0 0にしたい要素右から Hl を乗算 0 0 更新される要素 • 演算の特徴 – 並列粒度は O(n) • 並列アルゴリズム： R. Suda et al.（1999），G. Henry et al. (2002) – QR法の1ステップで，各行列要素は3回のみ更新データ再利用性が低く，キャッシュの有効利用が困難２. マルチシフト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） m = 6 の場合 – シフト s1, s2, …, sm を含む (m+1)×(m+1) のバルジを導入 – (m+1)×(m+1) のハウスホルダー変換を用いてこれを追跡 0 • 方式２：小さなバルジ複数個を追跡（Braman et al., 2003） – シフトを2つずつ組にし， 3×3 の小さなバルジ m/2 個を導入 – ダブルシフトQR法と同様にして，これらを追跡 – バルジ同士は3行離れていれば，干渉なしに計算が可能 • 密に詰めることで，データ参照の局所性を向上 – 得られる行列は，無限精度演算では方式１と同じ 0 方式１と方式２の収束性比較 • シフトの数と総演算量との関係 – DHSEQR：方式１（LAPACK） – TTQR：方式２ K. Braman et al.: “The Multishift QR Algorithm I” (2002) より引用有限精度演算では，方式２が有利（Cf. Watkins による理論的解析）以下では方式２について考えるレベル3 BLAS の利用 • 更新処理の分割 – m/2 個のバルジをそれぞれ r 行追跡する際，まず対角ブロックのみを更新 – 非対角ブロックは，更新に使ったハウスホルダー変換を1個の行列に蓄積し，後でまとめて更新非対角ブロックの更新でレベル3 BLAS が使用可能 0 バルジ（3×3）最初に更新まとめてBLAS3で更新 • r の決め方 – 演算量の点から，r ～ 3m とするのが最適 – このとき，演算量は方式１の約2倍（非ゼロ構造を利用した場合）蓄積されたハウスホルダー変換の非ゼロ構造３. 本研究の目的 • シフト数 m の最適化 – m を大きくすると，BLAS3 部分の性能は向上 – しかし，シフトの計算量は増加（O(m3)） – 最適な m の値は計算機と問題サイズに依存 • r （1回のバルジ追跡行数）の最適化 – 演算量最小化のためには r ～ 3m が最適 – BLAS3 の実行時間は必ずしも演算量に比例しないため，実行時間の最適値はこれとは異なる n m の最適値 1000～1999 60 2000～2499 116 2500～3999 150 4000～ 180 実験的に求めた m の最適値（Origin2000上，Braman et al. より引用） m，r に対する自動チューニングの必要性本発表では，性能予測モデルに基づくシフト数 m の自動最適化について報告する４. 性能予測手法 • 階層的な性能モデリング（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) 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 の双一次式を用いる fDGEMM_NN (m, n, k) = fDGEMM_NNn (m, k) = (a11n m + a10n) k + (a01n m + a00n) – 係数 a00n ，a01n ，a10n ，a11n は実測データから最小二乗法で決定 – 代表値以外の n に対する値は，一次補間により計算基本ルーチンの性能モデリング（続き） • DGEMM（’N’,’T’） – ‘N’,’N’の場合と同様にして実行時間の予測関数を構成 • BCHASE2，COPY1，COPY2 – サイズを決めるパラメータは m, r の2つ – r の代表的な値に対し，実行時間を m の多項式として近似 • BCHASE1，BCHASE3，HQR – サイズを決めるパラメータは m のみ – 実行時間を m の多項式（3次式）として近似マルチシフトQR法全体の性能モデリング • 基本的な考え方 – 基本演算の各ルーチンと同じ引数を持ち，演算を行う代わりに実行時間の予測値のみを計算して返すルーチンを作成 – マルチシフトQR法中の基本演算ルーチンをこれらのルーチンで置き換えることにより，帯行列化の実行時間を予測するプログラムを作成計算プログラム予測プログラム 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 • 本方式のメリット – 複雑な解析的モデルの構築が不要 – 予測に必要な時間は O(N2/m2) ５. 実験結果 • 評価環境 – 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に規格化）相対実行時間相対実行時間 1.6 1.6 1.4 1.4 1.2 1.2 1 m=30 m=60 m=90 m=120 0.8 0.6 1 0.6 0.4 0.4 0.2 0.2 0 1000 2000 4000 予測結果 8000 m=30 m=60 m=90 m=120 0.8 0 1000 2000 4000 8000 実測結果 • モデルは m を変えたときの相対実行時間の変化を定性的に再現 → シフト数の自動最適化に使える可能性あり • ただし，絶対時間では20～40%の誤差 PowerPC G5上での予測結果（1CPU，続き） • m を変えたときの相対実行時間（実測時間の最小値を1に規格化）相対実行時間相対実行時間 1.8 1.8 1.6 1.6 1.4 1.4 1.2 1.2 m=30 m=60 m=90 m=120 1 0.8 0.6 0.8 0.6 0.4 0.4 0.2 0.2 0 1000 2000 4000 予測結果 8000 m=30 m=60 m=90 m=120 1 0 1000 2000 4000 実測結果 • 予測時間は，N が小さいとき過大評価，大きいとき過小評価 8000 PowerPC G5上での予測結果（2CPU） • m を変えたときの相対実行時間（最小値を1に規格化）相対実行時間相対実行時間 1.8 1.6 1.6 1.4 1.4 1.2 1.2 m=30 m=60 m=90 m=120 1 0.8 0.6 m=30 m=60 m=90 m=120 0.8 0.6 0.4 0.4 0.2 0.2 0 1 1000 2000 4000 予測結果 8000 0 1000 2000 4000 8000 実測結果 • モデルは m を変えたときの相対実行時間の変化を定性的に再現 • 絶対時間の誤差は1CPUの場合と同程度 Opteron 上での予測結果（4CPU） • m を変えたときの相対実行時間（最小値を1に規格化）相対実行時間相対実行時間 1.6 1.6 1.4 1.4 1.2 1.2 1 m=30 m=60 m=90 m=120 0.8 0.6 1 0.8 0.6 0.4 0.4 0.2 0.2 0 1000 2000 4000 予測結果 8000 m=30 m=60 m=90 m=120 0 1000 2000 4000 8000 実測結果 • モデルは 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で更新次の対角ブロックの更新で必要