Transcript 名古屋大学 大学院工学研究科 宮田 考史

SWoPP 旭川 2007
2007 年 8 月 2 日
複素非対称行列向け固有値解法の
CSX600 による高速化
宮田 考史 1
1 名古屋大学
山本 有作 1
大学院工学研究科
中村 佳正 2
2 京都大学
大学院情報学研究科
はじめに
• 本研究で扱う問題
– 標準固有値問題
• A:
複素 非対称 密行列
• すべての固有値・固有ベクトルを計算
• 応用分野
– 電磁流体力学
– 構造計算
2
非対称行列の固有値計算
• 行列を扱いやすい形に相似変換
– 非対角成分の消去
高速化対象
対角成分 = 固有値
有限回演算
密行列
0
反復計算
ヘッセンベルグ行列
ハウスホルダー法
演算量 : （10/3）n3
0
右上三角行列
QR 法
演算量 : 10 n3 （経験値）
– 演算量の大部分を占める QR 法の高速化が必要
3
QR 法の高速化
• 在来研究
– マルチシフト QR 法 （K. Braman et al., 2002）
• 複数のシフトを導入
• 演算の大部分が行列乗算（Lv. 3 BLAS が利用可能）
– キャッシュの有効利用
– 並列性の向上
• 本研究
– 演算の大部分を占める行列乗算を高速化
• コプロセッサの利用
– 1 チップに多数の演算器を集積
– GRAPE-DR , CSX600
4
コプロセッサ利用に伴う問題点
12000
実行時間（秒）
10000
その他
行列乗算
8000
6000
• 乱数行列 （n = 6000）
• マルチシフト QR 法により
すべての固有値を計算
•
計算機環境
– Xeon 3.2 GHz, Memory 8 GB
– ClearSpeed advance board
• CSX600 × 2
4000
2000
0
CSX600 （無） / （有）
シフト数増加
行列乗算以外の
• シフト数増加に伴う実行時間の変化演算も高速化したい！
– 行列乗算
– 行列乗算以外
減少
増大 （ボトルネック）
5
発表順序
• はじめに
• 在来研究
– 非対称 固有値問題の数値解法
• ダブルシフト QR 法
• マルチシフト QR 法
• 本研究
– コプロセッサ向けの固有値解法
• マルチシフト QR 法に基づく提案法
• 数値実験
• まとめ
6
ダブルシフト QR 法
•
ダブルシフト QR 法
–
–
シフト s1, s2 : 右下 2×2 小行列の固有値
QR 法の 2 ステップを一度に実行
•
①
行列：
A ヘッセンベルグ
Q ユニタリ
R 右上三角
0
0
1 反復の演算パターン
の第 1 列を e1 の定数倍にする
ハウスホルダー変換
を求める
②
0
・
・
・
（bulge 導入）
③ 相似変換の反復 （chasing 演算） により，
をヘッセンベルグ行列に変形
0
7
マルチシフト QR 法
• マルチシフト QR 法
行列：
A ヘッセンベルグ
Q ユニタリ
R 右上三角
– シフト s1, …, sm : 右下 m×m 小行列の固有値
– QR 法の m ステップを一度に実行
• 1 反復の演算パターン
– （m / 2） 個の bulge 導入
– chasing 演算により，ヘッセンベルグ行列に変形
（ダブルシフト QR 法と同様）
– 更新処理の分割により，
Lv. 3 BLAS を用いたキャッシュの有効利用が可能
0
（例：シフト数 m = 4）
8
Lv. 3 BLAS の利用
シフト数 m 増加の影響
• 更新処理の分割
– （m/2） 個の bulge を k 行 chasing
• 最初に 対角ブロックを 逐次的に更新
– 更新に用いたハウスホルダー変換を
1 個の行列に累積
BLAS 3
• 最後に 非対角ブロックを まとめて更新
– 累積した行列 と 非対角ブロックの行列乗算
k
0
バルジ（3×3）
逐次的に更新
まとめて更新
（BLAS 3 利用）
• 演算量増加
– chasing 行数 k = （3m/2） のとき，増加は最小限（約 2 倍）
• 演算量最小化 ≠ 実行時間最小化
9
アルゴリズムの修正（1）
再帰型アルゴリズムへの変換
k/q
k
0
0
逐次的に更新
逐次的に更新
まとめて更新（BLAS 3 利用）
（例：再帰段数 d = 1）
(m / 2) / q 個の bulge を
k / q 行 chasing
まとめて更新（BLAS 3 利用）
10
アルゴリズムの修正（2）
• マルチシフト QR 法の反復
– 右下
小行列が分離 （
）
• 分離した小行列の固有値計算
0
– ダブルシフト QR 法を適用
– シフト数 m 増加により，小行列の次元増加
• 更新処理の分割
0
ブロック化
逐次的に更新
まとめて更新
（BLAS 3 利用）
– 対角ブロックの更新 （右上三角化するまで）
• 更新に用いたハウスホルダー変換を
1 個の行列に累積
– 非対角ブロックの更新 （最後に 1 度だけ）
• 累積した行列 と 非対角ブロックの行列乗算
演算量削減
11
パラメータ増加の影響
パラメータ
メリット
デメリット
シフト数 m
複数ステップをまとめて
実行可能
演算量増加
（シフト計算など）
chasing 行数 k
行列乗算の性能向上
（非対角ブロックの更新）
その他の演算増加
（対角ブロックの更新）
分割数 q
再帰段数 d
行列乗算導入
（対角ブロックの更新）
演算量増加
（対角ブロックの更新）
12
数値実験
13
数値実験
• テスト問題
– [0, 1] の成分を持つ
乱数行列
• ハウスホルダー法によりヘッセンベルグ化
– すべての固有値と固有ベクトルを計算
• 計算機環境
– Xeon 3.2 GHz , Memory 8 GB
– Fortran 77 倍精度演算
– ClearSpeed advance board
• CSX600 × 2
– BLAS 3（行列乗算）
• CSXL Library
• Intel Math Kernel Library
14
数値実験
• アルゴリズム
– 従来法（マルチシフト QR 法）
• 行列乗算
– 非対角ブロックの更新
CSX600
– 提案法（マルチシフト QR 法 + 再帰型・ブロック化）
• 行列乗算
– 非対角ブロックの更新
– 対角ブロックの更新
– 分離した小行列の固有値計算
– パラメータは最適に近い値
• シフト数 m , chasing 行数 k ： 従来法と提案法で共通の値
• 再帰段数 d = 1 , 分割数 q = m / 40
– 固有値ベクトル計算は共通 （性能評価では省略）
15
性能評価：アルゴリズムの効果
5000
実行時間（秒）
4000
3000
従来法，提案法ともに CSX600 を利用
その他
小行列の固有値計算
対角ブロックの更新
非対角ブロックの更新
2000
1000
0
従来法
提案法
従来法
提案法
（n = 3000, m = 160, q = 4） （n = 6000, m = 200, q = 5）
• 従来法に比べ，提案法は 約 1.4 倍高速
– 対角ブロックの更新：約 1.5 倍高速化（演算量増加）
– 小行列の固有値計算：10 倍以上高速化（演算量削減）
16
性能評価：CSX600 の効果
100000
実行時間（秒）
80000
60000
その他
小行列の固有値計算
対角ブロックの更新
非対角ブロックの更新
40000
n = 12000
20000
n = 6000
0
（無）
（有）
（有）
（無）
（有）
（有）
m = 100
q=5
m = 200
q=5
m = 100
q=5
m = 240
q=6
• CSX600 （+ 提案法） により
– n = 6000 のとき，約 3.5 倍高速化
– n = 12000 のとき，約 3.8 倍高速化
17
まとめ
• 本研究では コプロセッサ向けに QR 法の高速化を行った．
• 従来法に対して
– 再帰型アルゴリズム + ブロック化を導入
– 行列乗算以外の演算部分を高速化
– 提案法は 約 1.4 倍高速
提案法は，コプロセッサ向けの
有力な固有値解法となる可能性がある．
18
今後の課題
• パラメータ最適化
–
–
–
–
シフト数 m
chasing 行数 k
ブロック化数 q
再帰段数 d
（最適に近い値）
m = 200
d=1
5000
4000
– 演算量最小化 ≠
実行時間最小化
3000
2000
– 問題規模・計算機環境依存
1000
0
性能予測モデルによる最適化
q
1
5
10
20
25
500
400
300
200
100
k
19
20
CSX600 の性能（行列乗算）
30
K
N
•
GFLOPS
×
M
25
10
0
1000 2000 3000 4000 5000 6000
N
30
N
×
GFLOPS
25
K
M
15
5
長辺 ＞ 5000
– 性能飽和
短辺 ＞ 1000
– 高速化の効果大
•
M = K = 500
M = K = 1000
M = K = 1500
M = K = 2000
20
N = K = 500
N = K = 1000
N = K = 1500
N = K = 2000
20
15
10
5
0
1000 2000 3000 4000 5000 6000
M
21
数値実験
16000
n = 6000
実行時間（秒）
14000
12000
その他
10000
小行列の固有値計算
対角ブロックの更新
8000
非対角ブロックの更新
固有ベクトル計算
6000
ヘッセンベルグ化
4000
2000
従来法
0
m = 100 （無）
m = 120 （有）
従来法 提案法 従来法 提案法
+
CSX600
提案法
m = 100, q = 5 （無）
m = 200, q = 5 （有）
22