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

キャッシュマシン向け三重対角化
アルゴリズムの性能予測方式
名古屋大学 計算理工学専攻
山本有作
目次
1. はじめに
2. 三重対角化アルゴリズムとキャッシュ向け最適化
3. 本研究の目的
4. 性能予測手法
5. 実験結果
6. 関連研究
7. まとめと今後の課題
1. はじめに
• 本研究で対象とする問題
– 標準固有値問題 Ax = λx
• A : 実対称またはエルミートの密行列
• 応用分野
– 分子計算,統計計算など
• 分子計算では,10万元以上の超大規模行列の固有値計算も
需要あり
固有値計算の流れ
実対称行列 A
計算内容
計算手法
三重対角化
Q*AQ = T
(Q: 直交行列)
ハウスホルダー法
三重対角行列 T
三重対角行列の
固有値・固有ベクトル計算
Tyi =λi yi
Tの固有値 {λi },
固有ベクトル {yi }
逆変換
Aの固有ベクトル {xi }
xi = Qyi
QR法
分割統治法
二分法・逆反復法
MR3アルゴリズム
逆変換
固有値計算の処理時間の内訳 (P4クラスタ,16ノー
ド)
Paolo Bientinesi et al.: A Parallel
Eigensolver for Dense Symmetric
Matrices using Multiple Relatively
Robust Representations (2004) より
三重対角行列の固有値・固有ベクトル計算に最新のアルゴリズム
(MR3)を使った場合,処理時間の大部分を三重対角化が占める。
ハウスホルダー法の高速化が必要
2. 三重対角化アルゴリズムとキャッシュ向け最適化
• ハウスホルダー法による三重対角化
ベクトル
– 鏡像変換 H = I – a uut による列の消去
• Hは対称な直交行列
• 与えられたベクトルの第1成分以外をゼロにする。
左からH
を乗算
• 第 k ステップでの処理
0
0
0
0
0
左からH
を乗算
0
右からH
を乗算
k
非ゼロ要素
ゼロにしたい部分
影響を受ける部分
ハウスホルダー法のアルゴリズム
[Step 1] k = 1からN –2まで以下の[Step 2] ~ [Step 8]を繰り返す。
[Step 2] 内積 σ(k) = sqrt(d (k)t d (k)) の計算
[Step 3] 鏡像変換ベクトル作成:
u(k) = (d (k)1 – sgn(d (k)1)σ(k), d (k)2, … , d (k)N – k )
[Step 4] 規格化定数の計算: α(k) = 2 /∥u(k)∥2
[Step 5] 行列ベクトル積: p(k) =α(k)C (k)u(k)
[Step 6] ベクトルの内積: β(k) =α(k) u(k)tp(k) /2
[Step 7] ベクトルの加算: q(k) = p(k) –β(k)u(k)
[Step 8] 行列のrank-2更新: C (k) = C (k) – u(k)q(k)t – q(k)u(k)t
0
0
d (k )1
C (k )
d (k )
k
ピボット列
ピボット行
ハウスホルダー法の特徴と問題点
• 特徴
– 全演算量のほとんどが,level-2 BLAS である行列ベクトル積と
行列のrank-2更新で占められる。
• 全演算量:
約 (4/3)N3
• 行列ベクトル積の演算量:
約 (2/3)N3
• rank-2更新の演算量:
約 (2/3)N3
• 問題点
– level-2 BLASのため,行列データの再利用性なし
• キャッシュミス,SMPでのメモリ競合の影響大
– 最近のマイクロプロセッサ上での高性能が期待できない。
キャッシュマシン向けの三重対角化アルゴリズム
• 原理(Bishof et al., 1993, 1994)
– 密行列 A をまず半帯幅Lの帯行列 B に変換し,次にこの帯行列
を三重対角行列 T に変換する。
次数 N
半帯幅 L
帯行列化
約 (4/3)N3
A
0
村田法
O(N2L)
0
B
0
0
T
• 三重対角化を2段階で行うことの利点
– 帯行列への変換は, level-3 BLAS のみを使って実行可能。
– 帯行列から三重対角行列への変換の演算量はO(N2L)であり,
前半部に比べてずっと小さい。
帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – UαUt
• Hは直交行列
• 与えられたブロックベクトルの
第1ブロックを上三角行列にし,
それ以外をゼロにする。
左からH
を乗算
第kステップでの処理
0
0
左からH
を乗算
非ゼロ要素
0
0
0
ゼロにしたい部分
右からH
を乗算
0
影響を受ける部分
帯行列化のアルゴリズム(続き)
[Step 1] K = 1からN /L–1まで以下の[Step 2] ~ [Step 6]を繰り返す。
[Step 2] D(K) を第1ブロックが上三角行列で第2ブロック以下が
ゼロ行列であるブロックベクトルに変換するブロック鏡像変換
I – U(K) (K)U(K)tを求める。
[Step 3] 行列・ブロックベクトル積: P(K) = C (K)U(K)α(K)
[Step 4] ブロックベクトルの内積: β(K) = α(K)tU(K)tP(K) / 2
[Step 5] ブロックベクトルの加算: Q(K) = P(K) –U(K)β(K)
[Step 6] 行列のrank-2L更新: C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t
すべて level-3 BLAS(行列乗算)
• 本アルゴリズムの特徴
– 演算が level-3 BLAS 中心のため,キャッシュの有効利用が可能
– SMPにおけるメモリ競合の影響を低減可能
多段化による性能向上
• 原理(Wu et al., 1996)
– 帯行列化のアルゴリズムにおいて,複数のピボット列・行を溜め
ておき,後でまとめて更新
– M 本溜めることで,rank-2L 更新の部分を rank-2LM 更新にでき,
データ再利用性を更に向上可能
Opteron (1.6GHz) 上での三重対角化の性能
1800
Performance (GFLOPS)
1600
1400
L=24, M=4
Dongarra
本アルゴリズム
本アルゴリズム+多段化
L=48
1200
1000
800
M=32
600
400
200
0
480
960
1920
3840
Matrix size
本アルゴリズムは他の方法に比べて非常に高性能
N = 3840 のとき,ピークの50%以上の性能を達成
本アルゴリズムの課題
• 性能を大きく左右するパラメータが2つ存在
– 帯行列の半帯幅 L
– 多段化の段数 M
• 各パラメータに関するトレードオフ
– L が大 → BLAS性能は上がるが,村田法の演算量が増大
– M が大 → BLAS性能は上がるが,帯行列化の演算量が増大
対象とするプロセッサおよび行列サイズ N に応じて,L と M を
最適化することが必要
3. 本研究の目的
• キャッシュマシン向けの三重対角化アルゴリズムに対し,高精度な
性能予測モデルを構築する
最適な性能パラメータ L,M の値を実行前に推定する
自動チューニング型ライブラリのための基礎技術
• 性能予測モデルの他の応用
– グリッド上での最適マシン選択
– グリッド上でのスケジューリング
4. 性能予測手法
• 階層的な性能モデリング(Cuenca et al., 2004)
– 帯行列化のアルゴリズムは,多種類のlevel-3 BLASへのコール
から成る
– level-3 BLASの各ルーチンの性能を精度良くモデル化できれば,
その積み上げによりアルゴリズム全体の性能も精度良くモデル
化できるはず
帯行列化で用いるlevel-3 BLASルーチン
• アルゴリズム
[Step 1] K = 1からN /L–1まで以下の[Step 2] ~ [Step 6]を繰り返す。
[Step 2] ブロック鏡像変換 I – U(K) (K)U(K)t の計算
BlockHouse
[Step 3] 行列・ブロックベクトル積: P(K) = C (K)U(K)α(K)
DSYMM
[Step 4] ブロックベクトルの内積: β(K) = α(K)tU(K)tP(K) / 2
DGEMM
[Step 5] ブロックベクトルの加算: Q(K) = P(K) –U(K)β(K)
DGEMM
[Step 6] 行列のrank-2L更新: C (K) = C (K) – U(K)Q(K)t – Q(K)U(K)t
DSYR2K
• 必要なBLASルーチン
–
–
–
–
DGEMM 3種 (‘N’,’N’, ‘N’,’T’, ‘T’,’N’)
DSYMM
DSYR2K
BlockHouse (BLASではないがブロック鏡像変換を求めるルーチン)
BLAS性能のモデリング
• 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 を1つの多項式で近似すると,誤差が大きい
– 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 に対する値は,一次補間により計算
BLAS性能のモデリング(続き)
• DGEMM(’N’,’T’) および DGEMM(’T’,’N’)
– ‘N’,’N’の場合と同様にして 実行時間の予測関数を構成
• DSYMM
– サイズを決めるパラメータは n, m の2つのみ
– m の代表的な値に対し,実行時間を n の一次式として近似
• DSYR2K
– DSYMM の場合と同様
• BlockHouse
– サイズを決めるパラメータは,ベクトル長 n とブロック幅 L の2つ
– level-3 BLAS の場合と異なり,一次式/双一次式近似では実行時間が
うまく再現できない
– (n, L)平面上の格子点で実行時間を測定し,補間により実行時間を予測
帯行列化アルゴリズムの性能モデリング
• 基本的な考え方
– BLASの各ルーチンと同じ引数を持ち,演算を行う代わりに実行時間の
予測値のみを計算して返すルーチンを作成
– 帯行列化プログラム中のBLASをこれらのルーチンで置き換えることに
より,帯行列化の実行時間を予測するプログラムを作成
計算プログラム
予測プログラム
DO K=1, N/L
CALL BlockHouse(m,L,...)
CALL DSYMM(m,n,...)
CALL DGEMM(m,n,k,...)
CALL DGEMM(m,n,k,...)
CALL DSYR2K(m, k,...)
END DO
DO K=1, N/L
T1=T1+BlockHouse_TIME(m,L,...)
T1=T1+DSYMM_TIME(m,n,...)
T1=T1+DGEMM_TIME(m,n,k,...)
T1=T1+DGEMM_TIME(m,n,k,...)
T1=T1+DSYR2K_TIME(m, k,...)
END DO
• 本方式のメリット
– 複雑な解析的モデルの構築が不要
– 予測に必要な時間は O(N/L)
5. 数値実験
• 評価項目
– BLAS性能の予測
– 帯行列化アルゴリズムの性能の予測
– パラメータ L,M の最適化
• 評価環境
– Opteron(1.6GHz),g77,GOTO BLAS
– Alpha 21264A(750MHz),g77,GOTO BLAS
– Power PC G5(2.0GHz),g77,ATLAS 3.7.0
• 評価例題
– N=480 ~ 9600 の対称密行列の三重対角化
BLAS性能の予測
• 前述の方式で各 BLAS について実行時間の予測関数を構成し,各
プロセッサ上で様々な行列サイズに対して実測データと比較
• すべてのサイズパラメータ(m, n, k など)が6を超える領域では,どの
プロセッサ/BLASでも予測誤差は5~10%程度と高精度
予測時間/実測時間
予測時間/実測時間
1.20E+00
1.20E+00
1.00E+00
3
6
12
24
48
96
8.00E-01
6.00E-01
4.00E-01
1.00E+00
3
6
12
24
48
96
8.00E-01
6.00E-01
4.00E-01
2.00E-01
2.00E-01
0.00E+00
0.00E+00
3
6
3
12
k
100
24
200
48
96
400
800
m
DGEMM(‘N’,’N’), Opteron, n=12
6
12
k
100
24
200
48
96
400
800
m
DGEMM(‘N’,’N’), G5, n=12
帯行列化アルゴリズムの性能の予測
• L,M を変えたときの予測時間と実測時間(Opteron, N=1920)
時間(s)
時間(s)
18
18
16
L=3
L=6
L=12
L=24
L=48
L=96
14
12
10
8
6
16
L=3
L=6
L=12
L=24
L=48
L=96
14
12
10
8
6
4
4
1
2
4
0
L=3
L=6
16
L=12
L
L=24
L=48
64
L=96
予測時間
M
1
2
4
0
L=3
L=6
16
L=12
L
L=24
L=48
64
L=96
実測時間
• 予測誤差は10%以内
• モデルは L,M を変えたときの実行時間の変化をよく再現
M
帯行列化アルゴリズムの性能の予測(続き)
• L,M を変えたときの予測時間と実測時間(Alpha 21264A, N=3840)
時間(s)
時間(s)
400
400
350
L=3
L=6
L=12
L=24
L=48
L=96
300
250
200
150
100
350
L=3
L=6
L=12
L=24
L=48
L=96
300
250
200
150
100
1
50
4
0
L=3
L=6
16
L=12
L
L=24
L=48
64
L=96
予測時間
M
1
50
4
0
L=3
L=6
16
L=12
L
L=24
L=48
64
L=96
実測時間
• 予測誤差は5%以内
• モデルは L,M を変えたときの実行時間の変化をよく再現
M
帯行列化アルゴリズムの性能の予測(続き)
• L,M を変えたときの予測時間と実測時間(PowerPC G5, N=1920)
時間(s)
時間(s)
25
25
L=3
L=6
L=12
L=24
L=48
L=96
20
15
10
5
1
4
0
L=3 L=6
16
L=12
L
L=24
L=48
64
L=96
予測時間
M
L=3
L=6
L=12
L=24
L=48
L=96
20
15
10
5
1
4
0
L=3 L=6
16
L=12
L
L=24
L=48
64
L=96
実測時間
• 予測誤差は10%以内
• モデルは L,M を変えたときの実行時間の変化をよく再現
M
パラメータ L,M の最適化
• N が与えられたとき,(L,M)の各組について,三重対角化の実行時
間を次の2つの時間の和として予測
– 帯行列化の実行時間: モデルにより予測
– 村田法の実行時間:
実測値
この予測値を最小とする(L,M)を決定
パラメータ L,M の最適化(続き)
•
予測した(L,M)による実行結果(G5,N=960~9600)
3500
性能(MFLOPS)
3000
2500
2000
1500
(12,8)
(24,4)
(48,4)
auto-tuned
optimal
1000
500
0
960
•
•
1920 2880 3840 4800 5760 6720 7680 8640 9600
行列サイズ N
予測した(L,M)での性能は,最適な(L,M)での性能にほぼ一致
単一の(L,M)では,全行列サイズに対して最高性能を発揮することは不可能
本モデルを用いたパラメータの事前最適化は有効
6. 関連研究
• 三重対角化プログラムの自動チューニング(Katagiri et al., 2000)
– 分散メモリ向けの三重対角化プログラムにおいて,ループ展開の段数,
通信関数の種類などのパラメータを自動的に最適化する方式
– パラメータを変化させてプログラム全体を何回も実行することにより最適
値を求めるため,最適化に時間がかかるという問題点がある
• 階層的な性能モデリング(Cuenca et al., 2004,Dackland et al., 1996)
– 線形計算プログラムの自然な階層構造を利用して,下位のルーチン
(BLASなど)の性能モデルをまず構築し,それを用いて順次上位のルー
チンの性能モデルを構築していく方式
– 本研究のモデルもこの考え方に基づく
– ただし,従来の適用例は,LU分解やQR分解などの基本的分解,および
ヤコビ法などの単純なアルゴリズムに限られている
7. まとめと今後の課題
• 本研究のまとめ
– キャッシュマシン向けの三重対角化アルゴリズムに対し,階層的なモ
デリング手法に基づく性能予測モデルを開発した
– Opteron,Alpha 21264A,Power PC G5上で評価した結果,予測誤差
は5~10%であり,L,M を変えたときの実行時間の変化も良く再現で
きた
– 本モデルを用いて L,M の最適値を推定し,三重対角化を実行したと
ころ,真の最適値での性能とほぼ等しい性能が得られた
• 今後の課題
–
–
–
–
より多くの種類のプロセッサ上での検証
共有メモリ/分散メモリ型並列計算機向けの拡張
性能安定化手法の適用
自動チューニング型ライブラリへの適用
謝辞
• 日頃から有益な議論をして頂いている自動チューニング
研究会のメンバーに感謝いたします