並列有限要素法への道
ー並列データ構造ー
C言語編
中島 研吾
東京大学情報基盤センター

Intro pFEM

2並列計算の目的• 高速，大規模– 「大規模」の方が「新しい科学」という観点からのウェイトとしては高い。しかし，「高速」ももちろん重要である。
• 細かいメッシュ• +複雑• 理想：Scalable– N倍の規模の計算をN倍のCPUを使って，「同じ時間で」解く（大規模性の追求：Weak Scaling）• 実際はそうは行かない– 例：共役勾配法⇒問題規模が大きくなると反復回数が増加– 同じ問題をN倍のCPUを使って「1/Nの時間で」解く・・・という場
合もある（高速性の追求：Strong Scaling）• これも余り簡単な話では無い

Intro pFEM

3並列計算とは ?（1/2）• より大規模で複雑な問題を高速に解きたい

Homogeneous/Heterogeneous
Porous Media
Lawrence Livermore National Laboratory

このように不均質な場を模擬するには非常に細かいメッシュが必要となる
Homogeneous

HeterogeneousIntro pFEM

4並列計算とは ?（2/2）• 1GB程度のPC → <106メッシュが限界：FEM– 1000km×1000km×100kmの領域（西南日本）を1kmメッシュで切ると108メッシュになる• 大規模データ → 領域分割，局所データ並列処理• 全体系計算 → 領域間の通信が必要

領域分割
大規模データ

局所
局所
データ データ
局所
局所
データ データ

局所データ

局所データ
通信

局所データ

局所データ

Intro pFEM

5通信とは ?
• 並列計算とはデータと処理をできるだけ「局所的（local）」に実施すること。
– 要素マトリクスの計算– 有限要素法の計算は本来並列計算向けである• 「大域的（global）」な情報を必要とする場合に通信が生じる（必要となる）。
– 全体マトリクスを線形ソルバーで解く

Intro pFEM

6

並列有限要素法の処理：SPMD
巨大な解析対象 → 局所分散データ，領域分割（Partitioning）
有限要素コードは領域ごとに係数マトリクスを生成：要素単位の処理によって可能：シリアルコードと変わらない処理
グローバル処理，通信は線形ソルバーのみで生じる
内積，行列ベクトル積，前処理
局所データ

FEMコード

線形ソルバー

MPI
局所データ

FEMコード

線形ソルバー

MPI
局所データ

FEMコード

線形ソルバー

MPI
局所データ

FEMコード

線形ソルバーIntro pFEM

7

並列有限要素法プログラムの開発• 前頁のようなオペレーションを実現するためのデータ構
造が重要– アプリケーションの「並列化」にあたって重要なのは，適切な局所分散データ構造の設計である。

• 前処理付反復法• マトリクス生成：ローカルに処理

四角形要素
5671212583364711516

Intro pFEM

376228

これではマトリクス生成に必要な情報は不十分

321「節点ベース（領域ごとの節点数がバランスする）」の分割
自由度：節点上で定義

47232833

マトリクス生成のためには，オーバーラップ部分の要素と節点の情報が必要

48

並列有限要素法の局所データ構造
節点ベース：Node-based partitioning 局所データに含まれるもの :






その領域に本来含まれる節点
それらの節点を含む要素
本来領域外であるが，それらの要素に含まれる節点節点は以下の3種類に分類


内点：Internal nodes
外点：External nodes
境界点：Boundary nodes

その領域に本来含まれる節点
本来領域外であるがマトリクス生成に必要な節点
他の領域の「外点」となっている節点



領域間の通信テーブル



領域間の接続をのぞくと，大域的な情報は不要
Intro pFEM

有限要素法の特性：要素で閉じた計算
9Node-based Partitioning
internal nodes - elements - external nodes
PE#3
PE#12122

4

PE#0235171812117619132２PE#3Intro pFEM

３４PE#211310789111012510141345101238911121091112963１7PE#0159615201481225241166848465512PE#37712PE#2310

Node-based Partitioning
internal nodes - elements - external nodes
Partitioned nodes themselves (Internal Nodes) 内点
Elements which include Internal Nodes 内点を含む要素
External Nodes included in the Elements 外点
in overlapped region among partitions.

Info of External Nodes are required for completely localelement–based operations on each processor.
151413108
Intro pFEM

67451239 111211

Node-basedPartitioning
マトリクス生成時の通信は不要
internal nodes - elements - external nodes

Partitioned nodes themselves (Internal Nodes) 内点
Elements which include Internal Nodes 内点を含む要素
External Nodes included in the Elements 外点
in overlapped region among partitions.

Info of External Nodes are required for completely local
element–based operations on each processor.
151413108
Intro pFEM

67451239 111212

Parallel Computing in FEM
SPMD: Single-Program Multiple-Data
Local Data

FEM codeLinear Solvers
MPILocal Data

FEM codeLinear Solvers
MPILocal Data

FEM codeLinear Solvers
MPILocal Data
Intro pFEM

FEM codeLinear Solvers
13

Parallel Computing in FEM
SPMD: Single-Program Multiple-Data
Local Data

FEM codeLinear Solvers
MPILocal Data

FEM codeLinear Solvers
MPILocal Data

FEM codeLinear Solvers
MPILocal Data
Intro pFEM

FEM codeLinear Solvers
14

Parallel Computing in FEM
SPMD: Single-Program Multiple-Data
15614134510123911121091112465712345612127891011101252Local Data

FEM codeLinear Solvers

Local Data

FEM codeLinear Solvers
MPI841

Linear Solvers

963

FEM code

MPI113Local Data

MPI88

Intro pFEM

77Local Data

FEM codeLinear Solvers
15

Parallel Computing in FEM
SPMD: Single-Program Multiple-Data
156714134510

123Local Data

818

911121091112

123456128

91011101252

FEM codeLinear Solvers

961

Linear Solvers

MPILocalData
384
Intro pFEM

FEM code

MPI113Local2 Data
7

Linear Solvers
MPILocal4 Data657

FEM code

7

FEM codeLinear Solvers
16

Parallel Computing in FEM
SPMD: Single-Program Multiple-Data
Local Data

FEM codeLinear Solvers
MPILocal Data

FEM codeLinear Solvers
MPI

Local Data

FEM codeLinear Solvers
MPI

Local Data
Intro pFEM

FEM codeLinear Solvers
17通信とは何か ?
• 「外点」の情報を外部の領域からもらってくること• 「通信テーブル」にその情報が含まれている

Intro pFEM

18

Intro pFEM

19

一次元問題：11要素，12節点，3領域
0000111122223333444455556666777788889999101010101111Intro pFEM

20

一次元問題：11要素，12節点，3領域
001122334455667897891010110011223344556677889910101111

Intro pFEM

21

節点がバランスするよう分割：内点
0011223344556678978910101100112#0

2334455667#1

7889910101111#2

Intro pFEM

22

内点だけで分割するとマトリクス不完全
0123456789101101#02345#16789#21011

Intro pFEM

23

要素＋外点
00112233401#02333445556778899101011466778#1789#21011

Intro pFEM

24

一次元問題：11要素，12節点，3領域
0001112334344556778899101011122306778#02334455667#17889910101111#2

Intro pFEM

25

一次元問題：11要素，12節点，3領域
000000111111222222333333334444555566667777888899991010101011114445566777788899101011Intro pFEM

26

SPMD向け局所番号付け
内点が1~N（0~N）となっていれば，もとのプログラムと同じ
外点の番号は?
00

112233334445566770

123?

?

01277883?

?

0891910210113

MPIによる並列化：SPMD• Single Program/Instruction Multiple Data
• 基本的に各プロセスは「同じことをやる」が「データが違う」– 大規模なデータを分割し，各部分について各プロセス（プロセッサ）が計算する• 全体データと局所データ，全体番号と局所番号• 通信以外は単体CPUと同じ，というのが理想

PE #0

PE #1

PE #2

PE #M-1

Program

Program

Program

ProgramData #0

Data #1

Data #2

Data #M-1

27

Intro pFEM

28

SPMD向け局所番号付け
内点が1~N（0~N）となっていれば，もとのプログラムと同じ
外点の番号は?，N+1,N+2（N,N+1）0123440123

540

123

Intro-pFEM

29

並列有限要素法：マトリクス生成
各要素での積分，要素マトリクス⇒（局所）全体マトリクス生成の計算は，内点・外点，それらを含む要素の情報があれば可能#000112#1234334001122#2344350011223

Intro pFEM

30

有限要素法の処理：プログラム• 初期化：並列計算可– 制御変数読み込み– 座標読み込み⇒要素生成（N:節点数，NE：要素数）– 配列初期化（全体マトリクス，要素マトリクス）– 要素⇒全体マトリクスマッピング（Index，Item）• マトリクス生成：並列計算可– 要素単位の処理（do icel= 1, NE）• 要素マトリクス計算• 全体マトリクスへの重ね合わせ– 境界条件の処理• 連立一次方程式：?
– 共役勾配法（CG）

Intro pFEM

31

前処理付き共役勾配法
Preconditioned Conjugate Gradient Method （CG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|
end

 前処理
 対角スケーリング

 並列処理が必要なプロセス
 内積
 行列ベクトル積0D1 0 ... 00 D002M    ...
...
... 

00D0N 1
 0 0 ... 0DN Intro pFEM

32

前処理，ベクトル定数倍の加減
局所的な計算（内点のみ）が可能⇒並列処理
0/*//-- {z}= [Minv]{r}*/for(i=0;i<N;i++){W[Z][i] = W[DD][i] * W[R][i];}1234/*//-- {x}= {x} + ALPHA*{p}// {r}= {r} - ALPHA*{q}*/for(i=0;i<N;i++){U[i]+= Alpha * W[P][i];W[R][i] -= Alpha * W[Q][i];}567891011

Intro pFEM

33

内積全体で和をとる必要がある⇒通信?
0/*//-- ALPHA= RHO / {p}{q}*/C1 = 0.0;
for(i=0;i<N;i++){C1 += W[P][i] * W[Q][i];}Alpha = Rho / C1;

1234567891011

Intro pFEM

34
P#0 A0 B0 C0 D0

MPI_ReduceP#1 A1 B1 C1 D1P#0 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3Reduce

P#1P#2 A2 B2 C2 D2P#2P#3 A3 B3 C3 D3P#3• Reduces values on all processes to a single value– Summation, Product, Max, Min etc.

• MPI_Reduce (sendbuf,recvbuf,count,datatype,op,root,comm)– sendbuf– recvbuf

choicechoiceIO– count
int
I– datatype MPI_Datatype I
FORTRANC– op

starting address of send buffer
starting address receive buffer
type is defined by ”datatype”
number of elements in send/receive buffer
data type of elements of send/recive buffer

MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION, MPI_CHARACTER etc.
MPI_INT, MPI_FLOAT, MPI_DOUBLE, MPI_CHAR etc
MPI_OpI

reduce operation

MPI_MAX, MPI_MIN, MPI_SUM, MPI_PROD, MPI_LAND, MPI_BAND etc

Users can define operations by MPI_OP_CREATE– root– commintI
MPI_Comm I

rank of root process
communicatorCIntro pFEM

35
P#0 A0 B0 C0 D0

MPI_Bcast

P#1

P#0 A0 B0 C0 D0Broadcast

P#1 A0 B0 C0 D0

P#2

P#2 A0 B0 C0 D0

P#3

P#3 A0 B0 C0 D0• Broadcasts a message from the process with rank "root" to all otherprocesses of the communicator• MPI_Bcast (buffer,count,datatype,root,comm)– bufferchoiceI/O– count
int
I– datatype MPI_Datatype I
FORTRANC– root
– commstarting address of buffer
type is defined by ”datatype”
number of elements in send/receive buffer
data type of elements of send/recive buffer

MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION, MPI_CHARACTER etc.
MPI_INT, MPI_FLOAT, MPI_DOUBLE, MPI_CHAR etc.

int

I
MPI_Comm I

rank of root process
communicatorCIntro pFEM

36

MPI_Allreduce

P#0 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3

P#0 A0 B0 C0 D0P#1 A1 B1 C1 D1All reduceP#1 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3

P#2 A2 B2 C2 D2P#2 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3

P#3 A3 B3 C3 D3P#3 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3• MPI_Reduce + MPI_Bcast• Summation (of dot products) and MAX/MIN values are likely to utilized ineach process• call MPI_Allreduce(sendbuf,recvbuf,count,datatype,op, comm)– sendbuf– recvbuf

choicechoiceIO– countint
I– datatype MPI_Datatype I– op– commMPI_OpI
MPI_Comm I

starting address of send buffer
starting address receive buffer
type is defined by ”datatype”
number of elements in send/receive buffer
data type of elements of send/recive buffer
reduce operation
communicatorCIntro pFEM

37

“op” of MPI_Reduce/Allreduce
MPI_Reduce(sendbuf,recvbuf,count,datatype,op,root,comm)• MPI_MAX，MPI_MIN• MPI_SUM，MPI_PROD• MPI_LAND

Max, Min
Summation, ProductLogical AND

CIntro pFEM

38

前処理付き共役勾配法
Preconditioned Conjugate Gradient Method （CG）
Compute r(0)= b-[A]x(0)for i= 1, 2, …
solve [M]z(i-1)= r(i-1)
i-1= r(i-1) z(i-1)
if i=1p(1)= z(0)elsei-1= i-1/i-2p(i)= z(i-1) + i-1 p(i-1)
endifq(i)= [A]p(i)i = i-1/p(i)q(i)x(i)= x(i-1) + ip(i)r(i)= r(i-1) - iq(i)check convergence |r|end

 前処理
 対角スケーリング

 並列処理が必要なプロセス
 内積
 行列ベクトル積0
D1 0 ... 00 D002M    ...
...
... 
00D0N 1
 0 0 ... 0DN Intro pFEM

39

行列ベクトル積
外点の値が必要⇒通信?
/*//-- {q}= [A]{p}*/for(i=0;i<N;i++){W[Q][i] = Diag[i] * W[P][i];for(j=Index[i];j<Index[i+1];j++){W[Q][i] += AMat[j]*W[P][Item[j]];}}401235

Intro pFEM40行列ベクトル積：ローカルに計算実施可能
012345600112233445567891011=677889910101111Intro pFEM41行列ベクトル積：ローカルに計算実施可能
012345600112233445567891011=677889910101111Intro pFEM

42行列ベクトル積：ローカルに計算実施可能
0123012001122330011230123=23300112233Intro pFEM43行列ベクトル積：ローカル計算 #00120011=23334012324Intro pFEM

44行列ベクトル積：ローカル計算 #1012120112300233001123=3=2345401235

Intro pFEM

45行列ベクトル積：ローカル計算 #2013

120112200233001123=3=23440123

Intro pFEM

461D FEM: 12 nodes/11 elem’s/3 domains
000000111111222222333333334444555566667777888899991010101011114445566777788899101011

Intro pFEM

47

1D FEM: 12 nodes/11 elem’s/3 domains
Local ID: Starting from 0 for node and elem at each domain#0001

12#123433400

1

122#23443

500

1

1223

Intro pFEM

48

1D FEM: 12 nodes/11 elem’s/3 domains
Internal/External Nodes

#0001

12

#1234334001122#2344350011223

Intro pFEM49

1対1通信とは ?
• グループ通信：Collective Communication– MPI_Reduce, MPI_Scatter/Gather など– 同じコミュニケータ内の全プロセスと通信する– 適用分野• 境界要素法，スペクトル法，分子動力学等グローバルな相互作用のある手法• 内積，最大値などのオペレーション• 1対1通信：Point-to-Point#0– MPI_Send, MPI_Receive– 特定のプロセスとのみ通信がある

0• 隣接領域– 適用分野• 差分法，有限要素法などローカルな情報を使う手法

01

12

#1

234

3

3

400

1

122

#2

34

4

3

500

1

1223

Intro pFEM

50

SEND（送信）：境界点の送信
送信バッファの連続したデータを隣接プロセスに送る• MPI_Isend(sendbuf,count,datatype,dest,tag,comm,request)––––sendbuf 任意
count整数
datatype 整数
dest整数

IIII

送信バッファの先頭アドレス，メッセージのサイズ
メッセージのデータタイプ
宛先プロセスのアドレス（ランク）PE#3
456121567

PE#0
12789111012584

12

PE#1

1014134510

1238911121091112963113784

7

1652

PE#2

350

Intro pFEM

51

C

MPI_Isend•

送信バッファ「sendbuf」内の，連続した「count」個の送信メッセージを，タグ「tag」を付けて，コミュニケータ内の，「dest」に送信する。「MPI_Waitall」を呼ぶまで，送信バッファの内容を更新してはならない。

• MPI_Isend(sendbuf,count,datatype,dest,tag,comm,request)–––––

sendbuf 任意
count整数
datatype 整数
dest整数tag整数– comm– request

MPI_CommIIIIIIMPI_Request O

送信バッファの先頭アドレス，メッセージのサイズ
メッセージのデータタイプ
宛先プロセスのアドレス（ランク）メッセージタグ，送信メッセージの種類を区別するときに使用。
通常は「0」でよい。同じメッセージタグ番号同士で通信。
コミュニケータを指定する
通信識別子。MPI_Waitallで使用。
（配列：サイズは同期する必要のある「MPI_Isend」呼び出し数（通常は隣接プロセス数など））:C言語については後述

Intro pFEM

52

RECV（受信）：外点への受信受信バッファに隣接プロセスから連続したデータを受け取る• MPI_Irecv(recvbuf,count,datatype,dest,tag,comm,request)––––

recvbuf 任意
count整数
datatype 整数
dest整数

IIII

受信バッファの先頭アドレス，メッセージのサイズ
メッセージのデータタイプ
宛先プロセスのアドレス（ランク）PE#3456121567

PE#01278911101258412

PE#11014134510123891112109111296311378471652

PE#2352

Intro pFEM53

C

MPI_Irecv•

受信バッファ「recvbuf」内の，連続した「count」個の送信メッセージを，タグ「tag」を付けて，コミュニケータ内の，「dest」から受信する。「MPI_Waitall」を呼ぶまで，受信バッファの内容を利用した処理を実施してはならない。

• MPI_Irecv(recvbuf,count,datatype,dest,tag,comm,request)–––––

recvbuf 任意
count整数
datatype 整数
dest整数tag整数– comm– request

MPI_CommIIII
I

I
MPI_Request O

受信バッファの先頭アドレス，メッセージのサイズ
メッセージのデータタイプ
宛先プロセスのアドレス（ランク）メッセージタグ，受信メッセージの種類を区別するときに使用。
通常は「0」でよい。同じメッセージタグ番号同士で通信。
コミュニケータを指定する
通信識別子。MPI_Waitallで使用。
（配列：サイズは同期する必要のある「MPI_Irecv」呼び出し数（通常は隣接プロセス数など））：C言語については後述

Intro pFEM

54

C

MPI_Waitall• 1対1非ブロッキング通信関数である「MPI_Isend」と「MPI_Irecv」を使用した場合，プロセスの同期を取るのに使用する。
• 送信時はこの「MPI_Waitall」を呼ぶ前に送信バッファの内容を変更してはならない。
受信時は「MPI_Waitall」を呼ぶ前に受信バッファの内容を利用してはならない。
• 整合性が取れていれば， 「MPI_Isend」と「MPI_Irecv」を同時に同期してもよい。
– 「MPI_Isend/Irecv」で同じ通信識別子を使用すること•「MPI_Barrier」と同じような機能であるが，代用はできない。
– 実装にもよるが，「request」，「status」の内容が正しく更新されず，何度も「MPI_Isend/Irecv」を呼び出すと処理が遅くなる，というような経験もある。

• MPI_Waitall(count,request,status)– count– request

整数
整数

I
I/O

– status

MPI_Status

O

同期する必要のある「MPI_ISEND」 ，「MPI_RECV」呼び出し数。
通信識別子。「MPI_ISEND」，「MPI_IRECV」で利用した識別子名に対応。（配列サイズ：（count））
状況オブジェクト配列
MPI_STATUS_SIZE： “mpif.h”,”mpi.h”で定められるパラメータ：C言語については後述

Intro pFEM

55

Node-based Partitioning
internal nodes - elements - external nodes
PE#1
PE#1
21

224PE#0235171812117619132２PE#3３４PE#211310789111012510141345101238911121091112963１7PE#0159615201481225241166848465512

PE#37712

PE#23

Intro pFEM

56

各領域データ（局所データ）仕様
1091112• 内点，外点（internal/external nodes）– 内点～外点となるように局所番号をつける

84716 •523

隣接領域情報– オーバーラップ要素を共有する領域– 隣接領域数，番号• 外点情報– どの領域から，何個の，どの外点の情報を「受信：import」するか• 境界点情報– 何個の，どの境界点の情報を，どの領域に「送信：export」するか

Intro pFEM

57

External Nodes（外点）: RECEIVE
PE#2 : receive information for “external nodes”
1567PE#0111051284PE#31213451012389111210911129631478471652PE#23

Intro pFEM

58

Boundary Nodes（境界点）: SEND
PE#2 : send information on “boundary nodes”
1567

PE#0111051284PE#31345101238911121091112963148465

PE#21277123

Intro pFEM

59

並列計算向け局所（分散）データ構造• 差分法，有限要素法，有限体積法等係数が疎行列のアプリ
ケーションについては領域間通信はこのような局所（分散）データによって実施可能– SPMD– 内点～外点の順に「局所」番号付け– 通信テーブル：一般化された通信テーブル• 適切なデータ構造が定められれば，処理は非常に簡単。
– 送信バッファに「境界点」の値を代入– 送信，受信– 受信バッファの値を「外点」の値として更新