マルチコアPCクラスタ環境におけるBDD法のハイブリッド並列実装

MPI/Pthreadハイブリッド並列化による

領域分割型有限要素法の高速化

荻野正雄

九州大学 大学院工学研究院 機械工学部門

先駆的科学計算に関するフォーラム

2010

∼

GPGPU技術とその適応性を探る∼

2010年9月29日 九州大学伊都キャンパス

内容

• 研究背景・目的

• 領域分割法

(DDM)

–

DDMのハイブリッド並列実装

– 部分領域問題の解法

• 数値実験

– マルチコア

CPUにおけるCPU内コア間並列性能

– マルチコア

CPUクラスタにおける演算性能評価

• まとめ

– マルチコア

CPUとFEMについて

–

GPUとFEMについて

研究背景

•

CPUの

Multi-core

化

•

NUMA型設計などに基づく

Multiprocessing

•

CellやGPUなどのAccelerator活用, GPU統合型CPUなど

Heterogeneous computing

環境の普及

co re co re co re co re L3 cache Memory 4 core co re co re co re co re L3 cache Memory co re co re co re co re L3 cache Memory co re co re co re co re L3 cache Memory co re co re co re co re L3 cache Memory 4 x 4 core (NUMA) Memory coreL2 cache core core core core core core core 8 core Memory co re co re co re co re

Last level cache

GPU 4 + 1 core CPU + GPU CPU Memory GPU Memory

研究目的

•

Multi-core CPUの活用

– 将来的には

CPU当たり数十∼数百コア (Many core)

–

CPU内コア間の並列効率が重要となる

– ノード間など分散メモリ環境向けは

MPI並列

• 大規模問題解析向けのソルバ

―

– 反復法

– 反復回数

(全プロセッサの同期回数)を抑える前処理が重要

– 計算機アーキテクチャを意識した並列化アルゴリズム

領域分割法

(DDM)のMPI/POSIXスレッドによる

ハイブリッド並列実装を行い

, CPU内コア間並列

効率の改善を行う

反復法のマルチコア

CPUクラスタ向け並列化を考える

f

Ku =

Kp

q =

q

=

å

K

_i

p

例えば, CG法における 係数行列ベクトル積 - 解析領域を分割して並列プロセスに負荷分散 - 各プロセスが担当する部分行列とのベクトル積を計算

p

K

q

_i

=

_i 解くべき式(線形対称)

をマルチコア

CPU内で以下に高速に計算するか

線形代数ライブラリを OpenMPなどで並列化する

å

=

_K

_p

q

_{i j} さらに領域を分割し, MPI又はスレッド並列 こちらが一般的

領域分割法

(DDM)の導入

• 領域分割法

(Domain Decomposition Method)

– 解析領域を重なりがないように領域分割する

– 領域間境界上に縮約した方程式を反復法で解く

– 全自由度のまま解く場合に比べて反復法の収束性が良い

– さらに

, バランシング領域分割(BDD)法などの強力な前処理

を適用可能

DDMのマルチコアCPUクラスタ向け並列化を考える

DDM

f

Ku =

g

Su

_B

=

領域間境界

å

å

= =

=

=

=

N i i T i N i i i T i i i

R

u

,

K

R

K

R

,

f

R

f

u

1 1

1. N個に領域分割

2. 領域を領域内部Iと領域間境界Bに分割

ú û ù ê ë é = ú û ù ê ë é = ú û ù ê ë é = ú û ù ê ë é = i B i I i i B i I i i BB T i IB i IB i II i i B i I i _R R R , f f f , K K K K K , u u u 0 0

3. 領域間境界に縮約

(

)

å

å

-=

-=

=

i I i II T i IB i B T i B i IB i II T i IB i BB i i B i T i B

f

K

K

f

R

g

,

K

K

K

K

S

,

R

S

R

S

1 1

4. CG法などで解く

DDMにおける係数行列ベクトル積

þ

ý

ü

î

í

ì

-=

þ

ý

ü

î

í

ì

×

K

R

p

x

i B i i

0

i i II i

K

x

w

=

-1

þ

ý

ü

î

í

ì

=

þ

ý

ü

î

í

ì ×

p

R

w

K

q

_Bi i i i 基本境界条件処理 （行列ベクトル積） 求解処理 反力計算 （行列ベクトル積） for ( i = 0; i < N; i++)

å

å

= =

=

=

=

N i i T i B N i i Bi T i B

S

R

p

R

q

R

Sp

q

1 1

各領域で有限要素解析を行うことに相当する

p

R

S

q

_i

=

_{i Bi}

をマルチコア

CPU内で以下に高速に計算するか

DDMのマルチコアCPUクラスタ向け並列化を考える

g

Su

_B

=

å

=

=

_Sp

_R

_S

_R

_p

q

_BT_{i i Bi} 例えば, CG法における 係数行列ベクトル積 解くべき式(線形対称) - 解析領域を分割して並列プロセスに負荷分散 - 各プロセスが担当する部分行列とのベクトル積を計算 線形代数ライブラリを OpenMPなどで並列化する さらに領域を分割(総領域分割数を 増加)し, MPI又はスレッド並列 i i II i

K

x

w

=

-1 _{この求解処理が必要なことから}, 領域当たりの自由度数に強く依 存するため_{, 1領域は十分に小さ} くとった方が良い

経験的にこちらが最適

DDMのハイブリッド並列実装

þ

ý

ü

î

í

ì

-=

þ

ý

ü

î

í

ì

×

K

R

p

x

i B i i

0

i i II i

K

x

w

-1

=

þ

ý

ü

î

í

ì

=

þ

ý

ü

î

í

ì ×

p

R

w

K

y

_Bi i i i

å

=

=

N i i T i B

y

R

y

1

Sp

y =

部分: 1プロセッサの担当領域 部分領域： 計算の最小単位 基本境界条件処理 （行列ベクトル積） 求解処理 反力計算 （行列ベクトル積）

for ( i = 0; i < my_ndomain; i++)

係数行列ベクトル積

MPI/Pthreadハイブリッドプログラム例

void spartsolv(...) {

for (i = 0; i< n_domain; i++) { make_subdomain_matrices(...); domsolve(...); put_results(...); delete_subdomain_matrices(...); } }

void spartsolv_pthread(int threadid, int n_domain, ...) {

for (i = start_domid; i < end_domid; i++) { make_subdomain_matrices(...); domsolve(...); put_results(...); delete_subdomain_matrices(...); } }

void* spartsolv_callback(HDDM_MTArgs* mtarg) {

spartsolv_pthread(mtarg->threadid, mtarg->n_domain, ...); return(NULL);

}

void spartsolv(...) {

for (ithread = 0; ithread < NUM_THREADS; ithread++) { hddm_mt_args[ithread].threadid = ithread;

hddm_mt_args[ithread].n_domain = n_domain; ...

pthread_create(&hddm_pthr[ithread], NULL,

(void *)spartsolv_callback, &hddm_mt_args[ithread]); }

for (ithread = 0; ithread < NUM_THREADS; ithread++)

pthread_join(hddm_pthr[ithread], NULL);

for (ithread = 0; ithread < NUM_THREADS; ithread++) merge_results(...); } 元のFlat MPI向けプログラムにおける 部分領域ループをPthreadで並列化 - 元が関数呼び出し構造や引数が複雑な ので_{, POSIXスレッド並列を直接書いた方} が楽であった事例 - 引数の渡し方, メモリへの二重書き込み を行わないように注意が必要 - コンパイル時に -pthread をつける - プログラム中は #ifdef _REENTRANT ∼ #endifでFlat MPI用コードと共存

部分領域問題

w

_i

=K

_IIi

-1

x

_i

の求解方法

(1/2)

• 線形代数ソルバ

―

– 直接法

(Direct solver)

修正コレスキー分解, 前進/後退代入, スカイライン型記憶

– 反復法

(Iterative solver)

CG法, 前処理, 非零型記憶

• 主記憶装置への保存

– データ保存型

(Storage type)

メッシュ, 係数行列, LDLt分解行列, 前処理行列等

– データ非保存型

(Storage-Free type)

メッシュのみ

部分領域問題

w

_i

=K

_IIi

-1

x

_i

の求解方法

(2/2)

Kawai's classification

• 直接法ベース

(Direct-solver-based)

– データ保存型

(Storage type)

DS

• 係数行列とその分解結果をメモリーに保存 • 各CG反復では, 前進/後退代入を行う

– データ非保存型

(Storage-Free type)

DSF

• 各CG反復では, 係数行列の作成, LDLt分解, 前進/後退代入を行う

• 反復法ベース

(Iterative-solver-based)

– データ保存型

(Storage type)

IS

• 係数行列と前処理行列をメモリーに保存 • 各CG反復では, PCG法で解く

– データ非保存型

(Storage-Free type)

ISF

• 各CG反復では, 係数行列の作成, 前処理行列の作成, PCG法で解く

* Kawai et al., USNCCM 2009

メモリー使用量削減と

共に, メモリーアクセス

数値実験

1

•

100万自由度規模, 線形弾性問題

• 使用計算機

– Core2Duo: Intel Core2Duo E6600 (L2 4MB), DDR2 SDRAM 8GB – Core2Quad: Intel Core2Quad Q6600 (L2 8MB), DDR2 SDRAM 8GB – Core i7: Intel Core i7 920 (L3 8MB), DDR3 SDRAM 12GB

– Xeon: Intel Xeon E5345 (L2 8MB), FB-DIMM 32GB

• 実験条件

– パラメータ

: 領域自由度数, 使用コア数(1 or 全コア)

–

DDM実装:

DS

,

DSF

,

ISF

•

1領域計算時間のシングル/マルチコア比

•

DDM反復1回に要する計算時間

– 並列化

:

Flat MPI

領域自由度 vs 1部分領域の計算時間比(1コア使用/マルチコア使用)

Core2Duo Core2Quad

領域自由度 vs DDM反復1回の計算時間

Core2Duo Core2Quad

数値実験

2

• 線形弾性問題

–

1,200万自由度

–

5,000部分領域 (2,400自由度/部分領域)

•

DDM実装

– 部分領域問題解析

: DS, DSF, ISF

– 並列化

: Flat MPI, MPI/Pthread Hybrid

– 前処理

: BDD

• 計算機環境

–

Corei7クラスタ

Intel Corei7 920 (2.66GHz QC, L2 256KB/core, L3 8MB/CPU), Intel Compiler 11.1 and MPICH2-1.2

–

HA8000@東大ITC (T2K)

AMD Opteron 8356 (2.3GHz QC, L2 256KB/core, L3 2MB/CPU) x 4 / Node, Intel Compiler 11.0 and MPICH-1.2.7

実験方法

Memory core core core core Memory core core core core Memory core core core core Memory core core core core

①

①

①

①

②

②

②

②

③

④ ③

④

③

④ ③

④

CPU内コア間並列効率測定

-

CPU内の使用コア数を１つず

つ増やした場合の計算時間を

比較

- Flat MPI実装ではMPIプロセス

数

, Hybrid実装ではスレッド数を

増やす

HA8000ではNUMA最適化として"--cpunodebind"でCPUを, "--membind"でメモリを指定

実験結果

- Intel Corei7クラスタ

PC8台を使用, CPU内コアの使用数を1,2,3,4と変えた場合の結果 細線: Flat MPI, 太線: MPI-Pthread, DS, DSF, ISF

コア数と計算時間

コア数とスピードアップ率

Matrix Storage-Free型は85%以上のCPU内コア間並列効率

実際の計算時間でも

Matrix Storage型に近づく(まだ約5倍)

< 40% = 100% > 85%

実験結果

- HA8000@東大ITC (T2K)

コア数と計算時間

コア数とスピードアップ率

16ノードを使用, CPU内コアの使用数を1,2,3,4と変えた場合の結果 細線: Flat MPI, 太線: MPI-Pthread, DS, DSF, ISF

Corei7に比べて並列効率が低下しているのはBDD前処理の非並列化

部分による影響

, しかしNUMAアーキテクチャにおいても高い並列効率

が得られている

マルチコア

CPUにおけるまとめ

•

Multi-core CPU環境では, 演算量に比べてメモリアク

セス量が少ないアルゴリズムであれば

, 高いCPU内

コア間並列効率が得られる

.

• 再利用可能な演算結果をメモリにストアするよりも

,

キャッシュ上で再構築する方が速くなるかもしれな

い

.

•

CPU, メモリーの階層構造を意識した数値計算アル

ゴリズムの導入を検討することも重要である

.

• 将来の

Many Core環境(16コア以上)においては, ISF

実装が

DS実装よりも高速になる見通しである.

GPUとFEMの現状

• 工学分野

(構造解析)

– 非構造格子

+有限要素法

→ 疎行列, インデックス参照アクセス

•

GPUとFEM

– 性能を出しづらい

(数%, 良くて数十%の速度改善)

→ コストパフォーマンスが悪い(移植, 価格)

GPUや今後のCPUとの付き合い方

•

Heterogeneous computing, Memory hierarchy

を意識して数値計算アルゴリズムを見直す

→ CPUとGPUの役割分担, メモリ利用方法, 線形代数ライブラリ

への依存度を自由に調整できるように解法レベルで考える

• 例えば

DDMのCPU+GPU向け実装

– 部分領域演算は

GPU, 前処理はCPU

– 三角分解は用いずに

, 行列ベクトル積だけで全体を構成

• 線形代数ライブラリについて

– 密行列

, バンド行列, できれば疎行列

– 行列ベクトル積

, できれば三角分解

–

GPU向けに最適化

–

CPU向けにはCPU内コア全てや分散メモリ対応よりも適当なコ

ア数で最適化

–

DDMをやるので, 数百元程度までの最適化でも良い

896MB VRAM version2.1 CUDA NVIDIA GeForce GTX260 GPU gcc version4.2.4 compiler 6GB memory Intel Core i7 920 (2.67GHz) CPU Fedora 10 OS

実験環境

Comparison of data transfer speed

111.9GB/s 25.5GB/s

Data transfer

GPU CPU

Matrix-vector product benchmark

About 25 G flops (DP) about 30 % of peak

(80Gflops)

x6 faster than CPU

DDM hot spot benchmark on GPU

Kawai et al., WCCM 2010

DDMの主要演算を密行列-ベクトル積に落とし込んだと想定し,