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

キャッシュマシン向け三重対角化
アルゴリズムの性能予測方式
名古屋大学計算理工学専攻
山本有作
目次
１. はじめに
２. 三重対角化アルゴリズムとキャッシュ向け最適化
３. 本研究の目的
４. 性能予測手法
５. 実験結果
６. 関連研究
７. まとめと今後の課題
１. はじめに
• 本研究で対象とする問題
– 標準固有値問題 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）を使った場合，処理時間の大部分を三重対角化が占める。
ハウスホルダー法の高速化が必要
２. 三重対角化アルゴリズムとキャッシュ向け最適化
• ハウスホルダー法による三重対角化
ベクトル
– 鏡像変換 H = I – a uut による列の消去
• Hは対称な直交行列
• 与えられたベクトルの第１成分以外をゼロにする。
左から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
• 三重対角化を２段階で行うことの利点
– 帯行列への変換は， level-3 BLAS のみを使って実行可能。
– 帯行列から三重対角行列への変換の演算量はO(N2L)であり，
前半部に比べてずっと小さい。
帯行列化のアルゴリズム
ブロックベクトル
ブロック鏡像変換によるブロック列の消去
– ブロック鏡像変換 H = I – UαUt
• Hは直交行列
• 与えられたブロックベクトルの
第１ブロックを上三角行列にし，
それ以外をゼロにする。
左から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) を第１ブロックが上三角行列で第２ブロック以下が
ゼロ行列であるブロックベクトルに変換するブロック鏡像変換
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%以上の性能を達成
本アルゴリズムの課題
• 性能を大きく左右するパラメータが２つ存在
– 帯行列の半帯幅 L
– 多段化の段数 M
• 各パラメータに関するトレードオフ
– L が大 → BLAS性能は上がるが，村田法の演算量が増大
– M が大 → BLAS性能は上がるが，帯行列化の演算量が増大
対象とするプロセッサおよび行列サイズ N に応じて，L と M を
最適化することが必要
３. 本研究の目的
• キャッシュマシン向けの三重対角化アルゴリズムに対し，高精度な
性能予測モデルを構築する
最適な性能パラメータ L，M の値を実行前に推定する
自動チューニング型ライブラリのための基礎技術
• 性能予測モデルの他の応用
– グリッド上での最適マシン選択
– グリッド上でのスケジューリング
４. 性能予測手法
• 階層的な性能モデリング（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 を１つの多項式で近似すると，誤差が大きい
– 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)
５. 数値実験
• 評価項目
– 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）では，全行列サイズに対して最高性能を発揮することは不可能
本モデルを用いたパラメータの事前最適化は有効
６. 関連研究
• 三重対角化プログラムの自動チューニング（Katagiri et al., 2000）
– 分散メモリ向けの三重対角化プログラムにおいて，ループ展開の段数，
通信関数の種類などのパラメータを自動的に最適化する方式
– パラメータを変化させてプログラム全体を何回も実行することにより最適
値を求めるため，最適化に時間がかかるという問題点がある
• 階層的な性能モデリング（Cuenca et al., 2004，Dackland et al., 1996）
– 線形計算プログラムの自然な階層構造を利用して，下位のルーチン
（BLASなど）の性能モデルをまず構築し，それを用いて順次上位のルー
チンの性能モデルを構築していく方式
– 本研究のモデルもこの考え方に基づく
– ただし，従来の適用例は，LU分解やQR分解などの基本的分解，および
ヤコビ法などの単純なアルゴリズムに限られている
７. まとめと今後の課題
• 本研究のまとめ
– キャッシュマシン向けの三重対角化アルゴリズムに対し，階層的なモ
デリング手法に基づく性能予測モデルを開発した
– Opteron，Alpha 21264A，Power PC G5上で評価した結果，予測誤差
は5～10%であり，L，M を変えたときの実行時間の変化も良く再現で
きた
– 本モデルを用いて L，M の最適値を推定し，三重対角化を実行したと
ころ，真の最適値での性能とほぼ等しい性能が得られた
• 今後の課題
–
–
–
–
より多くの種類のプロセッサ上での検証
共有メモリ／分散メモリ型並列計算機向けの拡張
性能安定化手法の適用
自動チューニング型ライブラリへの適用
謝辞
• 日頃から有益な議論をして頂いている自動チューニング
研究会のメンバーに感謝いたします

Download Report