Transcript 近況報告 耐遅延ガウス消去

大域分散環境に適した
連立一次方程式の解法
東京大学 情報理工学系研究科
戦略ソフトウェア・電子情報 田浦研
遠藤敏夫
COEワークショップ
大域向け連立一次・遠藤敏夫
1
大規模連立一次方程式の重要性

流体解析，構造解析・・・
 a 11

 a 21
 

a
 n1


a 12
a 22

a n1

A
a 1 n   x1   b1 
   
a 2 n  x 2   b2 
 






   
a nn   x n   b n 
x= b
シミュレーションの詳細化
⇒ 方程式の大規模化
COEワークショップ
大域向け連立一次・遠藤敏夫
2
大規模計算へのアプローチ(1)

スーパーコンピュータ




地球シミュレータ，Blue Gene・・・
京速計算機(2010年)
高価，大物は世界に数台
クラスタ



比較的安価
予算・熱・電力の問題により台数限界 情報理工COE ISTBSクラスタ
低電力クラスタの動きも
高熱の問題
(RAID観測カメラより)
COEワークショップ
大域向け連立一次・遠藤敏夫
3
大規模計算へのアプローチ(2)
グリッドコンピューティング
大域分散した多数の計算機を活用
 複数組織のクラスタ・スパコンを協調
利用



Internet
Globus, Condor-G, Ninf-G, Grid MPI…
TeraGrid, NAREGI, Grid5000,
PlanetLab…
家庭等のPCを協調利用

SETI@home, Folding@home, BOINC…
Internet
COEワークショップ
大域向け連立一次・遠藤敏夫
4
本研究の目標
Ax=b
Ax=b
Interne
t
Gw=h
Cy=d
Interne
t
Ez=f
Ax=b
Interne
t
Interne
t

主流: 各サイトで別の
計算
COEワークショップ

単一の大規模計算を
大域分散で行いたい
大域向け連立一次・遠藤敏夫
5
連立一次を大域分散環境で
解くのはなぜ困難?

「大量に」「頻繁に」通信が必要な並列アルゴリズムだから
⇒ ネットワーク性能の低い大域分散環境では困難

「大量に」は，バンド幅向上により解決可能



SuperSINETの10Gbps化，一般家庭への光ファイバー普及…
「long fat pipe」のための通信技術
「頻繁に」は，待っていても解決されない

北海道-沖縄間の通信遅延は10ms以下にはならない(光速限界)
⇒ 大遅延に耐えられる並列アルゴリズムが重要に
COEワークショップ
大域向け連立一次・遠藤敏夫
6
連立一次の様々な解法
ガウス消去法
スカイライン・・・
直接法
ヤコビ，ガウスザイデル
SOR・・・
定常
反復法
非定常

今回の発表
共役勾配法(CG)
BiCG, CGS, BiCGStab・・・
GMRES, ORTHOMIN・・・
各手法につき，効率化・並列化の研究が超多数あり
COEワークショップ
大域向け連立一次・遠藤敏夫
7
発表のアウトライン

ガウス消去法
(SWoPP05, Grid05で発表)



Batched pivotingの提案
大遅延環境での実験
CG法

Block CG法を用いた実験
COEワークショップ
大域向け連立一次・遠藤敏夫
8
ガウス消去法(密行列用)とは


基本は中学で習う方法
応用分野



流体解析，構造解析 ・・・？
Top500スーパーコンピュータランキング(Linpack)
ピボット選択が計算精度の肝
COEワークショップ
大域向け連立一次・遠藤敏夫
9
主流: ガウス消去+Partial pivoting
(GE+PP)
for k = 1 to n
ピボット選択
第k列中で絶対値最大の値(ピボッ
ト)を選ぶ
行交換
n
更新
a ij  a ij  a ik a kj / a kk
ピボットは分母として使われるため，
絶対値の大きいものが良い
COEワークショップ
大域向け連立一次・遠藤敏夫
10
並列GE+PPの問題

ノード数p=6 (=2x3)
二次元ブロックサイクリック分割が
主流
利点：通信量が少
O(n 2 p )
毎回のピボット選択の度に
同期的通信
n
大遅延環境ではボトルネックに！
sb
COEワークショップ
大域向け連立一次・遠藤敏夫
11
大遅延があるときの
GE+PPの性能
Linuxクラスタ上で，遅延を
エミュレートして実験



全ノード間に同一遅延
+0ms, +2ms, +5ms,
+10ms
High performance Linpack
(HPL) 利用


行列サイズn=32768
ノード数64 (=8x8)
10msの遅延では
はるかに遅い!
1400
Execution time (sec)

1200
1000
800
600
400
200
やはりGE+PPは大遅延に弱い
COEワークショップ
大域向け連立一次・遠藤敏夫
0
+0
+2
+5
+10
Added latency (msec)
Partial
12
他のpivoting手法はどうか?
まじめ
Complete
Rook
[Neal92]
耐遅延でない
Partial
Threshold
[Malard91] etc.
Batched (提案手法)
Pairwise
さぼる
COEワークショップ
[Sorensen85] etc.
No pivoting
大域向け連立一次・遠藤敏夫
誤差の影響激しい
13
Batched Pivoting (BP)のねらい

複数回(d回)分のピボット選択を「まとめる」

同期通信回数を1/d に削減
耐遅延性の向上！
COEワークショップ
大域向け連立一次・遠藤敏夫
14
Batched Pivoting アルゴリズム (1)
d 回分のピボット選択アルゴリズム



右図ではd 行はP1とP2で分割
各ノードは自担当分を複製
複製に対して局所的，投機的に
GE+PP

d
d
その仮定でd個のピボット候補がみ
つかる
sb
P1
P1
P3
P2
P4
P1
P3
P2
P4
P2
sb
COEワークショップ
大域向け連立一次・遠藤敏夫
15
Batched Pivoting アルゴリズム (2)


ピボット候補の集合を通信で集める
「最良」候補集合を選ぶ


悪いピボットを避けたい・・・
図：d=4の例
P1
I recommend
4.8 on 50th row,
-2.5 on 241th row,
4.3 on 285th row,
-3.6 on 36th row
P2
比較
I recommend
-9.2 on 310th row,
6.8 on 121th row,
0.8 on 170th row,
-5.9 on 146th row
採用!

最良候補集合を用いて計算を続ける
COEワークショップ
大域向け連立一次・遠藤敏夫
16
Partial pivotingとの比較

選ばれるピボット



PPは各ステップで独立にピボットを選ぶ
BPでは，連続するdステップのピボットは単一ノード出身
選択の幅が狭まるため，PPより悪い可能性
計算量
PP: ( 2 / 3 ) n  O ( n )
3
2
2
(
2
/
3
)
n

O
(
n
)

O
(
dn
)
BP:
3


2
局所GEのため
d<<n なら，差は小さい
COEワークショップ
大域向け連立一次・遠藤敏夫
17
他の手法との比較

Threshold pivoting [Malard91] etc.




絶対値最大でなくても， a pk    max a ik
( 0    1)
Good: 行交換の通信量を削減
Bad: 耐遅延ではない
であればピボットになりうる
Pairwise pivoting [Sorensen85] etc.



次々に隣接した2行を取り出し，うち1行を消去 (cf. バブルソート)
Good: ピボット選択のパイプライン化可能⇒耐遅延
Bad: 計算精度悪い
COEワークショップ
大域向け連立一次・遠藤敏夫
18
並列実験環境

192ノードLinuxクラスタ

Dual Xeon 2.4/2.8GHz




ノードあたり1 CPU使用
Gigabit ethernet
クラスタ内片道遅延: 55—75 us
HPLの改造によりBPを実装


mpich 1.2.6
BLAS library by Kazushige Goto
情報理工COE
ISTBS クラスタ
COEワークショップ
大域向け連立一次・遠藤敏夫
19
Execution time (sec)
基本的な並列性能

400
350
300
250
200
150
100
50
0
クラスタでの実験(遅延は普通)




32
96 128 160
64
The number of nodes

同等のスケーラビリティ
BPはdの大きさに伴いオーバ
ヘッドあり

Partial
Batched(16)
COEワークショップ
Batched(4)
Batched(64)
PPとBP (d=4, 16, 64)比較
32 から 160 ノード
n=32768, sb=256
d=64 で7から15%
大域向け連立一次・遠藤敏夫
20
遅延が大きいときの性能
遅延が大きいとき，
BPがはるかに有利

Execution time (sec)
1400
遅延をエミュレート

1200
1000

800

600
400
200

0
+0
+2
+5
+10
Added latency (msec)
Partial
Batched(16)
COEワークショップ
全ノード間に+2ms, +5ms,
+10ms
64(=8x8) ノード
n=32768, sb=256
BPは遅延に耐久可能!

d が大きいほど耐久
Batched(4)
Batched(64)
大域向け連立一次・遠藤敏夫
21
計算精度の評価方法

Partial, Batched, Threshold, Pairwise, No pivoting を
比較



1ノードで実験
BPでは，サイズ64のブロックをノードと見なす
各条件につき，100個の乱数行列
行列サイズ 128 ～ 2048
 正規化残差 || A ~
x  b ||  / || A ||  || ~
x ||  n   を評価
~
 53
 x : 計算結果, ε: マシンイプシロン(= 2 )
 100回の実験の平均

COEワークショップ
大域向け連立一次・遠藤敏夫
22
正規化残差
計算精度の評価
0.1
10

1

0.01
0.1

0.01
0.001
0.001
100
100
PP は精度良好
No pivoting, Pairwise は
精度悪
BPはPPに匹敵する精度！

1000
1000
行列サイズ n
Partial
Partial
Batched(16)
Batched(16)
Threshold(0.5)
Threshold(0.5)
Pairwise
Pairwise
COEワークショップ
d の値により，精度と耐遅延
性のトレードオフ
10000
Batched(4)
Batched(4)
Batched(64)
Batched(64)
Threshold(0.25)
Threshold(0.25)
No
No
大域向け連立一次・遠藤敏夫
23
ガウス消去法についてまとめ

ガウス消去法を遅延に強くする，
batched pivotingを提案


耐遅延性と計算精度を両立
最速アルゴリズム，最適パラメータは状況に依存

遅延，ノード数，問題サイズ・・・
COEワークショップ
大域向け連立一次・遠藤敏夫
24
発表のアウトライン

ガウス消去法



Batched pivotingの提案
大遅延環境での実験
CG法

Block CG法を用いた実験
COEワークショップ
大域向け連立一次・遠藤敏夫
25
CG法とは

反復法


n次元ベクトルx が，次々に更新されて真の解へ近づく
CG法は，正定値対称行列向けの反復法のメジャー
な手法 [Hestenes&Stiefel 52]

非常に非常に多くの効率化手法が提案されている
x(2)
x(1)
x(3)
(4)
x
x(5)
COEワークショップ
大域向け連立一次・遠藤敏夫
26
CG法のアルゴリズム
x
(0)
 0; r
(0)
 b; p
r
(0)
) /( p
(k )
(0)
for k =0, 1, 2…

(k )
,r
(k )
 x
(k )

(k )
p
( k 1)
r
(k )

(k )
Ap

p
r
(k )
( k 1)
x
r
 (r
(k )
 (r
( k 1)
( k 1)
( k 1)
r
,r
( k 1)
0
( k 1)

, Ap
)
(k )
(k )
) /( r
(k )
(k )
（素人から見た）特徴：
 「あと探索すべき空間」が，
1ステップで1次元減っていく
p
(k )
(k )
ならば終了
,r
(k )
)

理論上はn 回の反復で真
の解へ到達


COEワークショップ
たいていはもっと早く収束
丸め誤差のため遅いことも
大域向け連立一次・遠藤敏夫
27
CG法と大遅延環境の関係
x
(0)
 0; r
(0)
 b; p
r
(0)
) /( p
(k )
(0)
ベクトル内積
行列・ベクトル積
for k =0, 1, 2…

(k )
,r
(k )
 x
(k )

(k )
p
( k 1)
r
(k )

(k )
Ap

p
r
(k )
( k 1)
x
r
 (r
(k )
 (r
( k 1)
( k 1)
( k 1)
r
,r
( k 1)
0
COEワークショップ
( k 1)

, Ap
)

(k )
(k )
) /( r
(k )
(k )
p
(k )
(k )
ならば終了
,r
(k )
行列・ベクトル積は耐遅
延化可能[Allen 01]

)
領域オーバラップ
内積の結果共有のため
に毎ステップ同期が必要
⇒遅延に弱い！

大域向け連立一次・遠藤敏夫
28
Block CG法の採用

Block CG法
[O’Leary 80]など

s 個の方程式を同時に解くことができる

あえて1個の方程式に利用

B の残り(s-1) 列には
Block CG法
CG法
A
COEワークショップ
ダミー値
×x = b
大域向け連立一次・遠藤敏夫
A
s
×X = B
29
Block CG法のアルゴリズム
X
(0)
 0; R
 B; P
(0)
(0)
R
(0)

for k =0, 1, 2…

X
R

P
R
(k )
 (P
( k 1)
( k 1)
(k )
 X
 R
 (R
( k 1)
( k 1)
( k )T
(k )
COEワークショップ
P
R
( k 1)
0
(k )
(k )
1
) (R
(k )
 AP
(k )
( k )T
 R
AP

(k )
) (R
P
(k )

R
(k )
)
(k )

1
( k )T

(k )
( k 1) T
R
( k 1)
)
要するに，n次元ベクトル
x, b, r, p が，n×s 行列
へ置き換わっている
依然，毎ステップ同期が
必要
(k )
ならば終了
大域向け連立一次・遠藤敏夫
30
Block CG法の特徴


1ステップの計算量はs 倍強に
理論上の反復回数がn/s


「あと探索すべき空間」が， 1ステップでs 次元減
CG法の1/s 倍
⇒ 合計計算量はやや増
3
3
2
O
(
n

sn
)
O
(
n
)
CG:
Block CG:
同期通信を減らせるので
大遅延環境で得ができそう
COEワークショップ
大域向け連立一次・遠藤敏夫
31
実験

CG, Block CGの実行時間比較



行列はUF sparse matrix collectionより



C++で実装，最適化の余地多
クラスタの16ノード利用，遅延エミュレート
msc10848 (n=10848)
cfd2 (n=123440)
実行条件




RCM(reverse Cuthill-Mckee)オーダリング
行を均等分散
スケーリング
8
残差が r 2 / r0 2  10
となるまで反復
COEワークショップ
大域向け連立一次・遠藤敏夫
32
並列実行時間の比較
実行時間(s)
250
200
大遅延の
とき有利
150
100
50
0
600
400
200
+5
+10
追加遅延(ms)
CG

800
0
+0

もともと
遅い
cfd2
実行時間(s)
msc10848
BCG(4)
BCG(16)
+0
CG
+10
+5
追加遅延(ms)
BCG(4)
BCG(16)
Block CGは遅延の増大に強い
もともとが遅すぎる場合も(行列に依存)
COEワークショップ
大域向け連立一次・遠藤敏夫
33
実際の反復回数は1/sになるか？
行列
前処理
CG
BCG(4)
BCG(16)
scaling
5541
1731
602
(1/3.2)
(1/9.2)
307
94
(1/4.7)
(1/15.3)
5400
2816
(1/1.6)
(1/3.0)
1867
1027
(1/1.6)
(1/2.8)
31560
8304
(1/4.0)
(1/15.0)
13080
2900
(1/3.7)
(1/16.8)
msc10848
(n=10848)
IC
scaling
1437
8411
cfd2
(n=123440) 加速IC(2)
scaling
2910
124827
ct20stif
(n=52329)
加速IC(2)
COEワークショップ
48855
大域向け連立一次・遠藤敏夫
まぁまぁ
不満
IC: 不完全コレ
スキー分解
34
残差の推移の比較
msc10848
1.E+02
1.E+00
相対残差
相対残差
1.E+00
1.E-02
1.E-04
1.E-02
1.E-04
1.E-06
1.E-06
1.E-08
1.E-08
1
1001 2001 3001 4001 5001
反復回数
CG

cfd2
1.E+02
BCG(4)
BCG(16)
1
2001
CG
4001
6001
反復回数
BCG(4)
8001
BCG(16)
CGのグラフと分岐してからの落ち方が違う？
COEワークショップ
大域向け連立一次・遠藤敏夫
35
CG法についてまとめ

CG法もそのままでは遅延に弱いため，Block CG法
[O’Leary 80]を用いて実験

同期通信回数を大幅に減らし，耐遅延性が高い
大域環境で構造解析，流体解析ができるかも

使いどころ，最適パラメータ選択は難しい


反復回数の減り方が激しく行列依存
より優れたオーダリング，前処理との関係は？
COEワークショップ
大域向け連立一次・遠藤敏夫
36
おまけ：本郷・柏間での実験


情報理工クラスタ(本郷)と近山研クラスタ(柏)
1クラスタ
2クラスタ
16ノード
8+8ノード
CG
27.5(s)
510(s)
BCG(4)
13.6
321
BCG(16)
20.4
198
• 片道遅延は2—3ms
• 行列はmsc10848
• MPI/GXP2[斎藤06]利用
エミュレーション通りに行かない・・・



物理的バンド幅の差？⇒1Gbpsのはずなので悪くない
クラスタ間TCPが1本あたり11Mbpsしか出ていない
分割や通信順序のチューニング不足，etc, etc.
COEワークショップ
大域向け連立一次・遠藤敏夫
37
おわりに

大域分散環境で大規模計算を実現するために
は，インフラ・ハード・ミドルウェア・アプリ全ての
層の研究が必要
Ax=b
Internet
COEワークショップ
Ax=b
Internet
大域向け連立一次・遠藤敏夫
38