Transcript 非対称行列向けマルチシフト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 }
逆変換
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で更新
次の対角ブロックの
更新で必要