Transcript 日本語版 - 東京大学

「地球シミュレータ」における非構造格子
アプリケーション向け前処理付反復法
Preconditioned Iterative Linear Solvers for
Unstructured Grids on the Earth Simulator
中島 研吾
[email protected]
東京大学情報基盤センター スーパーコンピューティング研究部門
海洋研究開発機構 計算システム計画・運用部
地球シミュレータ 利用講習会，2008年4月22日
08-APR22
2
概 要
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における最
適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
プロセッサの動向：スカラーとベクトル

スカラープロセッサ


08-APR22
クロック数とメモリバンド幅のギャップ
低い対ピーク性能比（<10%）
3
スカラープロセッサ
CPU-キャッシュ－メモリの階層構造
CPU
Register
FAST
Cache
小容量 (MB)
高価
大きい (1億以上のトランジスタ)
SLOW
Main Memory
08-APR22
大容量 (GB)
廉価
4
プロセッサの動向：スカラーとベクトル

スカラープロセッサ



クロック数とメモリバンド幅のギャップ
低い対ピーク性能比（<10%）
ベクトルプロセッサ

高い対ピーク性能比


そのためには・・・



08-APR22
例：地球シミュレータ，FEM型アプリケーション → >35 %
ベクトルプロセッサ用チューニング
充分長いベクトル長（問題サイズ）
比較的単純な問題に適している
5
ベクトルプロセッサ
ベクトルレジスタと高速メモリ
Vector Processor
Vector
Register
• 単純構造のDOループの並列処理
• 単純，大規模な演算に適している
do i= 1, N
A(i)= B(i) + C(i)
enddo
Very
FAST
Main Memory
08-APR22
6
プロセッサの動向：スカラーとベクトル

スカラープロセッサ



クロック数とメモリバンド幅のギャップ
低い対ピーク性能比 （<10%）
ベクトルプロセッサ

高い対ピーク性能比


そのためには・・・



08-APR22
例：地球シミュレータ，FEM型アプリケーション → >35 %
ベクトルプロセッサ用チューニング
充分長いベクトル長（問題サイズ）
比較的単純な問題に適している
7
08-APR22
8
Hitachi
Earth
SR8000/MPP
Simulator
(U.Tokyo)
IBM SP3
(LBNL)
Hitachi
SR11000/J2
(U.Tokyo)
Peak Performance
GFLOPS
8.00
1.80
1.50
9.20
Measured Memory
BW (STREAM)
(GB/sec/PE)
26.6
2.85
0.623
8.00
Estimated
Performance/Core
GFLOPS
(% of peak)
2.31-3.24
(28.8-40.5)
.291-.347
(16.1-19.3)
.072-.076
(4.8-5.1)
.880-.973
(9.6-10.6)
Measured
Performance/Core
2.93
(36.6)
.335
(18.6)
.122
(8.1)
1.34
(14.5)
BYTE/FLOP
3.325
1.583
0.413
0.870
12.3
1.60
1.00
12.0
5.6-7.7
6-20
16.3
4.7 *
Comm. BW
(GB/sec/Node)
MPI Latency
(msec)
* IBM p595
J.T.Carter
et al.
典型的な挙動
問題規模～GFLOPS：大きい程よい
3.00
1.0E+02
40 % of peak
2.00
GFLOPS
GFLOPS
2.50
1.50
1.00
1.0E+00
8 % of peak
0.50
0.00
1.0E-01
1.0E+04
1.0E+05
1.0E+06
1.0E+07
DOF: Problem Size
IBM-SP3:
問題サイズが小さい場合はキャッシュの
影響のため性能が良い
08-APR22
1.0E+01
1.0E+04
1.0E+05
1.0E+06
1.0E+07
DOF: Problem Size
Earth Simulator:
大規模な問題ほどベクトル長が長くな
り，性能が高い
9
並列計算
Ideal
PE#
08-APR22
Performance
Performance
Strong Scaling (全体問題規模固定)
Ideal
PE#
IBM-SP3:
Earth Simulator:
PE（Processing Element）数が少ない
場合はいわゆるスーパースカラー。
PE数が増加すると通信オーバーヘッド
のため性能低下。
PE数が増加すると，通信オーバーヘッ
ドに加え，PEあたりの問題規模が小さ
くなるため性能低下。
10
08-APR22
11
Improvement of Memory
Performance on
IBM SP3 ⇒ Hitachi SR11000/J2
■ Flat-MPI/DCRS
□ Hybrid/DCRS
Memory Performance
(BWDTH, latency etc.)
375 MHz
1.0 GB/sec
8M L2 cache/PE
SR11000/J2
POWER5+
3.00
15.0
2.00
10.0
GFLOPS
GFLOPS
IBM SP-3
POWER3
1.00
5.0
0.00
0.0
1.E+04
1.E+05
1.E+06
DOF
1.E+07
1.E+04
2.3 GHz
8.0 GB/sec
18M L3 cache/PE
1.E+05
1.E+06
DOF
1.E+07
08-APR22
12
中島研吾とベクトル計算機の付き合い
• Cray-1S (1985-1988)
– 三菱総合研究所
• Cray-YMP (1988-1995)
– 三菱総合研究所
– テキサス大学
• Fujitsu VP, VPP series (1985-1999)
– 日本原子力研究所（現JAEA）
– 動力炉・核燃料開発事業団（現JAEA ）
• NEC SX-4 (1997-2001)
– 日本原子力研究所 CCSE（現JAEA ）
• Hitachi SR2201 (1997-2004)
– 東京大学
– 日本原子力研究所 CCSE （現JAEA ）
• Hitachi SR8000 (2000-2007)
– 東京大学
• Earth Simulator (2002-)
08-APR22
13
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における最
適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
08-APR22
14
GeoFEM: FY.1998-2002
http://geofem.tokyo.rist.or.jp/
• 文部科学省「科学技術振興調整費総合研究」
– 「高精度の地球変動予測のための並列ソフトウェア開発に
関する研究」の一部
– リーダー：奥田洋司教授（東大・人工物）
• 固体地球シミュレーション用並列有限要素法プラッ
トフォーム
– 並列I/O，並列線形ソルバー，並列可視化をサポート
• HPCと自然科学の緊密な協力
08-APR22
15
System Configuration of
GeoFEM
Pluggable Analysis Modules
Utilities
One-domain mesh
Partitioner
Structure
Fluid
Wave
Comm.
I/F
構造計算（Static linear）
構造計算（Dynamic
構造計算（
linear）
Contact）
Solver
I/F
Vis.
I/F
Platform
Parallel Equation Visualizer
solvers
I/O
Partitioned mesh
PEs
GPPView
Visualization data
08-APR22
16
奥田，中島編「並列有限要素解析〔I〕
クラスタコンピューティング」
培風館，2004．
• 「GeoFEM」の成果のまとめ
• 「地球シミュレータ」上での最
適化，シミュレーション結果
を紹介
08-APR22
17
Results on Solid Earth Simulation
Complicated Plate Model around Japan Islands
Magnetic Field of the Earth : MHD code
Dh=5.00
Dh=1.25
Simulation of Earthquake Generation Cycle
in Southwestern Japan
T=100
T=200
T=300
T=400
T=500
Transportation by Groundwater Flow
through Heterogeneous Porous Media
TSUNAMI !!
08-APR22
18
Results by GeoFEM
08-APR22
19
FEM型アプリケーションの特徴（1/2）
• 要素単位のローカルなオペレーション
– 疎行列：差分も同じ
– 実は並列計算には向いている
• メモリーへの負担大きい
– 「不規則な」間接参照
do i= 1, N
jS= index(i-1)+1
jE= index(i)
do j= jS, jE
in = item(j)
Y(i)= Y(i) + AMAT(j)*X(in)
enddo
enddo
08-APR22
20
FEM型アプリケーションの特徴（2/2）
• 並列計算では・・・
– 基本的に通信は隣接領域とのみ（内積などは除く）
– 通信メッセージの量も比較的少ない：境界領域のみ
– MPI Latencyが結構効く
08-APR22
21
FEM型アプリケーションのベクトル化
• 非構造格子向け連立一次方程式の解法として，RCM
（Reverse Cuthil-Mckee，差分のHyperplaneに相当）が
1980年代から90年代にかけて開発。
• 浮動小数点演算部分「以外」についてはベクトル化は困
難，あるいは研究途上。
• 本講演における「ベクトル化」という用語について
08-APR22
22
本講演における「ベクトル化」という用語
• ベクトルプロセッサの特性を利用したベクトル処理，並列
処理
• このほか，「ベクトル化」による「最適化」という意味で使
用している，すなわち「最内ループ」において以下の要
件を達成していること：
– データ依存性が無い
– 分岐が無い
– ループ長がベクトルレジスタサイズ（256）より大きい
• 以下の場合でも「ベクトル化」はやってくれる（すなわち
ベクトル化率は高くなる可能性がある）が，性能は悪い：
– 分岐がある
– ループ長が短い
08-APR22
最内ループにおける分岐の有無による
性能比較：接触問題
1 SMP Node of ES, 64 GFLOPS Peak Performance
● 分岐無し, ○ 分岐有り：大きい程良い
西南日本モデル
Simple Block
25
25
分岐無し
15
分岐有り
10
15
分岐有り
10
5
5
1
23
分岐無し
20
GFLOPS
GFLOPS
20
10
100
Colors
1000
10
100
Colors
1000
08-APR22
24
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における最
適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
08-APR22
25
Earth Simulator (ES)
http://www.es.jamstec.go.jp/
• 640×8= 5,120 Vector Processors
–
–
–
–
SMP Cluster-Type Architecture
8 GFLOPS/PE
64 GFLOPS/Node
40 TFLOPS/ES
• 16 GB Memory/Node, 10 TB/ES
• 640×640 Crossbar Network
– 16 GB/sec×2
• Memory BWTH with 32 GB/sec.
• 35.6 TFLOPS for LINPACK (2002-March)
• 26 TFLOPS for AFES (Climate Simulation)
08-APR22
26
本研究の動機（motivation）
• GeoFEM Project (FY.1998-2002)
• FEM型の複雑形状に関するアプリケーション
（LINPACKや差分法ではなく）を「地球シ
ミュレータ（以下しばしばES）」上で動かす
には？
– 陰的解法（連立一次方程式）
– 並列プログラミングモデル：Flat-MPIとHybrid
08-APR22
27
Flat MPI vs. Hybrid
Flat-MPI：Each PE -> Independent
Memory
P
E
P
E
P
E
Memory
P
E
P
E
P
E
P
E
Hybrid：Hierarchal Structure
Memory
Memory
P P P P
E E E E
P P P P
E E E E
P
E
08-APR22
28
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における最
適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
08-APR22
29
科学技術計算における
大規模線形方程式の解法
• 多くの科学技術計算は，最終的に大規模線形方程式
Ax=bを解くことに帰着される。
– important, expensive
• アプリケーションに応じて様々な手法が提案されている
– 疎行列（sparse），密行列（dense）
– 直接法（direct），反復法（iterative）
• 密行列（dense）
– グローバルな相互作用：BEM，スペクトル法，MO，MD（気液）
• 疎行列（sparse）
– ローカルな相互作用：FEM，FDM，MD（固），高速多重極展開付
BEM
08-APR22
30
GeoFEMにおける戦略
• 反復解法のみが大規模並列計算のための唯一の手法
である。
• 前処理が重要
– ILU(0)/IC(0)（不完全LU分解，不完全コレスキー分解）などの
一般的手法によって広範囲の問題を解くことが可能である。
– 問題に特有の前処理手法が必要になる場合もある。
08-APR22
31
前処理（preconditioning）とは?
• 反復法の収束は係数行列の固有値分布に依存
– 固有値分布が少なく，かつ1に近いほど収束が早い（単位行列）
– 条件数（condition number）（対称正定）＝最大最小固有値の比
• 条件数が1に近いほど収束しやすい
08-APR22
32
前処理（preconditioning）とは?
共役勾配法の収束は係数行列の性質（固有値分布）に強く依存
するため，前処理を適用して固有値が1の周辺に集まるように
行列の性質を改善する。前処理行列を[M]とすると，
 
~
A  M
~
1
，
として，
  A  b   M  b 
A { }  {b }
1
~
~
という方程式を代わりに解くことになる。
[M]-1が[A]-1をよく近似した行列であれば，
 
~
A  M
 1  A 
は単位行列に近くなり，それだけ解きやすくなる。
08-APR22
33
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における
最適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
08-APR22
34
Block IC(0)-CG Solver
ブロックIC前処理付共役勾配法
• 三次元弾性問題
• GeoFEMにおける並列反復法ソルバー使用
– http://geofem.tokyo.rist.or.jp/
– 局所ブロックIC(0)前処理
• 正確には非対角成分を保存した不完全修正コレスキー分解
– Additive Schwartz Domain Decomposition（ASDD）
• Hybrid並列プログラミング
– OpenMP+MPI
– 1-SMPノード＝1領域
– 3層の並列構造
• ノード間：MPI，ノード内：OpenMP，各PE：ベクトル化
– ベクトル化，並列化のためのオーダリング手法
08-APR22
35
Flat MPI vs. Hybrid
Flat-MPI：各プロセッサ独立
Memory
P
E
P
E
P
E
Memory
P
E
P
E
P
E
P
E
Hybrid：階層構造
Memory
Memory
P P P P
E E E E
P P P P
E E E E
P
E
08-APR22
36
局所データ構造
節点ベース領域分割（領域間オーバーラップあり）
PE#1
21
22
PE#0
23
24
4
PE#1
5
17
12
11
18
13
7
6
19
２
PE#3
３
15
9
４
PE#2
6
7
PE#0
2
3
10
7
8
9
11
10
12
5
11
14
13
4
5
10
10
1
2
3
8
9
11
12
10
9
11
12
9
6
3
１
15
20
14
8
12
25
1
16
6
8
4
8
4
6
5
5
1
2
PE#3
7
7
1
2
PE#2
3
08-APR22
37
各SMPノード＝幾何学的な領域
MPIによる領域間通信
Node-0
Node-2
Node-1
Node-3
Node-0
PE
PE M
PE e
PE m
PE o
PE r
PE y
PE
Node-1
PE
PE M
PE e
PE m
PE o
PE r
PE y
PE
PE
PE M
PE e
PE m
PE o
PE r
PE y
PE
Node-2
PE
PE M
PE e
PE m
PE o
PE r
PE y
PE
Node-3
08-APR22
38
前処理付共役勾配法
Preconditioned Conjugate Gradient Method （PCG）
Compute r(0)= b-[A]x(0)
for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
ri-1= r(i-1) z(i-1)
if i=1
p(1)= z(0)
else
bi-1= ri-1/ri-2
p(i)= z(i-1) + bi-1 z(i)
endif
q(i)= [A]p(i)
ai = ri-1/p(i)q(i)
x(i)= x(i-1) + aip(i)
r(i)= r(i-1) - aiq(i)
check convergence |r|
end
08-APR22
39
前処理付き反復法の計算プロセス
•
•
•
•
内積
ベクトルの定数倍の加減
行列ベクトル積
前処理
08-APR22
40
前処理付き反復法の計算プロセス
「前処理」以外は容易にベクトル化可能
① ベクトル内積
RHO= 0.d0
do i= 1, N
RHO= RHO + r(i)*z(i)
enddo
② ベクトルの実数倍の加減
do i= 1, N
p(i)= z(i) + BETA*p(i)
enddo
③ 行列ベクトル積
do i= 1, N
VAL= D(i)*p(i)
do j= 1, INL(i)
VAL= VAL + AL(j,i)*p(IAL(j,i))
enddo
do j= 1, INU(i)
VAL= VAL + AU(j,i)*p(IAU(j,i))
enddo
q(i)= VAL
enddo
① ベクトル内積
RHO= 0.d0
!CDIR NODEP
do i= 1, N
RHO= RHO + r(i)*z(i)
enddo
② ベクトルの実数倍の加減
!CDIR NODEP
do i= 1, N
p(i)= z(i) + BETA*p(i)
enddo
③ 行列ベクトル積
do i= 1, N
VAL= D(i)*p(i)
!CDIR NODEP
do j= 1, INL(i)
VAL= VAL + AL(j,i)*p(IAL(j,i))
enddo
!CDIR NODEP
do j= 1, INU(i)
VAL= VAL + AU(j,i)*p(IAU(j,i))
enddo
q(i)= VAL
enddo
08-APR22
41
ベクトル計算
• ベクトルレジスタサイズが256の場合（e.g.地球シ
ミュレータ），下記において，N=1000とすると：
do i= 1, N
A(i)= B(i) + C(i)
enddo
• i=1~256のループは同時に計算される（並列性）。
– 以下同様に，i=257~512, i=513~768, i=769~1000も
同時に計算される。チューニングによって100倍変わる。
08-APR22
42
不完全修正コレスキー分解
i
l
ik
 d k  l jk  a ij
( j  1, 2 ,  , i  1)
k 1
i
l
ik
d k  l ik  a ii
k 1
ここで l ii  d i  1 とすると以下が導かれる
i  1, 2 ,  , n
j  1, 2 ,  , i  1
j 1
l ij  a ij 
l
ik
l ij  a ij
 d k  l jk
k 1

d i   a ii 


i 1

k 1
l ik
2

dk 


1
 l ii
1
対角成分のみがもとの
行列と変わる
08-APR22
43
不完全修正コレスキー分解を使用した
前進後退代入
 M z   LDL T z   r 
z   LDL
T
 r 
1
 L  y   r 
DL z   y 
T
d1

0

0

0
0

0
0
0
d2
0
0
0
d3
0
0
0
d4
0
0
0
0   l11

0
0

0  0

0  0
d 5   0
l 21
l 31
l 41
l 22
l 32
l 42
0
l 33
l 43
0
0
l 44
0
0
0
l 51   1
 
l 52
0
 
l 53    0
 
l 54   0
l 55   0
l 21 / l11
l 31 / l11
l 41 / l11
1
l 32 / l 22
l 42 / l 22
0
1
l 43 / l 33
0
0
1
0
0
0
l 51 / l11 

l 52 / l 22

l 53 / l 33 

l 54 / l 44 
1 
データ依存性：メモリの読込み書出し
が同時に発生し，並列化・ベクトル化困
難
不完全修正 do i= 1, N
08-APR22
44
コレスキー
分解
VAL= D(i)
do j= 1, INL(i)
VAL= VAL - (AL(j,i)**2) * W(IAL(j,i),DD)
enddo
W(i,DD)= 1.d0/VAL
enddo
前進代入
do i= 1, N
WVAL= W(i,Z)
do j= 1, INL(i)
WVAL= WVAL - AL(j,i) * W(IAL(j,i),Z)
enddo
W(i,Z)= WVAL * W(i,DD)
enddo
08-APR22
45
「地球シミュレータ」における最適化
• 仮説
– IC(0)前処理はFactorizationを含むため，何らかの陽的なデー
タ操作（e.g. Re-Ordering）をしないと，ベクトル化はできないだ
ろう（データ依存性の排除によるベクトル処理，最適化）。
• 並列・ベクトル計算のためのRe-Ordering手法
– 処理の局所化，依存性の排除
– 連続的なメモリアクセス
– 長いループ長（ベクトル化の場合）
08-APR22
46
要素の色づけ（coloring）
• 例えば以下のようなオペレーション（Gauss-Seidel）
を考慮した場合
do j= 1, N
do i= 1, M
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
(i,j+1)
(i-1,j)
(i,j)
(i,j-1)
(i+1,j)
08-APR22
47
要素の色づけ（coloring）
• 要素（i,j）と隣接する4要素はデータ依存性がある
（dependent）：このままではベクトル化できない
do j= 1, N
do i= 1, M
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
(i,j+1)
(i-1,j)
(i,j)
(i,j-1)
(i+1,j)
08-APR22
48
要素の色づけ（coloring）
• 無理にベクトル化すると正しい解が得られない。
do j= 1, N
!CDIR NODEP
do i= 1, M
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
(i,j+1)
(i-1,j)
(i,j)
(i,j-1)
(i+1,j)
08-APR22
49
要素の色づけ（coloring）
• 要素(i,j) と (i-1,j+1), (i+1,j+1), (i-1,j-1), (i+1,j+1) はデー
タ依存性が無い（independent）。
• したがって，これらの要素（下図の黄緑色の要素）は
同時に計算が可能である：並列性⇒ベクトル計算の可
能性
(i-1,j+1)
(i,j+1)
(i+1,j+1)
(i-1,j)
(i,j)
(i+1,j)
(i-1,j-1)
(i,j-1)
(i+1,j-1)
08-APR22
50
Red-Black Coloring
08-APR22
51
Red-Black Coloring
赤と黒：独立に計算できる
08-APR22
52
Red-Black Coloring
do j= 1, N
do i= 1, M
if (RED(i,j) then
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
endif
enddo
enddo
do j= 1, N
do i= 1, M
if (BLACK(i,j) then
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
endif
enddo
enddo
08-APR22
53
Red-Black Coloring
RED（i+jが偶数）：右辺はBLACK
do j= 1, N-1, 2
!CDIR NODEP
do i= 1, M-1, 2
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
do j= 2, N, 2
!CDIR NODEP
do i= 2, M, 2
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
08-APR22
54
Red-Black Coloring
BLACK（i+jが奇数）：右辺はRED
do j= 1, N-1, 2
!CDIR NODEP
do i= 2, M, 2
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
do j= 2, N, 2
!CDIR NODEP
do i= 1, M-1, 2
PHI(i,j)= (RHS(i,j)&
AW(i,j)*PHI(i-1,j)-AE(i,j)*PHI(i+1,j)&
AS(i,j)*PHI(i,j-1)-AN(i,j)*PHI(i,j+1))/AC(i,j
enddo
enddo
08-APR22
55
Red-Black Coloring
• Red要素を計算する場合
– 右辺に来るのはBlack要素
– Red要素を計算している間，
右辺のBlack要素は不変
• Black要素を計算する場合
– 右辺に来るのはRed要素
– Black要素を計算している間，
右辺のRed要素は不変
• データ依存性を排除するこ
とができる
08-APR22
56
Red-Black Ordering/Reordering
13 14 15 16
15
7
16
8
9
10 11 12
5
13
6
14
5
6
7
8
11
3
12
4
1
2
3
4
1
9
2
10
• 差分格子のような（i,j）型ではなく，1から順番に要
素番号が振られている場合，要素番号を「Red」⇒
「Black」の順にふり直す（reordering, ordering）と
プログラムが簡単になる（計算効率も高い）
08-APR22
57
ベクトル計算のためのオーダリング法
15 11 16 12
15
7
16
8
10 13 15 16
9
13 10 14
5
13
6
14
6
9
12 14
7
3
8
4
11
3
12
4
3
5
8
11
1
5
２
6
1
9
2
10
1
2
4
7
マルチカラーオーダリング
（4色，2色＝Red-Black）：
- 要素間結合グラフをもとにオーダリング
- 並列性，独立性重視
CM法（Cuthill-McKee）：
- 要素間結合グラフに基づき
オーダリング
- 色間の依存性も考慮
• いずれの手法も，要素間結合グラフを利用。
• 同じ「色」付けされた要素同士は，データ依存性が無い。
– 並列計算（すなわちベクトル計算）が可能
08-APR22
58
Re-Ordering Technique for
Vector/Parallel Architectures
Cyclic DJDS(RCM+CMC) Re-Ordering
(Doi, Washio, Osoda and Maruyama (NEC))
1. RCM (Reverse Cuthil-Mckee)
2. CMC (Cyclic Multicolor)
Vector
SMP Parallel
3. DJDS re-ordering
Vector
4. Cyclic DJDS for SMP unit
SMP Parallel
These processes can be substituted by
traditional multi-coloring (MC).
08-APR22
59
「地球シミュレータ」における最適化
• 仮説
– IC(0)前処理はFactorizationを含むため，何らかの陽的なデー
タ操作（e.g. Re-Ordering）をしないと，ベクトル化はできないだ
ろう（データ依存性の排除によるベクトル処理，最適化）。
• 並列・ベクトル計算のためのRe-Ordering手法
– 処理の局所化，依存性の排除
– 連続的なメモリアクセス
– 長いループ長（ベクトル化の場合）
この2項目はクリア
08-APR22
60
「地球シミュレータ」における最適化
• 仮説
– IC(0)前処理はFactorizationを含むため，何らかの陽的なデー
タ操作（e.g. Re-Ordering）をしないと，ベクトル化はできないだ
ろう（データ依存性の排除によるベクトル処理，最適化）。
• 並列・ベクトル計算のためのRe-Ordering手法
– 処理の局所化，依存性の排除
– 連続的なメモリアクセス
– 長いループ長（ベクトル化の場合）
次はこれ
08-APR22
61
非構造格子における疎行列格納法
1D-Storage (CRS)
memory saved, short vector length
2D-Storage
long vector length, many ZERO’s
08-APR22
62
各「色」内におけるRe-Ordering
非ゼロ非対角成分の数
各「色」内の要素は互いに独立，従って，色内で計算順序を
変更しても計算結果は変わらない。
DJDS : Descending-Order Jagged Diagonal Storage
08-APR22
63
「地球シミュレータ」における最適化
• 仮説
– IC(0)前処理はFactorizationを含むため，何らかの陽的なデー
タ操作（e.g. Re-Ordering）をしないと，ベクトル化はできないだ
ろう（データ依存性の排除によるベクトル処理，最適化）。
• 並列・ベクトル計算のためのRe-Ordering手法
– 処理の局所化，依存性の排除
– 連続的なメモリアクセス
– 長いループ長（ベクトル化の場合）
• Hybrid Parallel Programming：3層の並列構造
– ノード間 : MPI
– ノード内 : OpenMP
– 各PE
: Vectorization
08-APR22
64
Cyclic DJDS (MC/CM-RCM)
Cyclic Re-Ordering for SMP units
SMPノード内の各PEにおける負荷分散
iS+1
iv0+1
iE
do iv= 1, NCOLORS
!$omp parallel do
do ip= 1, PEsmpTOT
iv0= STACKmc(PEsmpTOT*(iv-1)+ip- 1)
do j= 1, NLhyp(iv)
iS= INL(npLX1*(iv-1)+PEsmpTOT*(j-1)+ip-1)
iE= INL(npLX1*(iv-1)+PEsmpTOT*(j-1)+ip )
!cdir nodep
do i= iv0+1, iv0+iE-iS
k= i+iS - iv0
kk= IAL(k)
(Important Computations)
enddo
enddo
enddo
enddo
npLX1= NLmax * PEsmpTOT
INL(0:NLmax*PEsmpTOT*NCOLORS)
08-APR22
65
Flat MPI & Hybridの違い
• ほとんどの労力はベクトル化に費やされる。
• 充分に長い依存性のないベクトルを生成できれば，そ
れを分割してSMPノード上の各PEにばらまくだけ。
• 従って，HybridとFlat MPIはプログラム上はほとんど変
りなくなる。
–Flat MPIはスレッド数＝１のHybridに相当する。
08-APR22
66
Cyclic DJDS (MC/CM-RCM)
for Forward/Backward Substitution in BILU Factorization
do iv= 1, NCOLORS
!$omp parallel do private (iv0,j,iS,iE… etc.)
do ip= 1, PEsmpTOT
iv0= STACKmc(PEsmpTOT*(iv-1)+ip- 1)
SMP
do j= 1, NLhyp(iv)
parallel
iS= INL(npLX1*(iv-1)+PEsmpTOT*(j-1)+ip-1)
iE= INL(npLX1*(iv-1)+PEsmpTOT*(j-1)+ip )
!CDIR NODEP
do i= iv0+1, iv0+iE-iS
k= i+iS - iv0
kk= IAL(k)
Vectorized
X(i)= X(i) - A(k)*X(kk)*DINV(i) etc.
enddo
enddo
enddo
enddo
08-APR22
67
Simple 3D Cubic Model
z
Uniform Distributed Force in
z-dirrection @ z=Zmin
Uy=0 @ y=Ymin
Ux=0 @ x=Xmin
Nz-1
Ny-1
Uz=0 @ z=Zmin
x
y
Nx-1
08-APR22
68
Effect of Ordering
08-APR22
69
Effect of Re-Ordering
Long Loops
Continuous Access
PDJDS/CM-RCM
Short Loops
Continuous Access
Short Loops
Irregular Access
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
08-APR22
70
Matrix Storage, Loops
• DJDS (Descending order
Jagged Diagonal Storage) は
最内ループ長が長くベクトル計
算向け
DJDS
do iv= 1, NVECT
iv0= STACKmc(iv-1)
do j= 1, NLhyp(iv)
iS= index_L(NL*(iv-1)+ j-1)
iE= index_L(NL*(iv-1)+ j )
do i= iv0+1, iv0+iE-iS
k= i+iS - iv0
kk= item_L(k)
Z(i)= Z(i) - AL(k)*Z(kk)
enddo
iS= STACKmc(iv-1) + 1
iE= STACKmc(iv )
do i= iS, iE
Z(i)= Z(i)/DD(i)
enddo
enddo
enddo
• DCRS （Compressed RowStorage）は，ローカルなオペ
レーションが特徴，キャッシュの
あるスカラープロセッサ向け
DCRS
do i= 1, N
SW= WW(i,Z)
isL= index_L(i-1)+1
ieL= index_L(i)
do j= isL, ieL
k = item_L(j)
SW= SW - AL(j)*Z(k)
enddo
Z(i)= SW/DD(i)
enddo
08-APR22
71
PDJDS/CM-RCM
Effect of Re-Ordering
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
Results on 1 SMP node
Color #: 99 (fixed)
Re-Ordering is REALLY required !!!
●
■
▲
1.E+02
GFLOPS
1.E+01
Effect of Vector Length
X10
1.E+00
+ Re-Ordering X100
22 GFLOPS,
34% of the Peak
1.E-01
1.E-02
1.E+04
1.E+05
1.E+06
DOF
Ideal Performance: 40%-45%
1.E+07 for Single CPU
08-APR22
72
Effect of Re-Ordering
Results on 1 SMP node
Color #: 99 (fixed)
Re-Ordering is REALLY required !!!
PDJDS/CM-RCM
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
●
■
▲
1.E+02
80x80x80 case (1.5M DOF)
GFLOPS
1.E+01
● 212 iter’s, 11.2 sec.
■ 212 iter’s, 143.6 sec.
▲ 203 iter’s, 674.2 sec.
1.E+00
1.E-01
1.E-02
1.E+04
1.E+05
1.E+06
DOF
1.E+07
08-APR22
PDJDS/CM-RCM
73
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
3D Elastic Simulation
Problem Size~GFLOPS
Earth Simulator
1 SMP node (8 PE’s)
■
▲
1.0E+02
1.0E+02
1.0E+01
1.0E+01
GFLOPS
GFLOPS
●
1.0E+00
1.0E+00
1.0E-01
1.0E+04
1.0E-01
1.0E+05
1.0E+06
1.0E+07
DOF: Problem Size
1.0E+04
1.0E+05
Flat-MPI is better.
Nice Intra-Node MPI.
1.0E+06
1.0E+07
DOF: Problem Size
Flat-MPI
Hybrid (OpenMP)
23.4 GFLOPS, 36.6 % of Peak
21.9 GFLOPS, 34.3 % of Peak
08-APR22
74
Hitachi
Earth
SR8000/MPP
Simulator
(U.Tokyo)
IBM SP3
(LBNL)
Hitachi
SR11000/J2
(U.Tokyo)
Peak Performance
GFLOPS
8.00
1.80
1.50
9.20
Measured Memory
BW (STREAM)
(GB/sec/PE)
26.6
2.85
0.623
8.00
Estimated
Performance/Core
GFLOPS
(% of peak)
2.31-3.24
(28.8-40.5)
.291-.347
(16.1-19.3)
.072-.076
(4.8-5.1)
.880-.973
(9.6-10.6)
Measured
Performance/Core
2.93
(36.6)
.335
(18.6)
.122
(8.1)
1.34
(14.5)
BYTE/FLOP
3.325
1.583
0.413
0.870
12.3
1.60
1.00
12.0
5.6-7.7
6-20
16.3
4.7 *
Comm. BW
(GB/sec/Node)
MPI Latency
(msec)
* IBM p595
J.T.Carter
et al.
08-APR22
75
3D Elastic Simulation
Problem Size~GFLOPS
PDJDS/CM-RCM
Hitachi-SR8000-MPP with
Pseudo Vectorization
Hybrid is better.
1 SMP node (8 PE’s)
Low Intra-Node MPI.
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
■
▲
3.00
3.00
2.50
2.50
2.00
2.00
GFLOPS
GFLOPS
●
1.50
1.00
1.50
1.00
0.50
0.50
0.00
0.00
1.0E+04
1.0E+05
1.0E+06
1.0E+07
DOF: Problem Size
1.0E+04
1.0E+05
1.0E+06
1.0E+07
DOF: Problem Size
Flat-MPI
Hybrid (OpenMP)
2.17 GFLOPS, 15.0 % of Peak
2.68 GFLOPS, 18.6 % of Peak
08-APR22
PDJDS/CM-RCM
76
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
3D Elastic Simulation
Problem Size~GFLOPS
IBM-SP3 (NERSC)
1 SMP node (8 PE’s)
■
▲
3.00
3.00
2.50
2.50
2.00
2.00
GFLOPS
GFLOPS
●
1.50
1.00
1.50
1.00
0.50
0.50
0.00
0.00
1.0E+04
Cache is well-utilized
in Flat-MPI.
1.0E+05
1.0E+06
DOF: Problem Size
Flat-MPI
1.0E+07
1.0E+04
1.0E+05
1.0E+06
DOF: Problem Size
Hybrid (OpenMP)
1.0E+07
08-APR22
PDJDS/CM-RCM
77
PDCRS/CM-RCM CRS no re-ordering
short innermost loop
3D Elastic Simulation
Problem Size~GFLOPS
Hitachi SR11000/J2 (U.Tokyo)
1 SMP node (8 PE’s)
■
▲
15.0
15.0
10.0
10.0
GFLOPS
GFLOPS
●
5.0
5.0
0.0
0.0
1.E+04
1.E+05
1.E+06
DOF
Flat-MPI
1.E+07
1.E+04
1.E+05
1.E+06
DOF
Hybrid (OpenMP)
1.E+07
08-APR22
78
SMP node # > 10
up to 176 nodes (1408 PEs)
Problem size for each SMP node is fixed.
PDJDS-CM/RCM, Color #: 99
08-APR22
79
3D Elastic Model (Large Case)
256x128x128/SMP node, up to 2,214,592,512 DOF
Parallel Work Ratio
GFLOPS rate
4000
100
GFLOPS
3000
●Flat MPI
2000
1000
0
Parallel Work Ratio: %
○Hybrid
90
80
70
60
50
40
0
16
32
48
64
80
96 112 128 144 160 176 192
NODE#
0
16
32
48
64
80
96 112 128 144 160 176 192
NODE#
3.8TFLOPS
for 2.2G DOF
●：Flat MPI, ○：Hybrid
176 nodes (33.8% of peak)
08-APR22
80
3D Elastic Model (Small Case)
64x64x64/SMP node, up to 125,829,120 DOF
Parallel Work Ratio
GFLOPS rate
100
GFLOPS
3000
○Hybrid
2000
1000
●Flat MPI
0
Parallel Work Ratio: %
4000
90
80
70
60
50
40
0
16
32
48
64
80
96 112 128 144 160 176 192
0
16
32
48
64
80
96
112 128 144 160 176 192
NODE#
NODE#
●：Flat MPI, ○：Hybrid
08-APR22
81
Hybrid vs. Flat-MPI
• Hybridが優位に立つのは・・・
– SMPノードの数が多い場合
– 各ノードあたりの問題規模が小さい場合
• なぜなら，Flat-MPIは
– 通信プロセス数が8倍
– 通信/計算 比が2倍（通信の割合が高い）
• 通信の影響はノード数が多くなると顕著になる
• D.Kerbyson (LANL)による性能予測
– LA-UR-02-5222
– ESの通信 Latency の影響
08-APR22
Flat-MPI and Hybrid
Flat-MPI
Hybrid
N
N
N
82
Flat MPI
Hybrid
Problem size/each MPI Process
(N=number of FEM nodes in one
direction of cube geometry
3N3
38N3
Size of messages on each
surfaces with each neighboring
domain
3N2
34N2
Ratio of
communication/computation
1/N
1/(2N)
08-APR22
83
Hitachi
Earth
SR8000/MPP
Simulator
(U.Tokyo)
IBM SP3
(LBNL)
Hitachi
SR11000/J2
(U.Tokyo)
Peak Performance
GFLOPS
8.00
1.80
1.50
9.20
Measured Memory
BW (STREAM)
(GB/sec/PE)
26.6
2.85
0.623
8.00
Estimated
Performance/Core
GFLOPS
(% of peak)
2.31-3.24
(28.8-40.5)
.291-.347
(16.1-19.3)
.072-.076
(4.8-5.1)
.880-.973
(9.6-10.6)
Measured
Performance/Core
2.93
(36.6)
.335
(18.6)
.122
(8.1)
1.34
(14.5)
BYTE/FLOP
3.325
1.583
0.413
0.870
12.3
1.60
1.00
12.0
5.6-7.7
6-20
16.3
4.7 *
Comm. BW
(GB/sec/Node)
MPI Latency
(msec)
* IBM p595
J.T.Carter
et al.
08-APR22
84
通信オーバーヘッドは何故生じる?
• 立ち上がりのLatency
• 有限の通信バンド幅
• 同期オーバーヘッド
– SEND/RECV, ALLREDUCE etc.
• 実はメモリバンド幅も影響
– 送受信の際にメモリーコピーが発生する
08-APR22
85
Domain-to-Domain Communication
Exchange Boundary Information (SEND/RECV)
subroutine SOLVER_SEND_RECV
(N, NEIBPETOT, NEIBPE,
IMPORT_INDEX, IMPORT_NODE,
EXPORT_INDEX, EXPORT_NODE,
WS, WR, X, SOLVER_COMM, my_rank)
&
&
&
&
implicit REAL*8 (A-H,O-Z)
include 'mpif.h'
parameter (KREAL= 8)
integer IMPORT_INDEX(0:NEIBPETOT), IMPORT_NODE(N)
integer EXPORT_INDEX(0:NEIBPETOT), EXPORT_NODE(N)
integer SOLVER_COMM, my_rank
integer req1(NEIBPETOT), req2(NEIBPETOT)
integer sta1(MPI_STATUS_SIZE, NEIBPETOT)
integer sta2(MPI_STATUS_SIZE, NEIBPETOT)
real(kind=KREAL) X(N), NEIBPE(NEIBPETOT), WS(N), WR(N)
do neib= 1, NEIBPETOT
SEND
istart= EXPORT_INDEX(neib-1)
inum = EXPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
WS(k)= X(EXPORT_NODE(k))
enddo
call MPI_ISEND
(WS(istart+1), inum, MPI_DOUBLE_PRECISION,
&
NEIBPE(neib), 0, SOLVER_COMM,
&
req1(neib), ierr)
enddo
do neib= 1, NEIBPETOT
RECEIVE
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
call MPI_IRECV
(WR(istart+1), inum, MPI_DOUBLE_PRECISION, &
NEIBPE(neib), 0, SOLVER_COMM,
&
req2(neib), ierr)
enddo
call MPI_WAITALL (NEIBPETOT, req2, sta2, ierr)
do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
X(IMPORT_NODE(k))= WR(k)
enddo
enddo
call MPI_WAITALL (NEIBPETOT, req1, sta1, ierr)
return
end
08-APR22
86
Domain-to-Domain Communication
Exchange Boundary Information (SEND/RECV)
subroutine SOLVER_SEND_RECV
(N, NEIBPETOT, NEIBPE,
IMPORT_INDEX, IMPORT_NODE,
EXPORT_INDEX, EXPORT_NODE,
WS, WR, X, SOLVER_COMM, my_rank)
&
&
&
&
implicit REAL*8 (A-H,O-Z)
include 'mpif.h'
parameter (KREAL= 8)
integer IMPORT_INDEX(0:NEIBPETOT), IMPORT_NODE(N)
integer EXPORT_INDEX(0:NEIBPETOT), EXPORT_NODE(N)
integer SOLVER_COMM, my_rank
integer req1(NEIBPETOT), req2(NEIBPETOT)
integer sta1(MPI_STATUS_SIZE, NEIBPETOT)
integer sta2(MPI_STATUS_SIZE, NEIBPETOT)
real(kind=KREAL) X(N), NEIBPE(NEIBPETOT), WS(N), WR(N)
do neib= 1, NEIBPETOT
SEND
istart= EXPORT_INDEX(neib-1)
inum = EXPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
WS(k)= X(EXPORT_NODE(k))
enddo
call MPI_ISEND
(WS(istart+1), inum, MPI_DOUBLE_PRECISION,
&
NEIBPE(neib), 0, SOLVER_COMM,
&
req1(neib), ierr)
enddo
do neib= 1, NEIBPETOT
RECEIVE
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
call MPI_IRECV
(WR(istart+1), inum, MPI_DOUBLE_PRECISION, &
NEIBPE(neib), 0, SOLVER_COMM,
&
req2(neib), ierr)
enddo
call MPI_WAITALL (NEIBPETOT, req2, sta2, ierr)
do neib= 1, NEIBPETOT
istart= IMPORT_INDEX(neib-1)
inum = IMPORT_INDEX(neib ) - istart
do k= istart+1, istart+inum
X(IMPORT_NODE(k))= WR(k)
enddo
enddo
call MPI_WAITALL (NEIBPETOT, req1, sta1, ierr)
return
end
08-APR22
87
Communication Overhead
Memory
Copy
Comm.
Latency
Comm.
BandWidth
08-APR22
88
Communication Overhead
Memory
Copy
depends on
message size
Comm.
Latency
Comm.
BandWidth
depends on
message size
08-APR22
89
Communication Overhead
Earth Simulator
Comm.
Latency
08-APR22
90
Communication Overhead
Hitachi SR11000, IBM SP3 etc.
Memory
Copy
Comm.
BandWidth
08-APR22
91
Communication Overhead
= Synchronization Overhead
08-APR22
92
Communication Overhead
= Synchronization Overhead
Memory
Copy
Comm.
Latency
Comm.
BWTH
08-APR22
93
Communication Overhead
= Synchronization Overhead
Earth Simulator
Comm.
Latency
08-APR22
94
Communication Overhead
Weak Scaling: Earth Simulator
Comm. Overhead (sec.)
0.06
●○ 3x503 DOF/PE
▲△ 3x323 DOF/PE
●▲ Flat-MPI ○△ Hybrid
0.04
Effect of message size
is small. Effect of latency
is large.
0.02
0.00
10
100
1000
PE#
10000
Memory-copy is so fast.
08-APR22
95
Communication Overhead
Weak Scaling: IBM SP-3
Comm. Overhead (sec.)
0.40
●○ 3x503 DOF/PE
▲△ 3x323 DOF/PE
●▲ Flat-MPI ○△ Hybrid
0.30
Effect of message size
is more significant.
0.20
0.10
0.00
10
1000
100
PE#
10000
08-APR22
96
Communication Overhead
Weak Scaling:
Hitachi SR11000/J2 (8cores/node)
Comm. Overhead (sec.)
0.10
●○ 3x503 DOF/PE
▲△ 3x323 DOF/PE
●▲ Flat-MPI ○△ Hybrid
0.08
0.06
0.04
0.02
0.00
10
100
cores
1000
08-APR22
97
まとめ
• 「地球シミュレータ」における非構造格子向け
Hybridプログラミングモデルを使用した前処理付き
反復法を開発した。
• 高い並列，ベクトル性能を達成した。22億自由度の
三次元弾性問題，BIC(0)-CG法により176ノード
（1408 PE）を使用して3.8TFLOPS（ピーク性能
の33.8%）を達成した。
– 櫛田氏 (奥田研学生，現JAEA) は512ノードを使用して30
億自由度の問題で10TFLOPSを達成した。
• Re-Orderingはやっぱり必要である。
08-APR22
98
まとめ（Flat MPI vs. Hybrid）
• 少ないノード数ではFlat MPI有利。
• 大きなノード数ではHybrid有利。
– 特に1ノードあたりの問題規模が小さい場合
• 言い換えると・・・バンド幅の問題
– 通信支配：Hybrid
– メモリ支配：Flat MPI
– これはHW条件，問題規模ばかりでなく，アプリケーションの傾
向でもある。
• Hybrid はFlat-MPIと比較して色数に対して敏感
– Hybrid is much more sensitive to color numbers than flat
MPI due to synchronization overhead of OpenMP.
• 地球シミュレータのノード間ネットワーク
– 今にして有難さを実感している
08-APR22
99
• 背景
–
–
–
–
スカラープロセッサとベクトルプロセッサ
GeoFEMプロジェクト
地球シミュレータ
前処理付き反復法
• FEMコードの「地球シミュレータ」における
最適化戦略
– BIC(0)-CG法による並列反復法ソルバー
– 係数行列作成部分
• まとめ
08-APR22
100
“CUBE” Benchmark
• 3D linear elastic applications on cubes for a wide
range of problem size.
• 1CPUによる比較
– Earth Simulator
– AMD Opteron (1.8GHz)
z
Uniform Distributed Force in
z-dirrection @ z=Zmin
Uy=0 @ y=Ymin
Ux=0 @ x=Xmin
Nz-1
Ny-1
Uz=0 @ z=Zmin
x
y
Nx-1
08-APR22
101
Time for 3x643=786,432 DOF
ES
8.0 GFLOPS
Opteron
1.8GHz
DCRS
sec.
(MFLOPS)
DJDS original
sec.
(MFLOPS)
Matrix
28.6
(291)
34.2
(240)
Solver
360
(171)
21.7
(3246)
Total
389
55.9
Matrix
10.2
(818)
12.4
(663)
Solver
225
(275)
271
(260)
Total
235
283
3.6 GFLOPS
DCRS
DJDS
08-APR22
102
Matrix+Solver
60
50
sec.
40
DJDS on ES
original
30
20
10
0
41472
98304
192000
375000
786432
DOF
300
sec.
200
DJDS on Opteron
original
100
0
41472
98304
192000
DOF
375000
786432
Matrix
Solver
08-APR22
103
Computation Time vs. Problem Size
40
300
Matrix
Solver
30
sec.
sec.
200
20
100
10
0
0
41472
98304
192000
375000
786432
DOF
41472
98304
192000
375000
786432
DOF
300
Total
ES (DJDS original)
sec.
200
Opteron
(DJDS original)
Opteron (DCRS)
100
0
41472
98304
192000
DOF
375000
786432
08-APR22
104
係数マトリクス生成部のほうが実は計
算時間を要している
• この部分も「ベクトル化」の必要がある
• 非線形計算（例えば弾塑性解析，カップルした
Navier-Stokes方程式）では，非線形反復のたびに
係数マトリクスを更新する必要がある。
• 係数マトリクス生成はアプリケーションに強く依
存するため，線形ソルバーなどの違って「ライブ
ラリ化」もしにくい
– 特にFEMの場合，ベクトル化が困難な複雑なプロセス
を含む可能性もある
08-APR22
105
FEMにおける係数マトリクスの計算
• ガラーキン法を各要素に適用して積分
– 要素マトリクスを得る
• 要素マトリクスを重ね合わせて「全体マトリクス
（Global Matrix）」を得る
• 要素マトリクスの計算そのものは，要素ごとに独
立して実施できるため，いわゆるEmbarassingly
Parallelなプロセスである
08-APR22
106
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
• Linear equations for each node
19
20
21
22
23
24
13
14
15
16
17
18
7
8
9
10
11
12
1
2
3
4
5
6
Elements
08-APR22
107
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
• Linear equations for each node
29
30
19
22
20
23
13
15
24
16
8
21
12
13
5
5
28
20
12
Nodes
18
11
4
4
27
19
11
35
24
17
10
3
3
26
18
10
34
23
16
9
2
2
25
17
9
33
22
15
8
1
32
21
14
7
1
31
14
6
6
7
08-APR22
108
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
• Linear equations for each node
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
X X
08-APR22
109
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
E
E
E
E
• Linear equations for each node
11

E
 21
 E 31

 E 41
12
13
E 22
E 23
E 32
E 33
E 42
E 43


E 24

E 34 

E 44 
14
08-APR22
110
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
• Linear equations for each node
29
30
31
32
33
34
35
22
23
24
25
26
27
28
15
16
17
18
19
20
21
8
9
10
11
12
13
14
1
2
3
4
5
6
7
08-APR22
111
Element-by-Element Operations
• Integration over each element => element-matrix
• Element matrices are accumulated to each node
=> global-matrix
• Linear equations for each node
 a 1 ,1

a
 2 ,1








a1, 2
a 2,2
a 3,3
...
a 33 , 33
a 34 , 34
a 35 , 34
  u1   f1 

 

u2
f2

 

 u3   f3 

 

  ...    ... 
 u   f 
33
  33 

a 34 , 35   u 34   f 34 

 

a 35 , 35   u 35   f 35 
08-APR22
112
Element-by-Element Operations
22
23
13
15
14
16
7
8
24
17
8
9
10
 a 1 ,1

a
 2 ,1








a1, 2
a 2,2
a 3,3
...
a 33 , 33
a 34 , 34
a 35 , 34
  u1   f1 

 

u2
f2

 

 u3   f3 

 


...
...





 u   f 
33
  33 

a 34 , 35   u 34   f 34 

 

a 35 , 35   u 35   f 35 
• a23,16 と a16,23を計算するためには，13番と14番
要素の寄与を考慮する必要がある
08-APR22
113
Current Approach
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4
do je= 1, 4
- assemble element-matrix
- accumulated element-matrix to global-matrix
enddo
enddo
enddo
08-APR22
114
Current Approach
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4 Local Node ID
do je= 1, 4
- assemble element-matrix
- accumulat element-matrix to global-matrix
enddo
enddo
enddo
4
3
Local Node ID
for each bi-linear 4-node
element
1
2
08-APR22
115
Current Approach
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4 Local Node ID
do je= 1, 4
- assemble element-matrix
- accumulat element-matrix to global-matrix
enddo
enddo
enddo
• スカラープロセッサ向け
– キャッシュの再利用に適している
• ベクトルプロセッサ向けで無い
– a16,23 ・ a23,16 ⇒正しく計算されな
い可能性がある
– 最内ループ短い
– 分岐が多い
ES
8.0 GFLOPS
Opteron
1.8GHz
DCRS
sec.
(MFLOPS)
DJDS original
sec.
(MFLOPS)
Matrix
28.6
(291)
34.2
(240)
Solver
360
(171)
21.7
(3246)
Total
389
55.9
Matrix
10.2
(818)
12.4
(663)
Solver
225
(275)
271
(260)
Total
235
283
3.6 GFLOPS
08-APR22
116
Inside the loop: integration at
Gaussian quadrature points
do jpn= 1, 2
do ipn= 1, 2
coef= dabs(DETJ(ipn,jpn))*WEI(ipn)*WEI(jpn)
PNXi= PNX(ipn,jpn,ie)
PNYi= PNY(ipn,jpn,ie)
PNXj= PNX(ipn,jpn,je)
PNYj= PNY(ipn,jpn,je)
a11=
a22=
a12=
a21=
enddo
enddo
a11
a22
a12
a21
+
+
+
+
(valX*PNXi*PNXj
(valX*PNYi*PNYj
(valA*PNXi*PNYj
(valA*PNYi*PNXj
+
+
+
+
valB*PNYi*PNYj)*coef
valB*PNXi*PNXj)*coef
valB*PNXj*PNYi)*coef
valB*PNYj*PNXi)*coef
08-APR22
117
ベクトルプロセッサ向け対策
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4
do je= 1, 4
- assemble element-matrix
- accumulat element-matrix to global-matrix
enddo
enddo
enddo
• a16,23 ・ a23,16 ⇒（無理にベクトル化すると）正
しく計算されない可能性あり
– 再び「色づけ（coloring）」
– 同じ節点を共有する要素は違う「色」に属する
08-APR22
118
Coloring of Elements
29
30
31
32
33
34
35
22
23
24
25
26
27
28
15
16
17
18
19
20
21
8
9
10
11
12
13
14
1
2
3
4
5
6
7
08-APR22
119
Coloring of Elements
Elements sharing the 16th node are assigned to different colors
29
30
31
32
33
34
35
22
23
24
25
26
27
28
15
16
17
18
19
20
21
8
9
10
11
12
13
14
1
2
3
4
5
6
7
08-APR22
120
ベクトルプロセッサ向け対策
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4
do je= 1, 4
- assemble element-matrix
- accumulat element-matrix to global-matrix
enddo
enddo
enddo
• a16,23 ・ a23,16 ⇒正しく計算されない可能性あり
– 再び「色づけ（coloring）」
– 同じ節点を共有する要素は違う「色」に属する
• 短い最内ループ
– ループ交換
08-APR22
121
ベクトルプロセッサ向け対策
22
23
13
15
14
16
7
8
24
17
8
9
10
do icel= 1, ICELTOT
do ie= 1, 4
do je= 1, 4
- assemble element-matrix
- accumulat element-matrix to global-matrix
enddo
enddo
enddo
• a16,23 ・ a23,16 ⇒正しく計算されない可能性あり
– 再び「色づけ（coloring）」
– 同じ節点を共有する要素は違う「色」に属する
• 短い最内ループ
– ループ交換
• 分岐が多い
– 要素～全体行列参照配列（Element-to-Matrix array）
08-APR22
122
Define ELEMENT-to-MATRIX array
22
23
13
15
24
7
23
④
15
①
③
23
13
17
8
9
③
22
14
16
8
④
24
14
16
16
②
①
Local Node ID
17
②
10
if
if
kkU= index_U(16-1+k) and
item_U(kkU)= 23
then
ELEMmat(13,2,3)= +kkU
ELEMmat(14,1,4)= +kkU
endif
ELEMmat(icel, ie, je)
Element ID
Local Node ID
kkL= index_L(23-1+k) and
item_L(kkL)= 16
then
ELEMmat(13,3,2)= -kkL
ELEMmat(14,4,1)= -kkL
endif
08-APR22
123
Define ELEMENT-to-MATRIX array
22
23
13
15
24
7
23
④
15
①
③
23
13
17
8
9
③
22
14
16
8
④
24
14
16
16
②
①
Local Node ID
17
②
10
if
if
kkU= index_U(16-1+k) and
item_U(kkU)= 23
then
ELEMmat(13,2,3)= +kkU
ELEMmat(14,1,4)= +kkU
endif
kkL= index_L(23-1+k) and
item_L(kkL)= 16
then
ELEMmat(13,3,2)= -kkL
ELEMmat(14,4,1)= -kkL
endif
“ELEMmat”は各要素の各節点と全体マトリクスにおけるアド
レスとの関係を示す
08-APR22
124
Optimized Procedure
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- define “ELEMmat” array
enddo
enddo
enddo
enddo
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- assemble element-matrix
enddo
enddo
enddo
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- accumulate element-matrix to global-matrix
enddo
enddo
enddo
enddo
Extra Computation for
• ELEMmat
Extra Storage for
• ELEMmat array
• element-matrix components
for elements in each color
• < 10% increase
08-APR22
125
Optimized Procedure
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- define “ELEMmat” array
enddo
enddo
enddo
enddo
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- assemble element-matrix
enddo
enddo
enddo
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- accumulate element-matrix to global-matrix
enddo
enddo
enddo
enddo
PART I
“Integer” operations for
“ELEMmat”
In nonlinear cases, this part
should be done just once
(before initial iteration), as long
as mesh connectivity does not
change.
08-APR22
126
Optimized Procedure
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- define “ELEMmat” array
enddo
enddo
enddo
enddo
do icol= 1, NCOLOR_E_tot
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- assemble element-matrix
enddo
enddo
enddo
do ie= 1, 4
do je= 1, 4
do ic0= index_COL(icol-1)+1, indexCOL_(icol)
icel= item_COL(ic0)
- accumulate element-matrix to global-matrix
enddo
enddo
enddo
enddo
PART I
“Integer” operations for
“ELEMmat”
In nonlinear cases, this part
should be done just once
(before initial iteration), as long
as mesh connectivity does not
change.
PART II
“Floating” operations for matrix
assembling/accumulation.
In nonlinear cases, this part is
repeated for every nonlinear
iteration.
08-APR22
127
Time for 3x643=786,432 DOF
ES
8.0 GFLOPS
Opteron
1.8GHz
DCRS
sec.
(MFLOPS)
DJDS original
sec.
(MFLOPS)
DJDS improved
sec.
(MFLOPS)
Matrix
28.6
(291)
34.2
(240)
12.5
(643)
Solver
360
(171)
21.7
(3246)
21.7
(3246)
Total
389
55.9
34.2
Matrix
10.2
(818)
12.4
(663)
21.2
(381)
Solver
225
(275)
271
(260)
271
(260)
Total
235
283
292
3.6 GFLOPS
08-APR22
128
Time for 3x643=786,432 DOF
ES
8.0 GFLOPS
Opteron
1.8GHz
DCRS
sec.
(MFLOPS)
DJDS original
sec.
(MFLOPS)
DJDS improved
sec.
(MFLOPS)
Matrix
28.6
(291)
34.2
(240)
12.5
(643)
Solver
360
(171)
21.7
(3246)
21.7
(3246)
Total
389
55.9
34.2
Matrix
10.2
(818)
12.4
(663)
21.2
(381)
Solver
225
(275)
271
271
Opteron上ではオリジナル版より
(260)
(260)
遅い：
3.6 GFLOPS
Total
235
・計算量増大
283・データのローカル性喪失，
292
08-APR22
129
Matrix
Solver
Matrix+Solver
60
60
sec.
40
DJDS on ES
original
50
40
sec.
50
30
30
20
20
10
10
0
0
41472
98304
192000
375000
DJDS on ES
improved
786432
41472
98304
DOF
375000
786432
DOF
300
300
DJDS on Opteron
original
200
DJDS on Opteron
improved
sec.
sec.
200
192000
100
100
0
0
41472
98304
192000
DOF
375000
786432
41472
98304
192000
375000
786432
08-APR22
130
Computation Time vs. Problem Size
300
40
Matrix
Solver
30
sec.
sec.
200
20
100
10
0
0
41472
98304
192000
375000
786432
41472
98304
192000
375000
786432
DOF
DOF
300
Total
ES (DJDS improved)
sec.
200
Opteron
(DJDS improved)
Opteron (DCRS)
100
0
41472
98304
192000
375000
786432
08-APR22
131
“Matrix” computation time for
improved version of DJDS
25
Integer
(Part I)
Floating
(Part II)
ES
20
15
sec.
sec.
20
25
15
10
10
5
5
0
Opteron
0
41472
98304
192000
DOF
375000
786432
41472
98304
192000
DOF
375000
786432
08-APR22
132
係数マトリクス生成部の
ベクトルプロセッサ向け最適化
• DJDS用のものは3倍ほど速くなったが，それでも
Opetron版DCRSより遅い
• 整数演算部（PART I）がボトルネック
– 計算の準備部分
– ベクトル長短くせざるを得ない，分岐が結構多い
• 実数演算部（PART II）は非常に速くなっている
– 整数演算部はメッシュが変わらない限り，最初に一回
計算するだけで良い
08-APR22
133
Suppose “virtual” mode where …
• On scalar processor
– “Integer” operation part
• On vector processor
– “floating” operation part
– linear solvers
• Scalar performance of ES (500MHz) is smaller
than that of Pentium III
08-APR22
134
Time for 3x643=786,432 DOF
ES
8.0 GFLOPS
Opteron
1.8GHz
DCRS
sec.
(MFLOPS)
DJDS improved
sec.
(MFLOPS)
DJDS virtual
sec.
(MFLOPS)
Matrix
28.6
(291)
12.5
(643)
1.88
(4431)
Solver
360
(171)
21.7
(3246)
21.7
(3246)
Total
389
34.2
23.6
Matrix
10.2
(818)
21.2
(381)
Solver
225
(275)
271
(260)
Total
235
292
3.6 GFLOPS
08-APR22
135
まとめ：有限要素法のベクトル化
• 決して簡単な作業ではない
• 有限要素法の特徴である局所性は，ベクトル計算
機の特性と相反する場合もある
– 全体マトリクスを対象とする「線形方程式ソルバー」
ではなんとかなるが
• 大幅なプログラム書換必要性のある場合もある
– 概して，必要記憶容量，行数，オペレーションは増加
• ベクトルプロセッサ向けに最適化すると，スカ
ラープロセッサで著しく性能が低下する場合もあ
る（マトリクス生成部の例）
• ベクトル化できるところを切り分けて重点的に
チューニングする（PART I/II）ことも重要
– 最後まで性能の悪いところは残る可能性あり