分散メモリ型スパースソルバの開発と評価

応用数理工学特論
線形計算と
ハイパフォーマンスコンピューティング
第8回
名古屋大学計算理工学専攻
山本有作
１. はじめに
• 本研究で対象とする問題
– 標準固有値問題 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法の基礎
• アルゴリズム（Francis, 1961 & Kublanovskaya, 1961）
– 行列 A0 から始めて次のように QR 分解と相似変換を繰り返す。
• A0 = Q1R1
• A1 = R1Q1 (= Q1–1A0 Q1 )
• A1 = Q2R2
• A2 = R2Q2 (= Q2–1A1 Q2 = Q2–1Q1–1A0 Q1Q2 )
• 収束定理
– 適当な条件の下で，Ak は（ブロック）上三角行列に収束
– A0 の固有値を絶対値の大きい順に l1, l2, … , ln とすると，対角要素 aij
は li に１次収束
– 非対角要素 aij （j < i）は収束率 rij ≡ |li| / |lj| で 0 に１次収束
ダブルシフトQR法
• シフトの導入
– Ak から固有値 li の近似値 s を引いた行列に対して QR 法を適用
• Ak – s I = Qk Rk
• Ak+1 = Rk Qk + s I (= Qk–1Ak Qk )
rij ≡ |li – s| / |lj – s| ≒ 0 より収束が加速
• ダブルシフトQR法
– 共役複素数の固有値ペアに対し，シフト s, s による2反復をまとめて行う
• (Ak – s I)(Ak – s I) = Qk Rk
• Ak+2 = Qk–1Ak Qk
– シフトは右下隅の 2×2 行列の固有値に取る
→ 局所的に2次収束
– 複素数の固有値を持つ場合でも，実数演算だけで QR 法を実行可能
Hessenberg 形と Implicit Q 定理
• Hessenberg 形の利用
– Ak が Hessenberg 形のとき，
• QR法の1ステップは O(n2) で実行可能
• Ak+1 も Hessenberg 形
• Implicit Q 定理の利用により，QR 分解を陽に行う
ことなく1ステップの計算を実行可能
0
Hessenberg 形
aij = 0 if i > j+1
A0 を Hessenberg 形に相似変換してからQR法を適用
• Implicit Q 定理
– U，V が直交行列で，G = UtAU， H = VtAV が共に既約な Hessenberg 行
列であるとする。このとき，もし U と V の第1列が等しいならば，±1 を
要素に持つ対角行列 D が存在して，V = UD，H = DGD が成り立つ
– すなわち，行列 A を既約な Hessenberg 行列に相似変換する直交行列
は，第1列目だけが与えられれば（実質的に）一意に定まる
Implicit shift QR法
• 1ステップの計算
• (Ak – s2I)(Ak – s1I) の第1列を e1の定数倍にするハウスホル
ダー変換H0 を求める
• Ak’ = H0tAk H0
• 直交行列による相似変換を繰り返すことにより， Ak’ を再び
Hessenberg 行列に変形する（bulge-chasing）
Ak
H0 AkH0
0
0
0
• この変換の性質
– ある直交行列 H による相似変換
– H の第1列は，Qk の第1列に等しい
• 上記第3のステップは第2行・第2列以降のみに影響
– H は Ak を Hessenberg 形に変換
Implicit Q 定理より，Ak+1 = H tAk H はQR法の1ステップと等価
0
Ak+1
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倍
（非ゼロ構造を利用した場合）
蓄積されたハウス
ホルダー変換の非
ゼロ構造
マルチシフトQR法の性能
• 流体力学の問題に対する計算時間の例
– PowerPC G5 2.0GHz （1CPU）
– 行列は金田研究室水野氏提供
1600
QR法
45000
QR法
1400
Hessenberg 化
40000
Hessenberg 化
行列生成
35000
行列生成
1200
30000
1000
25000
800
20000
600
15000
400
10000
200
0
5000
8 (DHSEQR)
8 (HQR_M)
N=3645
0
11 (DHSEQR)
11 (HQR_M)
N=11232
３. シフト数の最適化
• シフト数 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 の自動最適化を行うこ
とを考える
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上での予測結果（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 の場合と同程度

Download Report