Transcript スライド - 計算基盤研究室

応用数理工学特論
線形計算と
ハイパフォーマンスコンピューティング
2008年6月19日
名古屋大学 計算理工学専攻
山本有作
期末発表について
• 概要
– 日時： 7/10 および 7/17
– １～３人のグループで１つの発表を行う。
– 論文を読んでまとめるか，あるいは講義で勉強したことをもとに
自分たちで数値実験を行った結果を発表
• グループの決定
– 本日の講義後に，グループを申告してください。
– 未決定の人は，今週末までに決めてメールで申告してください。
• 発表テーマの決定
– テーマの候補をホームページに掲載します。
– 来週の講義までに，自分たちのテーマを決めて報告してください。
前回の概要
３. 非対称行列の固有値計算のための前処理
・ ハウスホルダー法によるヘッセンベルグ化
４. ヘッセンベルグ行列の固有値計算
・ べき乗法
・ 直交化付き同時反復法
・ QR法
今回の内容
４. ヘッセンベルグ行列の固有値計算（続き）
・ シフト付きQR法
・ 減次（デフレーション）
・ QR法の演算量
・ ヘッセンベルグ行列に対するQR法
今回の内容
５. ヘッセンベルグ行列の固有値計算の高性能化
・ QR法の問題点
・ マルチシフトQR法
・ level-3 BLAS の利用
６. 三重対角行列の固有値計算の高性能化
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 に１次収束
各部分の実行時間
• 全固有値を求める場合の演算量
– ヘッセンベルグ 化： (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
レベル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