並列有限要素法による一次元定常熱伝導解析プログラム
Fortran編
中島 研吾
東京大学情報基盤センター
お試しアカウント付き講習会「MPI応用編：並列有限要素法」

S2-ref• 問題の概要，実行方法• プログラムの説明
• 計算例

2

FEM1D

3

対象とする問題：一次元熱伝導問題
体積当たり一様発熱Qɺ∂  ∂T  ɺ
λ+Q = 0∂x  ∂x 
x= xmaxx=0 (xmin)• 一様な：断面積A，熱伝導率λ• 体積当たり一様発熱（時間当たり）〔QL-3T-1〕• 境界条件– x=0 ：T= 0– x=xmax：

（固定）

∂T
=（断熱）
0
∂x

Qɺ

FEM1D

4

対象とする問題：一次元熱伝導問題
体積当たり一様発熱Qɺ∂  ∂T  ɺ
λ+Q = 0∂x  ∂x 
x= xmaxx=0 (xmin)• 一様な：断面積A，熱伝導率λ• 体積当たり一様発熱（時間当たり）〔QL-3T-1〕• 境界条件– x=0 ：T= 0– x=xmax：

（固定）

∂T
=（断熱）
0
∂x

Qɺ

FEM1D5

解析解
体積当たり一様発熱Qɺ∂  ∂T  ɺ
λ+Q = 0∂x  ∂x 

x=0 (xmin)

x= xmax

T = 0@ x = 0∂T= 0 @ x = x max∂x

λ T ′′ = − Qɺ
λ T ′ = − Qɺ x + C1  C1 = Qɺ x max , T ′ = 0 @ x = x max
1 ɺ 2
λ T = − Q x + C1 x + C 2  C 2 = 0 , T = 0 @ x = 02
1 ɺ 2 Qɺ x max∴ T =−Qx +x2λλS2-ref6

ファイルコピー，コンパイル
Fortran>$ cd /work/gt00/t00XXX/pFEM>$ cp /work/gt00/z30088/pFEM/F/s2r-f.tar .
>$ tar xvf s2r-f.tar
C>$ cd /work/gt00/t00XXX/pFEM>$ cp /work/gt00/z30088/pFEM/C/s2r-c.tar .
>$ tar xvf s2r-c.tar
Confirm/Compile>$ cd mpi/S2-ref>$ mpiifort -align array64byte -O3 -axCORE-AVX512 1d.f -o 1da>$ mpiifort -O3 1d.f -o 1db>$ mpiicc –align -O3 -axCORE-AVX512 1d.c –o 1da>$ mpiicc -O3 1d.c –o 1db<$O-S2r> = <$O-TOP>/mpi/S2-refS2-ref

7

制御ファイル：input.dat
制御ファイル4
1.0 1.0100
1.e-8

input.dat

1.0

NE（要素数）∆x（要素長さL），Q, A, λ反復回数（CG法後述）
CG法の反復打切誤差

1.0

∆x=1
1
x=0

2
x=1

3
x=2

5

4

3

2

1

4
x=3

x=4

要素番号
節点番号（全体）

S2-ref

g16.sh: 8-nodes, 256-cores, 16x2x8#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=8#PJM --mpi proc=256#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test.lst

256/8= 32 cores/node

mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db
export I_MPI_PIN_PROCESSOR_LIST=0-15,28-43
mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

32-cores are randomly
selected from 56-coreson the node
32-cores on eachsocket are assigned.
A little bit more stable8

S2-ref

g24.sh: 8-nodes, 384-cores, 24x2x8#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=8#PJM --mpi proc=384#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test.lst

384/8= 48 cores/node

mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db
export I_MPI_PIN_PROCESSOR_LIST=0-23,28-51
mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

48-cores are randomly
selected from 56-coreson the node
24-cores on eachsocket are assigned.
A little bit more stable9S2-ref

g28.sh: 8-nodes, 448-cores, 28x2x8#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=8#PJM --mpi proc=448#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test.lst

448/8= 56 cores/node

mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db
export I_MPI_PIN_PROCESSOR_LIST=0-55
mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

10

11

NUMA Architecture
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPSIntel® Xeon®
Platinum 8280

UPI

(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS

UPI
UPI

Soc. #0: 0th-27th cores

2933 MHz×6ch
140.8 GB/sec

Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores• Oakbridge-CX (OBCX)– 2 Sockets (CPU’s) of Intel CLX– 各ソケットは28コア，合計56コア• NUMAアーキテクチャ（Non-Uniform Memory Access）– メモリは各CPUに搭載されていて独立，異なるCPUのローカルメ
モリ上のデータをアクセスすることは可能– ローカルメモリ上のデータを使って計算するのが効率的• numactl –l : ローカルメモリ使用，本オプションによる遅くなる場合有り

S2-ref

12「並列計算」の手順• 制御ファイル，「全要素数」を読み込む• 内部で「局所分散メッシュデータ」を生成する• マトリクス生成• 共役勾配法によりマトリクスを解く• 元のプログラムとほとんど変わらない

S2-ref• 問題の概要，実行方法• プログラムの説明
• 計算例

13

S2-ref

14

並列有限要素法の局所データ構造• 節点ベース：Node-based partitioning• 局所データに含まれるもの :– その領域に本来含まれる節点– それらの節点を含む要素– 本来領域外であるが，それらの要素に含まれる節点• 節点は以下の3種類に分類– 内点：Internal nodes– 外点：External nodes– 境界点：Boundary nodes

その領域に本来含まれる節点
本来領域外であるがマトリクス生成に必要な節点
他の領域の「外点」となっている節点• 領域間の通信テーブル• 領域間の接続をのぞくと，大域的な情報は不要– 有限要素法の特性：要素で閉じた計算
14

S2-ref

15

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

22

4

PE#0
23

5

17

18

12

11

7

6

19

13

15

9

２

PE#3

３

４

PE#2

6

7

PE#0
2

11

3

10

7

8

9

11

10

12

5

10

14

13

4

5

10

1

2

3

8

9

11

12

10

9

11

12

9

6

3
１

15

20

14

8

12

25

24

1
16

6

8

4

8

4

6

5

5
1

2

PE#1

7

7

1

2

PE#2

3
15

S2-ref

16

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.

外点の情報は各プロセスで独立に要素マトリクス生成を実施するのに必要⇒マトリクス生成時の通信は不要となる
15

8

6

7

14

13

4

5

10

1

2

3

9

11

12
16

S2-ref17

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

2
1

1

3
2

2
1

1

3

4

3
2

5
4

3

2

2
1

4
3

5
4

4
3

6
5

7
6

6
5

8
7

7
6

9
8

8
7

10
9

9
8

11
10

10
9

12
11

11
10

12
11

5
4

4

5
4

6
5

7
6

8
7

9
8

8

9
8

10
9

11
10

12
11

Intro pFEM

18

各領域データ（局所データ）仕様
10

9

8

11

4

12

6

5

• 内点，外点（internal/external
nodes）– 内点～外点となるように局所番号をつける• 隣接領域情報
7

1

2

3

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

S2-ref

19

SPMD向け局所番号付け
局所節点・要素番号は「1」からスタート

#0

1

2
1

3

4
3

2

#1

5
4

5

1
4

2
1

3
2

4
3

#2

6
5

5

1
4

2
1

3
2

4
3

Intro pFEM

20

一次元問題：11要素，12節点，3領域
内点・外点・境界点
境界点は，他の領域の外点

#0

1

2
1

3
3

2

#1

X

4

5
4

5

1
4

2
1

3
2

4
3

外点：External Nodes
#2

Y

6
5

境界点：Boundary Nodes

5

1
4

2
1

3
2

4
3

S2-ref

21

1D FEM: Numbering of Local ID
1

2
1

N+1

1
N

N+1

2
1

1
N

…

…
2

2
1

N
N-1

2

N

#PETot-1:
N+1 nodes
N elements

N

Others (General):N+2 nodes
N+1 elements

N-1

…
2

N+1
N

N-1

#0:
N+1 nodes
N elements

N+2
N+1

S2-ref22

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

1

2
1

3

4
3

2

#1

5
4

5

1
4

2
1

3
2

4
3

#2

6
5

5

1
4

2
1

3
2

4
3

Intro pFEM

23

前処理付き共役勾配法
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 D

0
0
2


[M ] =  ...
...
... 


D
0
0
0
N −1


 0 0 ... 0
DN 

Intro pFEM

24

前処理，ベクトル定数倍の加減
局所的な計算（内点のみ）が可能⇒並列処理
1!C!C-- {z}= [Minv]{r}do i= 1, N
W(i,Z)= W(i,DD) * W(i,R)
enddo

!C
!C-- {x}= {x} + ALPHA*{p}!C{r}= {r} - ALPHA*{q}do i= 1, N
PHI(i)= PHI(i) + ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo

1

2

2

3

3

4

4
5

1

6

2

7

3

8

4

9

1

10

2

11

3

12

4

Intro pFEM

25

内積
全体で和をとる必要がある⇒集団通信
1!C!C-- ALPHA= RHO / {p}{q}C1= 0.d0
do i= 1, N
C1= C1 + W(i,P)*W(i,Q)
enddo
ALPHA= RHO / C1

1

2

2

3

3

4

4
5

1

6

2

7

3

8

4

9

1

10

2

11

3

12

4

S2-ref26

行列ベクトル積
外点の値が必要⇒1対1通信!C!C-- {q}= [A]{p}do i= 1, N
W(i,Q) = DIAG(i)*W(i,P)
do j= INDEX(i-1)+1, INDEX(i)W(i,Q) = W(i,Q) + AMAT(j)*W(ITEM(j),P)
enddo
enddo

5

1

2

3

4

6

Intro pFEM

27

一次元問題：11要素，12節点，3領域
内点・外点・境界点
境界点は，他の領域の外点

#0

1

2
1

3
3

2

#1

X

4

5
4

5

1
4

2
1

3
2

4
3

外点：External Nodes
#2

Y

6
5

境界点：Boundary Nodes

5

1
4

2
1

3
2

4
3

S2-ref28

行列ベクトル積：ローカルに計算実施可能
1
2
3
4
5
6

1

1

2

2

3

3

4

4

5

5

6
7
7
9
10
11
12

=

6

7

7

8

8

9

9

10

10

11

11

12

12

29

係数行列が疎であるため，局所行列の和集合が全体行列となる
1

1

1

1

2

2

2

2

3

3

#0

3

2
4

4

4

4

4

4

5

5

5

5

5
6

6

6

6
7

7

7

#1

7
8

8
8

8

8

8
9

9

9

9

9
10
10
11
11
12

10

10
11

11
12

12

#2

S2-ref30

行列ベクトル積：ローカルに計算実施可能
1
2
3
4
5
6

1

1

2

2

3

3

4

4

5

5

6
7
8
9
10
11
12

=

6

7

7

8

8

9

9

10

10

11

11

12

12

S2-ref31

行列ベクトル積：ローカルに計算実施可能
1
2
3
4
5
6

1

1

2

2

3

3

4

4

5

5

6
7
8
9
10
11
12

=

6

7

7

8

8

9

9

10

10

11

11

12

12

S2-ref32

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

1
2

2
3
4

1
2

4

=

2

3

3

4

4

1

1

2
3

1

=

2

3

3

4

4

5
6

5

1

2

3

4

6

Intro pFEM

33

各領域データ（局所データ）仕様
10

9

8

11

4

12

6

5

• 内点，外点（internal/external
nodes）– 内点～外点となるように局所番号をつける• 隣接領域情報
7

1

2

3

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

34

Intro pFEM

Boundary Nodes（境界点）: SEND
PE#2 : 境界点（Boundary Nodes）の値を送信
15

6

7

PE#0

11

10

5

12

13

4

5

10

1

2

3

8

9

11

12

10

9

11

12

9

6

3

14

8

4

8

4

6

5

PE#3

PE#2
1

2

7

7

1

2

3

S2-ref

35

SEND: MPI_Isend/Irecv/Waitall
SENDbuf

neib#1neib#2neib#3neib#4BUFlength_e

BUFlength_e

BUFlength_e

BUFlength_e

export_index(0)+1

export_index(1)+1 export_index(2)+1

export_index(3)+1

export_index(4)do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= VAL(kk)enddo
Copied to sending buffers
enddo
do neib= 1, NEIBPETOTiS_e= export_index(neib-1) + 1
iE_e= export_index(neib )BUFlength_e= iE_e + 1 - iS_e
call MPI_ISEND&&(SENDbuf(iS_e), BUFlength_e, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_send(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_send, stat_send, ierr)36

Intro pFEM

External Nodes（外点）: RECEIVE
PE#2 : 外点（External Nodes）の値を受信
15

6

7

PE#0

11

10

5

12

13

4

5

10

1

2

3

8

9

11

12

10

9

11

12

9

6

3

14

8

4

8

4

6

5

PE#2

PE#3
1

2

7

7

1

2

3

S2-ref37

RECV: MPI_Isend/Irecv/Waitall
do neib= 1, NEIBPETOTiS_i= import_index(neib-1) + 1
iE_i= import_index(neib )BUFlength_i= iE_i + 1 - iS_i
call MPI_IRECV&&(RECVbuf(iS_i), BUFlength_i, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_recv(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_recv, stat_recv, ierr)do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)VAL(kk)= RECVbuf(k)enddo
Copied from receiving buffer
enddo

RECVbuf

neib#1neib#2neib#3neib#4BUFlength_i

BUFlength_i

BUFlength_i

BUFlength_i

import_index(0)+1

import_index(1)+1 import_index(2)+1

import_index(3)+1

import_index(4)S2-ref

プログラム：1d.f（1/11）
諸変数
program heat1Dp
implicit REAL*8 (A-H,O-Z)
include 'mpif.h'
integer :: N, NPLU, ITERmax
integer :: R, Z, P, Q, DD
real(kind=8) :: dX, RESID, EPS
real(kind=8) :: AREA, QV, COND
real(kind=8), dimension(:), allocatable :: PHI, RHS
real(kind=8), dimension(: ), allocatable :: DIAG, AMAT
real(kind=8), dimension(:,:), allocatable :: Wreal(kind=8), dimension(2,2) :: KMAT, EMAT
integer, dimension(:), allocatable :: ICELNOD
integer, dimension(:), allocatable :: INDEX, ITEM
integer(kind=4) :: NEIBPETOT, BUFlength, PETOT
integer(kind=4), dimension(2) :: NEIBPE
integer(kind=4), dimension(0:2) :: import_index, export_index
integer(kind=4), dimension( 2) :: import_item , export_item
real(kind=8), dimension(2) :: SENDbuf, RECVbuf
integer(kind=4), dimension(:,:), allocatable :: stat_send
integer(kind=4), dimension(:,:), allocatable :: stat_recv
integer(kind=4), dimension(: ), allocatable :: request_send
integer(kind=4), dimension(: ), allocatable :: request_recv38

S2-ref

プログラム：1d.f（2/11）
制御データ読み込み!C!C +-------+!C | INIT. |!C +-------+!C===!C!C-- MPI init.
call MPI_Init(ierr)
call MPI_Comm_size (MPI_COMM_WORLD, PETOT, ierr )call MPI_Comm_rank (MPI_COMM_WORLD, my_rank, ierr )MPI初期化：必須
全プロセス数：PETOT
自分のランク番号（0～PETOT-1）：my_rank!C!C-- CTRL data
if (my_rank.eq.0) thenopen (11, file='input.dat', status='unknown')
read (11,*) NEg
read (11,*) dX, QV, AREA, COND
read (11,*) ITERmax
read (11,*) EPS
close (11)endifcall MPI_Bcast (NEg
, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (ITERmax, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (dX, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (QV, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (AREA, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (COND, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (EPS
, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)

39

S2-ref

プログラム：1d.f（2/11）
制御データ読み込み!C!C +-------+!C | INIT. |!C +-------+!C===!C!C-- MPI init.
call MPI_Init(ierr)
call MPI_Comm_size (MPI_COMM_WORLD, PETOT, ierr )call MPI_Comm_rank (MPI_COMM_WORLD, my_rank, ierr )!C!C-- CTRL data
if (my_rank.eq.0) thenopen (11, file='input.dat', status='unknown')
read (11,*) Neg
read (11,*) dX, QV, AREA, COND
read (11,*) ITERmax
read (11,*) EPS
close (11)endif

MPI初期化：必須
全プロセス数：PETOT
自分のランク番号（0～PETOT-1）：my_rank

my_rank=0 のとき制御データを読み込む
Neg: 「全」要素数

call MPI_Bcast (NEg
, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (ITERmax, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (dX, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (QV, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (AREA, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (COND, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (EPS
, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)40

S2-ref

プログラム：1d.f（2/11）
制御データ読み込み!C!C +-------+!C | INIT. |!C +-------+!C===!C!C-- MPI init.
call MPI_Init(ierr)
call MPI_Comm_size (MPI_COMM_WORLD, PETOT, ierr )call MPI_Comm_rank (MPI_COMM_WORLD, my_rank, ierr )!C!C-- CTRL data
if (my_rank.eq.0) thenopen (11, file='input.dat', status='unknown')
read (11,*) Neg
read (11,*) dX, QV, AREA, COND
read (11,*) ITERmax
read (11,*) EPS
close (11)endif

MPI初期化：必須
全プロセス数：PETOT
自分のランク番号（0～PETOT-1）：my_rank

my_rank=0 のとき制御データを読み込む
Neg: 「全」要素数

call MPI_Bcast (NEg
, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) 0番から各プロセスにデータ送信
call MPI_Bcast (ITERmax, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (dX, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (QV, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (AREA, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (COND, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)call MPI_Bcast (EPS
, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)

41

S2-ref42

MPI_BCASTP#0 A0 B0 C0 D0
P#1

P#0 A0 B0 C0 D0BroadcastP#1 A0 B0 C0 D0

P#2

P#2 A0 B0 C0 D0

P#3

P#3 A0 B0 C0 D0• コミュニケーター 「comm」内の一つの送信元プロセス「root」のバッファ「buffer」から，その他全てのプロセスのバッファ「buffer」にメッセージを送信。
• call MPI_BCAST (buffer,count,datatype,root,comm,ierr)– buffer

任意I/O– count

整数I

バッファの先頭アドレス，
タイプは「datatype」により決定
メッセージのサイズ– datatype 整数

I

メッセージのデータタイプ

FORTRANC– root– comm– ierr

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

整数
整数
整数

I
I
O

送信元プロセスのID（ランク）
コミュニケータを指定する
完了コード

S2-ref

プログラム：1d.f（3/11）
局所分散メッシュデータ!C!C-- Local Mesh Size
Ng= NEg + 1
N = Ng / PETOT

総節点数
局所節点数

nr = Ng - N*PETOT
if (my_rank.lt.nr) N= N+1

NgがPETOTで割り切れない場合

NE= N - 1 + 2
NP= N + 2
if (my_rank.eq.0) NE= N - 1 + 1
if (my_rank.eq.0) NP= N + 1
if (my_rank.eq.PETOT-1) NE= N - 1 + 1
if (my_rank.eq.PETOT-1) NP= N + 1
if (PETOT.eq.1) NE= N-1
if (PETOT.eq.1) NP= N!C!C– ARRAYs
allocate (PHI(NP), DIAG(NP), AMAT(2*NP-2), RHS(NP))
allocate (ICELNOD(2*NE))
allocate (INDEX(0:NP), ITEM(2*NP-2), W(NP,4))
PHI= 0.d0
AMAT= 0.d0
DIAG= 0.d0
RHS= 0.d043

S2-ref

44

プログラム：1d.f（3/11）
局所分散メッシュデータ，各要素→一様!C!C-- Local Mesh Size
Ng= NEg + 1
N = Ng / PETOT

総節点数
局所節点数

nr = Ng - N*PETOT
if (my_rank.lt.nr) N= N+1

NgがPETOTで割り切れない場合

NE= N - 1 + 2
NP= N + 2

NE：局所要素数
NP：局所節点数（内点＋外点）

if (my_rank.eq.0) NE= N - 1 + 1
if (my_rank.eq.0) NP= N + 1
if (my_rank.eq.PETOT-1) NE= N - 1 + 1
if (my_rank.eq.PETOT-1) NP= N + 1
if (PETOT.eq.1) NE= N-1
if (PETOT.eq.1) NP= N

N+1

1
N

2
1

…
2

!C
!C– ARRAYs
allocate (PHI(NP), DIAG(NP), AMAT(2*NP-2), RHS(NP))
allocate (ICELNOD(2*NE))
allocate (INDEX(0:NP), ITEM(2*NP-2), W(NP,4))
PHI= 0.d0
AMAT= 0.d0
DIAG= 0.d0
RHS= 0.d0

N
N-1

N+2
N+1

一般の領域：
N+2節点，N+1要素

S2-ref

プログラム：1d.f（3/11）
局所分散メッシュデータ，各要素→一様!C!C-- Local Mesh Size
Ng= NEg + 1
N = Ng / PETOT

総節点数
局所節点数

nr = Ng - N*PETOT
if (my_rank.lt.nr) N= N+1

NgがPETOTで割り切れない場合

NE= N - 1 + 2
NP= N + 2

NE：局所要素数
NP：局所節点数（内点＋外点）

if (my_rank.eq.0) NE= N - 1 + 1
if (my_rank.eq.0) NP= N + 1
if (my_rank.eq.PETOT-1) NE= N - 1 + 1
if (my_rank.eq.PETOT-1) NP= N + 1
if (PETOT.eq.1) NE= N-1
if (PETOT.eq.1) NP= N

1

2
1

…
2

!C
!C– ARRAYs
allocate (PHI(NP), DIAG(NP), AMAT(2*NP-2), RHS(NP))
allocate (ICELNOD(2*NE))
allocate (INDEX(0:NP), ITEM(2*NP-2), W(NP,4))
PHI= 0.d0
AMAT= 0.d0
DIAG= 0.d0
RHS= 0.d0

N
N-1

N+1
N

#0: N+1節点，N要素

45

S2-ref

プログラム：1d.f（3/11）
局所分散メッシュデータ，各要素→一様!C!C-- Local Mesh Size
Ng= NEg + 1
N = Ng / PETOT

総節点数
局所節点数

nr = Ng - N*PETOT
if (my_rank.lt.nr) N= N+1

NgがPETOTで割り切れない場合

NE= N - 1 + 2
NP= N + 2

NE：局所要素数
NP：局所節点数（内点＋外点）

if (my_rank.eq.0) NE= N - 1 + 1
if (my_rank.eq.0) NP= N + 1
if (my_rank.eq.PETOT-1) NE= N - 1 + 1
if (my_rank.eq.PETOT-1) NP= N + 1
if (PETOT.eq.1) NE= N-1
N+1
if (PETOT.eq.1) NP= N

1
N

2
1

…
2

N
N-1

!C
!C– ARRAYs
allocate (PHI(NP), DIAG(NP), AMAT(2*NP-2), RHS(NP))
allocate (ICELNOD(2*NE))
allocate (INDEX(0:NP), ITEM(2*NP-2), W(NP,4))
PHI= 0.d0
AMAT= 0.d0
DIAG= 0.d0
RHS= 0.d0#PETot-1: N+1節点，N要素

46

S2-ref

プログラム：1d.f（3/11）
局所分散メッシュデータ!C!C-- Local Mesh Size
Ng= NEg + 1
N = Ng / PETOT

総節点数
局所節点数

nr = Ng - N*PETOT
if (my_rank.lt.nr) N= N+1

NgがPETOTで割り切れない場合

NE= N - 1 + 2
NP= N + 2

NE：局所要素数
NP：局所節点数（内点＋外点）

if (my_rank.eq.0) NE= N - 1 + 1
if (my_rank.eq.0) NP= N + 1
if (my_rank.eq.PETOT-1) NE= N - 1 + 1
if (my_rank.eq.PETOT-1) NP= N + 1
if (PETOT.eq.1) NE= N-1
if (PETOT.eq.1) NP= N!C!C– ARRAYs
allocate (PHI(NP), DIAG(NP), AMAT(2*NP-2), RHS(NP))
allocate (ICELNOD(2*NE))
allocate (INDEX(0:NP), ITEM(2*NP-2), W(NP,4))
PHI= 0.d0
AMAT= 0.d0
DIAG= 0.d0
RHS= 0.d0

NでなくNPで配列を定義している点に注意

47

S2-ref

プログラム：1d.f（4/11）
配列初期化，要素～節点
do icel= 1, NE
ICELNOD(2*icel-1)= icel
ICELNOD(2*icel )= icel + 1
enddo
if (PETOT.gt.1) thenif (my_rank.eq.0) thenicel= NEICELNOD(2*icel-1)= N
ICELNOD(2*icel )= N + 1
else if (my_rank.eq.PETOT-1) thenicel= NEICELNOD(2*icel-1)= N + 1
ICELNOD(2*icel )= 1
else
icel= NE - 1
ICELNOD(2*icel-1)= N + 1
ICELNOD(2*icel )= 1
icel= NEICELNOD(2*icel-1)= N
ICELNOD(2*icel )= N + 2
endif
endif

icel
Icelnod(2*icel-1)=icel

Icelnod(2*icel)=icel+148

S2-ref

プログラム：1d.f（4/11）
配列初期化，要素～節点
do icel= 1, NE
ICELNOD(2*icel-1)= icel
ICELNOD(2*icel )= icel + 1
enddo「1-2」の要素を「1」とする

if (PETOT.gt.1) thenif (my_rank.eq.0) thenicel= NEICELNOD(2*icel-1)= N
ICELNOD(2*icel )= N + 1

12
1

else if (my_rank.eq.PETOT-1) thenicel= NEICELNOD(2*icel-1)= N + 1
N+11
ICELNOD(2*icel )= 1
N

else
icel= NE - 1
ICELNOD(2*icel-1)= N + 1
ICELNOD(2*icel )= 1
N+1
icel= NEN
ICELNOD(2*icel-1)= N
ICELNOD(2*icel )= N + 2
endif
endif

…

N
N-1

2

2
1

1

…
2

2
1

N+1
N

N

#PETot-1: N+1節点，N要素

N-1

…
2

#0: N+1節点，N要素

N
N-1

N+2
N+1

一般の領域：
N+2節点，N+1要素

49

S2-ref

プログラム：1d.f（5/11）
Index定義

KMAT(1,1)= +1.d0
KMAT(1,2)= -1.d0
KMAT(2,1)= -1.d0
KMAT(2,2)= +1.d0!C===!C!C +--------------+!C | CONNECTIVITY |!C +--------------+!C===
INDEX = 2

1
1
N+1

1
N

N+1

2

1
N

INDEX(N+1)= 1
INDEX(NP )= 1
if (my_rank.eq.0)
INDEX(1)= 1
if (my_rank.eq.PETOT-1) INDEX(N)= 1
do i= 1, NP
INDEX(i)= INDEX(i) + INDEX(i-1)
enddo

…

…
2

2
1

N
N-1

2

1

INDEX(0)= 0

NPLU= INDEX(NP)ITEM= 0

2

N
N

N
N-1

#0: N+1節点，N要素#PETot-1: N+1節点，N要素

N-1

…
2

N+1

N+2
N+1

一般の領域：
N+2節点，N+1要素

50

S2-ref

プログラム：1d.f（6/11）
Item定義do i= 1, N
jS= INDEX(i-1)
if (my_rank.eq.0.and.i.eq.1) then
ITEM(jS+1)= i+1
else if (my_rank.eq.PETOT-1.and.i.eq.N) then
ITEM(jS+1)= i-1
else
ITEM(jS+1)= i-1
ITEM(jS+2)= i+1
if (i.eq.1) ITEM(jS+1)= N + 1
if (i.eq.N) ITEM(jS+2)= N + 2
if (my_rank.eq.0.and.i.eq.N) ITEM(jS+2)= N + 1
endif
enddo
1

2
1

i = N + 1
jS= INDEX(i-1)
if (my_rank.eq.0) then
ITEM(jS+1)= N
else
ITEM(jS+1)= 1
endif

N+1

1
N

N+1

2

…
2

2
1

N
N-1

2

1

1
N

…

N

i = N + 2
if (my_rank.ne.0.and.my_rank.ne.PETOT-1) thenjS= INDEX(i-1)
ITEM(jS+1)= N
endif

N
N-1

#0: N+1節点，N要素#PETot-1: N+1節点，N要素

N-1

…
2

N+1
N

N+2
N+1

一般の領域：
N+2節点，N+1要素

51

S2-ref

プログラム：1d.f（7/11）
通信テーブル定義!C!C-- COMMUNICATIONNEIBPETOT= 2
if (my_rank.eq.0) NEIBPETOT= 1
if (my_rank.eq.PETOT-1) NEIBPETOT= 1
if (PETOT.eq.1)
NEIBPETOT= 0
NEIBPE(1)= my_rank - 1
NEIBPE(2)= my_rank + 1
if (my_rank.eq.0) NEIBPE(1)= my_rank + 1
if (my_rank.eq.PETOT-1) NEIBPE(1)= my_rank - 1
BUFlength= 1
import_index(1)= 1
import_index(2)= 2
import_item (1)= N+1
import_item (2)= N+2
export_index(1)= 1
export_index(2)= 2
export_item (1)= 1
export_item (2)= N
if (my_rank.eq.0) thenimport_item (1)= N+1
export_item (1)= N
endif

1

2
1

N+1

1
N

N+1

2

…
2

2
1

N
N-1

2

1

1
N

…

N

N
N-1

#0: N+1節点，N要素#PETot-1: N+1節点，N要素

N-1

…
2

N+1
N

N+2
N+1

一般の領域：
N+2節点，N+1要素!C!C-- INIT. arrays for MPI_Waitall
allocate (stat_send(MPI_STATUS_SIZE,NEIBPETOT), stat_recv(MPI_STATUS_SIZE,NEIBPETOT))
allocate (request_send(NEIBPETOT), request_recv(NEIBPETOT))52

S2-ref53

MPI_ISEND•

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

• call MPI_ISEND(sendbuf,count,datatype,dest,tag,comm,request, ierr)sendbuf 任意
count
整数
datatype 整数
dest
整数
tag
整数

I
I
I
I
I

– comm– request

整数
整数

I
O

– ierr

整数

O

–
–
–
–
–

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

S2-ref54

MPI_IRECV•

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

• call MPI_IRECV(recvbuf,count,datatype,dest,tag,comm,request, ierr)recvbuf 任意
count
整数
datatype 整数
dest
整数
tag
整数

I
I
I
I
I

– comm– request

整数
整数

I
O

– ierr

整数

O

–
–
–
–
–

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

S2-ref55

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」を呼び出すと処理が遅くなる，というような経験もある。

• call MPI_WAITALL(count,request,status,ierr)– count– request

整数
整数

I
I/O

– status

整数O– ierr

整数O

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

S2-ref56

一般化された通信テーブル：送信• 送信相手– NEIBPETOT，NEIBPE(neib)• それぞれの送信相手に送るメッセージサイズ– export_index(neib), neib= 0, NEIBPETOT• 「境界点」番号– export_item(k), k= 1, export_index(NEIBPETOT)• それぞれの送信相手に送るメッセージ– SENDbuf(k), k= 1, export_index(NEIBPETOT)S2-ref57

送信（MPI_Isend/Irecv/Waitall）
SENDbuf

neib#1neib#2neib#3neib#4BUFlength_e

BUFlength_e

BUFlength_e

BUFlength_e

export_index(0)+1

export_index(1)+1 export_index(2)+1

export_index(3)+1

export_index(4)do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= VAL(kk)送信バッファへの代入
enddo
enddo
温度などの変数を直接送信，受信に使
do neib= 1, NEIBPETOTiS_e= export_index(neib-1) + 1
iE_e= export_index(neib )BUFlength_e= iE_e + 1 - iS_e

うのではなく，このようなバッファへ一回代入して計算することを勧める。

call MPI_ISEND&&(SENDbuf(iS_e), BUFlength_e, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_send(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_send, stat_recv, ierr)S2-ref58

送信：一次元問題• 送信相手

N+1

1

2

N-1

SENDbuf(1)= VAL(1)• NEIBPETOT=2, NEIB(1)= my_rank-1, NEIB(2)= my_rank+1• それぞれの送信相手に送るメッセージサイズ– export_index(neib), neib= 0, NEIBPETOT• export_index(0)=0, export_index(1)= 1, export_index(2)= 2– export_item(k), k= 1, export_index(NEIBPETOT)• export_item(1)= 1, export_item(2)= N• それぞれの送信相手に送るメッセージ– SENDbuf(k), k= 1, export_index(NEIBPETOT)• SENDbuf(1)= BUF(1), SENDbuf(2)= BUF(N)

N+2

SENDbuf(2)= VAL(N)– NEIBPETOT，NEIBPE(neib)• 「境界点」番号N

S2-ref59

一般化された通信テーブル：受信• 受信相手– NEIBPETOT，NEIBPE(neib)• それぞれの受信相手から受け取るメッセージサイズ– import_index(neib), neib= 0, NEIBPETOT• 「外点」番号– import_item(k), k= 1, import_index(NEIBPETOT)• それぞれの受信相手から受け取るメッセージ– RECVbuf(k), k= 1, import_index(NEIBPETOT)S2-ref60

受信（MPI_Isend/Irecv/Waitall）
do neib= 1, NEIBPETOTiS_i= import_index(neib-1) + 1
iE_i= import_index(neib )BUFlength_i= iE_i + 1 - iS_i
call MPI_IRECV&&(RECVbuf(iS_i), BUFlength_i, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_recv(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_recv, stat_recv, ierr)do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)VAL(kk)= RECVbuf(k)
enddo
enddo

受信バッファから代入

RECVbuf

neib#1neib#2neib#3neib#4BUFlength_i

BUFlength_i

BUFlength_i

BUFlength_i

import_index(0)+1

import_index(1)+1 import_index(2)+1

import_index(3)+1

import_index(4)S2-ref61

受信:一次元問題• 受信相手

N+1

1

2

N-1

VAL(N+1)=
RECVbuf(1)N

VAL(N+2)=
RECVbuf(2)– NEIBPETOT，NEIBPE(neib)• NEIBPETOT=2, NEIB(1)= my_rank-1, NEIB(2)= my_rank+1• それぞれの受信相手から受け取るメッセージサイズ– import_index(neib), neib= 0, NEIBPETOT• import_index(0)=0, import_index(1)= 1, import_index(2)= 2• 「外点」番号– import_item(k), k= 1, import_index(NEIBPETOT)• import_item(1)= N+1, import_item(2)= N+2• それぞれの受信相手から受け取るメッセージ– RECVbuf(k), k= 1, import_index(NEIBPETOT)• BUF(N+1)=RECVbuf(1), BUF(N+2)=RECVbuf(2)N+2

S2-ref62

一般化された通信テーブル：Fortran
N+1

1

2

SENDbuf(1)= VAL(1)N+1

1

VAL(N+1)=
RECVbuf(1)

2

N-1

N

N+2

SENDbuf(2)= VAL(N)

N-1

N

N+2

VAL(N+2)=
RECVbuf(2)NEIBPETOT= 2
NEIBPE(1)= my_rank - 1
NEIBPE(2)= my_rank + 1
import_index(1)= 1
import_index(2)= 2
import_item (1)= N+1
import_item (2)= N+2
export_index(1)= 1
export_index(2)= 2
export_item (1)= 1
export_item (2)= N
if (my_rank.eq.0) thenimport_item (1)= N+1
export_item (1)= N
NEIBPE(1)= my_rank+1
endifS2-ref63

一般化された通信テーブル：C言語
N

0

1

N-2 N-1 N+1

SENDbuf[0]= VAL[0]SENDbuf[1]= VAL[N-1]0

N-2 N-1 N+1N

VAL[N]=
RECVbuf[0]1

VAL[N+1]=
RECVbuf[1]NEIBPETOT= 2
NEIBPE[0]= my_rank - 1
NEIBPE[1]= my_rank + 1
import_index[1]= 1
import_index[2]= 2
import_item [0]= N
import_item [1]= N+1
export_index[1]= 1
export_index[2]= 2
export_item [0]= 0
export_item [1]= N-1
if (my_rank.eq.0) thenimport_item [0]= N
export_item [0]= N-1
NEIBPE[0]= my_rank+1
endifS2-ref64

プログラム：1d.f（8/11）
全体マトリクス生成：1CPUのときと全く同じ：各要素→一様!C!C +-----------------+!C | MATRIX ASSEMBLE |!C +-----------------+#0!C===
do icel= 1, NE
in1= ICELNOD(2*icel-1)
in2= ICELNOD(2*icel )DL = dX
cK= AREA*COND/DL
EMAT(1,1)= Ck*KMAT(1,1)
EMAT(1,2)= Ck*KMAT(1,2)
EMAT(2,1)= Ck*KMAT(2,1)
EMAT(2,2)= Ck*KMAT(2,2)

1

2
1

3

DIAG(in1)= DIAG(in1) + EMAT(1,1)
DIAG(in2)= DIAG(in2) + EMAT(2,2)
if (my_rank.eq.0.and.icel.eq.1) thenk1= INDEX(in1-1) + 1
else
k1= INDEX(in1-1) + 2
endif
k2= INDEX(in2-1) + 1
AMAT(k1)= AMAT(k1) + EMAT(1,2)
AMAT(k2)= AMAT(k2) + EMAT(2,1)!C===

QN= 0.50d0*QV*AREA*DL
RHS(in1)= RHS(in1) + QN
RHS(in2)= RHS(in2) + QN
enddo

4
3

2

#1

5
4

5

1
4

2
1

3
2

4
3

#2

6
5

5

1
4

2
1

3
2

4
3

FEM3D

65

Local Matrix：各プロセスにおける係数行列

NP
N

NP
N

NP
N

internal

external

FEM3D66

本当に必要なのはこの部分

NP
N

NP
N

NP
N

internal

external

pFEM3D-2

67

MAT_ASS_MAIN: Overview
Gaussian Quad. points in ζ-direction
Gaussian Quad. points in η-direction
Gaussian Quad. Pointe in ξ-direction
Define Shape Function at Gaussian Quad. Points (8-points)
Its derivative on natural/local coordinate is also defined.

do kpn= 1, 2
do jpn= 1, 2
do ipn= 1, 2
enddo
enddo
enddo

do icel= 1, ICELTOT
Loop for Element
Jacobian and derivative on global coordinate of shape functions at
Gaussian Quad. Points are defined according to coordinates of 8 nodes.（JACOBI）
do ie= 1, 8
do je= 1, 8

Local Node ID
Local Node ID
Global Node ID: ip, jp
Address of Aip,jp in “item”: kk

do kpn= 1, 2
Gaussian Quad. points in ζ-direction
do jpn= 1, 2
Gaussian Quad. points in η-direction
do ipn= 1, 2 Gaussian Quad. points in ξ-direction
integration on each element
coefficients of element matrices
accumulation to global matrix
enddo
enddo
enddo
enddo
enddo
enddo

je

ie

pFEM3D-2

68

全ての要素の計算を実施する
外点を含むオーバーラップ領域の要素の計算も実施
PE#1
PE#1

21

22

4

PE#0

23

5

17

18

12

11

7

6

19

13

15

9

２
PE#3

３

４
PE#2

6

7

PE#0

2

11

3

10

7

8

9

11

10

12

5

10

14

13

4

5

10

1

2

3

8

9

11

12

10

9

11

12

9

6

3
１

15

20

14

8

12

25

24

1
16

6

8

4

8

4

6

5

5
1

2
PE#3

7

7

1

2
PE#2

3

FEM3D69

従って結果的にはこのような行列を得るが

NP
N

NP
N

NP
N

internal

external

FEM3D70

黒枠で囲んだ部分の行列は不完全
しかし，計算には使用しないのでこれで良い

NP
N

NP
N

NP
N

internal

external

S2-ref

71

プログラム：1d.f（9/11）
境界条件：1CPUのときとほとんど同じ!C!C +---------------------+!C | BOUNDARY CONDITIONS |!C +---------------------+
!C===

#0

1

2
1

!C
!C-- X=Xmin
if (my_rank.eq.0) theni = 1
jS= INDEX(i-1)
AMAT(jS+1)= 0.d0
DIAG(i)= 1.d0
RHS (i)= 0.d0
do k= 1, NPLU
if (ITEM(k).eq.1) AMAT(k)= 0.d0
enddo
endif
!C===

3

4
3

2

#1

5
4

5

1
4

2
1

3
2

4
3

#2

6
5

5

1
4

2
1

3
2

4
3

S2-ref

プログラム：1d.f（10/11）
共役勾配法

!C
!C +---------------+!C | CG iterations |!C +---------------+
!C===
R = 1
Z = 2
Q = 2
P = 3
DD= 4
do i= 1, N
W(i,DD)= 1.0D0 / DIAG(i)enddo!C!C-- {r0}= {b} - [A]{xini} |
!C-

init
do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= PHI(kk)
enddo
enddo

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|end72

S2-ref73

共役勾配法• 行列ベクトル積• 内積• 前処理：1CPUのときと同じ• DAXPY：1CPUのときと同じ

S2-ref

74

前処理，DAXPY!C!C-- {z}= [Minv]{r}do i= 1, N
W(i,Z)= W(i,DD) * W(i,R)
enddo!C!C-- {x}= {x} + ALPHA*{p}!C{r}= {r} - ALPHA*{q}do i= 1, N
PHI(i)= PHI(i) + ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddoS2-ref75

行列ベクトル積（1/2）
通信テーブル使用，{p}の最新値を計算前に取得!C!C-- {q}= [A]{p}do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= W(kk,P)
enddo
enddo
do neib= 1, NEIBPETOT
is= export_index(neib-1) + 1
len_s= export_index(neib) - export_index(neib-1)
call MPI_Isend (SENDbuf(is), len_s, MPI_DOUBLE_PRECISION,&&
NEIBPE(neib), 0, MPI_COMM_WORLD, request_send(neib), ierr)enddo
do neib= 1, NEIBPETOT
ir= import_index(neib-1) + 1
len_r= import_index(neib) - import_index(neib-1)
call MPI_Irecv (RECVbuf(ir), len_r, MPI_DOUBLE_PRECISION,&&
NEIBPE(neib), 0, MPI_COMM_WORLD, request_recv(neib), ierr)enddo
call MPI_Waitall (NEIBPETOT, request_recv, stat_recv, ierr)do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)W(kk,P)= RECVbuf(k)
enddo
enddo

S2-ref76

行列ベクトル積（2/2）{q}= [A]{p}call MPI_Waitall (NEIBPETOT, request_send, stat_send, ierr)
do i= 1, N
W(i,Q) = DIAG(i)*W(i,P)
do j= INDEX(i-1)+1, INDEX(i)W(i,Q) = W(i,Q) + AMAT(j)*W(ITEM(j),P)
enddo
enddo

S2-ref77

内積
各プロセスで計算した値を，MPI_Allreduceで合計!C!C-- RHO= {r}{z}RHO0= 0.d0
do i= 1, N
RHO0= RHO0 + W(i,R)*W(i,Z)
enddo
call MPI_Allreduce (RHO0, RHO, 1, MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)&S2-ref

78
P#0 A0 B0 C0 D0

MPI_ALLREDUCE••

P#1 A1 B1 C1 D1P#0 op.A0-A3 op.B0-B3 op.C0-C3 op.D0-D3
All 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
総和，最大値を計算したら，各プロセスで利用したい場合が多い• call MPI_ALLREDUCE(sendbuf,recvbuf,count,datatype,op, comm,ierr)
任意
任意

I
O

count
整数
datatype 整数
op
整数
comm
整数
ierr
整数

I
I
I
I
O

– sendbuf– recvbuf
–
–
–
–
–

送信バッファの先頭アドレス，
受信バッファの先頭アドレス，
タイプは「datatype」により決定
メッセージのサイズ
メッセージのデータタイプ
計算の種類
コミュニケータを指定する
完了コード

S2-ref79

CG法（1/5）!C!C-- {r0}= {b} - [A]{xini}
do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= PHI(kk)
enddo
enddo

&
&
&

&
&
&

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=1do neib= 1, NEIBPETOT
p(1)= z(0)is = export_index(neib-1) + 1
len_s= export_index(neib) - export_index(neib-1)
else
call MPI_Isend (SENDbuf(is), len_s,
βi-1= ρi-1/ρi-2
MPI_DOUBLE_PRECISION,
NEIBPE(neib), 0, MPI_COMM_WORLD,& (i)(i-1) + β
(i-1)
p
=
z
p
i-1
request_send(neib), ierr)
enddo
endif
q(i)= [A]p(i)do neib= 1, NEIBPETOT(i)q(i)ir = import_index(neib-1) + 1
α
=
ρ
/p
i
i-1
len_r= import_index(neib) - import_index(neib-1)
x(i)= x(i-1) + αip(i)call MPI_Irecv (RECVbuf(ir), len_r,
MPI_DOUBLE_PRECISION,(i-1) - α q(i)r(i)
iNEIBPE(neib), 0, MPI_COMM_WORLD,& = r
request_recv(neib), ierr)check convergence |r|enddo
call MPI_Waitall (NEIBPETOT, request_recv, stat_recv, ierr)end
do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)PHI(kk)= RECVbuf(k)
enddo
enddo
call MPI_Waitall (NEIBPETOT, request_send, stat_send, ierr)S2-ref80

CG法（2/5）
do i= 1, N
W(i,R) = DIAG(i)*PHI(i)do j= INDEX(i-1)+1, INDEX(i)W(i,R) = W(i,R) + AMAT(j)*PHI(ITEM(j))
enddo
enddo

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)
BNRM20= 0.0D0
if i=1
do i= 1, N
p(1)= z(0)BNRM20 = BNRM20 + RHS(i) **2
W(i,R) = RHS(i) - W(i,R)
else
enddo
βi-1= ρi-1/ρi-2
call MPI_Allreduce (BNRM20, BNRM2, 1,&
MPI_DOUBLE_PRECISION,
p(i)= z(i-1) + βi-1 p(i-1)&
MPI_SUM, MPI_COMM_WORLD, ierr)endif!C********************************************************************(i)q = [A]p(i)do iter= 1, ITERmax
αi = ρi-1/p(i)q(i)!C!C-- {z}= [Minv]{r}x(i)= x(i-1) + αip(i)r(i)= r(i-1) - αiq(i)do i= 1, N
W(i,Z)= W(i,DD) * W(i,R)
check convergence |r|
enddo
end
!C
!C-- RHO= {r}{z}RHO0= 0.d0
do i= 1, N
RHO0= RHO0 + W(i,R)*W(i,Z)
enddo
call MPI_Allreduce (RHO0, RHO, 1, MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)S2-ref81

CG法（3/5）!C!C-- {p} = {z} if
ITER=1!C BETA= RHO / RHO1 otherwise
if ( iter.eq.1 ) thendo i= 1, N
W(i,P)= W(i,Z)
enddo
else
BETA= RHO / RHO1
do i= 1, N
W(i,P)= W(i,Z) + BETA*W(i,P)
enddo
endif
!C
!C-- {q}= [A]{p}do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= W(kk,P)
enddo
enddo

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

do neib= 1, NEIBPETOT
is= export_index(neib-1) + 1
len_s= export_index(neib) - export_index(neib-1)
call MPI_Isend (SENDbuf(is), len_s, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD,&
request_send(neib), ierr)
enddo

&
&

S2-ref82

CG法（4/5）do neib= 1, NEIBPETOT
ir = import_index(neib-1) + 1
Compute r(0)= b-[A]x(0)len_r= import_index(neib) - import_index(neib-1)
call MPI_Irecv (RECVbuf(ir), len_r,
for i= 1, 2, …
&
MPI_DOUBLE_PRECISION,(i-1)= r(i-1)&
NEIBPE(neib), 0, MPI_COMM_WORLD,&solve [M]z&
request_recv(neib), ierr)ρi-1= r(i-1) z(i-1)
enddo
if i=1
call MPI_Waitall (NEIBPETOT, request_recv, stat_recv, ierr)p(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|do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)W(kk,P)= RECVbuf(kk)
enddo
enddo
call MPI_Waitall (NEIBPETOT, request_send, stat_send, ierr)
do i= 1, N
W(i,Q) = DIAG(i)*W(i,P)
do j= INDEX(i-1)+1, INDEX(i)W(i,Q) = W(i,Q) + AMAT(j)*W(ITEM(j),P)
enddo
enddo
!C
!C-- ALPHA= RHO / {p}{q}endC10= 0.d0
do i= 1, N
C10= C10 + W(i,P)*W(i,Q)
enddo
call MPI_Allreduce (C10, C1, 1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr)
ALPHA= RHO / C1

S2-ref83

CG法（5/5）!C!C-- {x}= {x} + ALPHA*{p}!C{r}= {r} - ALPHA*{q}do i= 1, N
PHI(i)= PHI(i) + ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo
DNRM20 = 0.0
do i= 1, N
DNRM20= DNRM20 + W(i,R)**2
enddo
call MPI_Allreduce (DNRM20, DNRM2, 1,&
MPI_DOUBLE_PRECISION,&MPI_SUM, MPI_COMM_WORLD, ierr)
RESID= dsqrt(DNRM2/BNRM2)
if (my_rank.eq.0.and.mod(iter,1000).eq.0) thenwrite (*, '(i8,1pe16.6)') iter, RESID
endif
if ( RESID.le.EPS) goto 900
RHO1 = RHO
enddo

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S2-ref

プログラム：1d.f（11/11）
結果書き出し：各プロセスごとに実施!C!C-- OUTPUT
if (my_rank.eq.0) thenwrite (*,'(2(1pe16.6))') E1Time-S1Time, E2Time-E1Time
endif
write (*,'(/a)') '### TEMPERATURE'
do i= 1, N
write (*,'(2i8, 2(1pe16.6))') my_rank, i, PHI(i)enddocall MPI_FINALIZE (ierr)end program heat1Dp

84

S2-ref• 問題の概要，実行方法• プログラムの説明
• 計算例

85

S2-ref86

6
計算結果（1次元）：CG法部分，10

1,000回反復に要する時間，Strong Scaling
1コア計算時間を基準，2ノード以上では56コア/ノード使用
5回計測，最速のものを採択
up to 1 node,56cores

up to 16 nodes,
896 cores

70

900

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

60

700
600

Speed-Up

Speed-Up50

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

800

40
30

500
400
300

20

200

10

100

0

0

0

10

20

30

Core #

40

50

60

0

100

200

300

400

500

Core #

600

700

800

900

S2-ref

1-node, 24-cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=24#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test.lst

24-cores are randomly
selected from 56-coreson the node
RANDOMmpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db
export I_MPI_PIN_PROCESSOR_LIST=0-23
mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

24-cores on Socket #0
are assigned.
SPECIFIED
Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores87

S2-ref

1-node, 28-cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=28#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o test.lst

28-cores are randomly
selected from 56-coreson the node
RANDOMmpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db
export I_MPI_PIN_PROCESSOR_LIST=0-27
mpiexec.hydra -n ${PJM_MPI_PROC} ./1da
mpiexec.hydra -n ${PJM_MPI_PROC} ./1db

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

28-cores on Socket #0
are assigned.
SPECIFIED
Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores88

S2-ref89

7
計算結果（1次元）：CG法部分，10

1,000回反復に要する時間，Strong Scaling
1コア計算時間を基準，2ノード以上では56コア/ノード使用
5回計測，最速のものを採択
up to 1 node,56cores

up to 16 nodes,
896 cores

70

900

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

60

700
600

Speed-Up

Speed-Up50

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

800

40
30

500
400
300

20

200
10

100
0

0
0

10

20

30

Core #

40

50

60

0

100

200

300

400

500

Core #

600

700

800

900

S2-ref90

大規模並列計算：理想値からのずれ• MPI通信そのものに要する時間– データを送付している時間– ノード間においては通信バンド幅によって決まる• Gigabit Ethernetでは 1Gbit/sec.（理想値）– 通信時間は送受信バッファのサイズに比例• MPIの立ち上がり時間– latency– 送受信バッファのサイズによらない• 呼び出し回数依存，プロセス数が増加すると増加する傾向– 通常，数～数十µsecのオーダー• MPIの同期のための時間– プロセス数が増加すると増加する傾向• 計算時間が小さい場合はこれらの効果を無視できない。
– 特に，送信メッセージ数が小さい場合は，「Latency」が効く。

S2-ref91

大規模並列計算：理想値からのずれ• 問題サイズが比較的小さいと，この効果は無視できない– 通信メッセージのサイズが小さいと，Latencyの効果大

1000
N=10^7
N=10^8N=10^9N=2x10^9
Ideal

Speed-UP100

1056-cores/node

1
1

10

100

Core #

1000

70
60
50

Speed-Up

Up to 56 cores(single node)Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

40
30
20

N=106

10
0
0

10

20

40

50

Core #

70

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

60
50

Speed-Up

30

N=107

40
30
20
10
0
0

10

20

30

Core #

40

50

92• ノード内（56コアまで）では問題
規模が小さいと性能良い• メモリ競合（通信の影響ではない）
60• ノード内メモリ転送性能は8コア
以上ではほとんど変化しない（Stream）• 問題サイズが小さいとキャッシュを有効利用できるため，メモリ転送性能の影響少ない• 2ソケット利用時のメモリ・キャ
シュは「Random」の方が有効に利用されている60

S2-ref93

問題サイズが小さいとキャッシュを有効利用できる• スカラープロセッサでは，全般に問題規模が小さいほど性能が高い。
– キャッシュの有効利用
3.00

CPU
Register
FAST

Performance of CG Solver
on IBM SP-3

Cache
(MB)

GFLOPS

2.50
2.00

SLOW

1.50
1.00

Main Memory
(GB)

0.50
0.00
1.0E+04

1.0E+05

1.0E+06

DOF: Problem Size

1.0E+07

Intro94

Category

Capacity

X-Way
Set Associative

L1$Data

32 KB/core

8-Way

64B

L1$Instruction

32 KB/core

8-Way

64B

L2

1.00 MB/core

16-Way64B

L338.5 MB/socket

11-Way64B

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Cache
Line

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

70
60
50

Speed-Up

Up to 56 cores(single node)Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

40
30

• Required Memory• N=106: 80MB

20

N=106

10

– Very Close to Cache Size

0
0

10

20

30

40

50

60

Core #

70

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

60
50

Speed-Up95

N=107• Memory throughput on eachsocket is constant for 8+ cores(Stream)• If problem size is small, cache canbe well-utilized. Therefore, effectof memory bandwidth is small.

40
30
20
10
0
0

10

20

30

Core #• N=107: 800MB

40

50

60

96

STREAM benchmark
http://www.cs.virginia.edu/stream/• Benchmarks for Memory Bandwidth– Copy:– Scale:– Add:– Triad:c(i)= a(i)c(i)= s*b(i)c(i)= a(i) + b(i)c(i)= a(i) + s*b(i)---------------------------------------------Double precision appears to have 16 digits of accuracy
Assuming 8 bytes per DOUBLE PRECISION word---------------------------------------------Number of processors =16
Array size =2000000
Offset=
0
The total memory requirement is732.4 MB(
45.8MB/task)
You are running each test 10 times-The *best* time for each test is used
*EXCLUDING* the first and last iterations------------------------------------------------------------------------------------------------------Function
Rate (MB/s) Avg time Min time Max timeCopy:
18334.1898
0.02800.02790.0280Scale:
18035.1690
0.02840.02840.0285
Add:
18649.4455
0.04120.04120.0413
Triad:
19603.8455
0.03940.03920.039897

マイクロプロセッサの動向
CPU性能，メモリバンド幅のギャップ

http://www.cs.virginia.edu/stream/

98

Sparse/Dense Matrices（疎・密）for (i=0;i<N;i++ {Y[i]= D[i]*X[i];for (k=index[i];k<index[i+1];k++) {Y[k]~ += AMAT[k]*X[item[k]];}}for (i=0;i<N;i++){for (j=0;j<N;j++) {Y[i] += A[i][j]*X[j];}}• “X” in RHS– 密行列：連続アクセス，キャッシュを有効利用– 疎行列：連続性は保証されていない，キャッシュ有効
利用は困難• より”memory-bound”

99

GeoFEM Benchmark
ICCG in FEM for Solid Mechanics
SR11K/J2

SR16K/M1

T2K

FX10

京

Core #/Node

16

32

16

16

8

Peak Performance(GFLOPS)

147.2

980.5

147.2

236.5

128.0

STREAM Triad (GB/s)

101.0

264.2

20.0

64.7

43.3

B/F

0.686

0.269

0.136

0.274

0.338

GeoFEM (GFLOPS)

19.0

72.7

4.69

16.0

11.0

% to Peak

12.9

7.41

3.18

6.77

8.59

LLC/core (MB)

18.0

4.00

2.00

0.75

0.75

Sparse Linear Solver: Memory-Bound100

Copy, Compile and Run>$ cd /work/gt00/t00XXX/pFEM>$ cp /work/gt00/z30088/pFEM/F/stream.tar .
>$ tar xvf stream.tar>$ cd mpi/stream>$ mpiifort -align array64byte -O3 -axCORE-AVX512 stream.f –o stream>$ pjsub XXX.shMemory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores101

s01.sh: Use 1 core#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=1#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o s01.lst
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream Cores are randomly
selected
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
export I_MPI_PIN_PROCESSOR_LIST=0
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Core are specified

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores102

s16.sh: Use 16 cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=16#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o s16.lst
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream Cores are randomly
selected
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
export I_MPI_PIN_PROCESSOR_LIST=0-15
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Core are specified

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores103s32.sh: Use 32 cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=32#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o s32.lst
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream Cores are randomly
selected
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
export I_MPI_PIN_PROCESSOR_LIST=0-15,28-43
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Core are specified

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores104

s48.sh: Use 48 cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=48#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o s48.lst
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream Cores are randomly
selected
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
export I_MPI_PIN_PROCESSOR_LIST=0-23,28-51
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Core are specified

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores105

s56.sh: Use 56 cores#!/bin/sh#PJM -N "test"#PJM -L rscgrp=tutorial#PJM -L node=1#PJM --mpi proc=56#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o s56.lst
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream Cores are randomly
selected
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
export I_MPI_PIN_PROCESSOR_LIST=0-55
mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./stream
mpiexec.hydra -n ${PJM_MPI_PROC} ./stream
Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Core are specified

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

106

Triad on a Single Node of OBCX
Peak is 281.57 GB/sec.
メモリ・キャッシュは「Random」の方がより効果的に利用されている（12-28コアの場合）
250○RadomRandom●Specified
Specified #0-#55#0-#55

GB/sec

200

150

100

50

0
1

10

Core #

100

107

Triad on a Single Node of OBCX
Peak is 281.57 GB/sec.
メモリ・キャッシュは「Random」の方がより効果的に利用されている（12-28コアの場合）
20

Speed-Up

15

10

5

Radoml○ Random
Specified#0-#55#0-#55● Specified
ー Ideal
Ideal

0
0

10

20

30

Core #

40

50

60

108
DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280

UPI

(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS

UPI
UPI

Soc. #0: 0th-27th cores

1

2(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

Socket #Intel® Xeon®
Platinum 8280

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

Triad on a Single
Node of OBCX
1コアに対する性能向上（Speed-Up）

Core #

RandomSpecified

1

1.000

1.000

2

1.998

1.963

4

3.912

6

5.792

8

7.615

7.177

12

11.0897.844

1612.2147.968

24

12.2787.945– 各ソケットに6枚のメモリを装着できる– 1ソケット6コア（6プロセス）まではほぼ比例28

15.9627.705 •

12-28 cores32

15.90115.86248

15.45215.49756

15.37015.358– Random：2ソケット利用– Specified/Compact：1ソケット• ソケット当たりメモリ
3.809
チャンネル数=6
5.682

109

Triad on a Single Node of OBCX
Peak is 281.57 GB/sec.
メモリ・キャッシュは「Random」の方がより効果的に利用されている（12-28コアの場合）
2-Sockets
Random250○RadomRandom●Specified
Specified #0-#55#0-#55

GB/sec

200

150

100

1-SocketSpecified

50

0
1

10

Core #

100

110

Triad on a Single Node of OBCX
Peak is 281.57 GB/sec.
●：メモリバンド幅は8-28コアでほぼ一定（飽和）32-56コアでは倍になる
250○RadomRandom●Specified
Specified #0-#55#0-#55

GB/sec

200

150

100

50

0
1

10

Core #

100

111

Exercises• Running the code
• Try various number of processes (1-56)• OpenMP-version and Single PE version areavailable– Fortran，C– Web-site of STREAM– http://www.cs.virginia.edu/stream/

900

112

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

800
700

Speed-Up

600

Up to 896 cores(16-nodes)

N=106

500
400
300
200
100
0
0

100

200

300

400

500

600

700

800

900

Core #
900

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

800
700

Speed-Up

600
500
400
300
200

N=107

100
0
0

100

200

300

400

500

Core #

600

700

800

900

• Performance at a Single
Core= 1.00
• ノード数が増すとN=106のケースは効率低下（台形積分と同様）• N=107のケースは徐々に
idealに近づき，672コア（12ノード）のケースではsuper linearになる

S2-ref

113

Strong-Scalingにおける「Super-Linear」
Speed-Up

ideal
super-linear

actualPE#• 問題規模を固定して，使用PE数を増加させて行った場合，通常は通信の影響のために，効率は理想値（m個のPEを使用した場合，理想的にはm倍の性能になる）よりも低くなるのが普通である。
• しかし，スカラープロセッサ（PC等）の場合，逆に理想値よりも，高い性能が出る場合がある。このような現象を「Super-Linear」と呼ぶ。
– ベクトル計算機では起こらない。

S2-ref

114

Super-Linearの生じる理由• キャッシュの影響• スカラープロセッサでは，全般に問題規模が小さいほど性能が高い。
– キャッシュの有効利用

CPU
Register
FAST

Cache
(MB)

3.00

Node # is smaller
Larger problem size
per MPI process

GFLOPS

2.50
2.00

SLOW

1.50
1.00
0.50

Main Memory(GB)

Node # is larger
Smaller problem size
per MPI process

0.00
1.0E+04

1.0E+05

1.0E+06

DOF: Problem Size

1.0E+07

S2-ref115

メモリーコピーも意外に時間かかる（1/2）
SENDbuf

neib#1neib#2neib#3neib#4BUFlength_e

BUFlength_e

BUFlength_e

BUFlength_e

export_index(0)+1

export_index(1)+1 export_index(2)+1

export_index(3)+1

do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= VAL(kk)enddo
Copied to sending buffers
enddo
do neib= 1, NEIBPETOTiS_e= export_index(neib-1) + 1
iE_e= export_index(neib )BUFlength_e= iE_e + 1 - iS_e
call MPI_ISEND&&(SENDbuf(iS_e), BUFlength_e, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_send(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_send, stat_send, ierr)

export_index(4)S2-ref116

メモリーコピーも意外に時間かかる（2/2）
do neib= 1, NEIBPETOTiS_i= import_index(neib-1) + 1
iE_i= import_index(neib )BUFlength_i= iE_i + 1 - iS_i
call MPI_IRECV&&(RECVbuf(iS_i), BUFlength_i, MPI_INTEGER, NEIBPE(neib), 0,&&
MPI_COMM_WORLD, request_recv(neib), ierr)enddo
call MPI_WAITALL (NEIBPETOT, request_recv, stat_recv, ierr)do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)VAL(kk)= RECVbuf(k)enddo
Copied from receiving buffer
enddo

RECVbuf

neib#1neib#2neib#3neib#4BUFlength_i

BUFlength_i

BUFlength_i

BUFlength_i

import_index(0)+1

import_index(1)+1 import_index(2)+1

import_index(3)+1

import_index(4)

900

117

800

28x2

24x2

16x2

Ideal

Up to 16-nodes

N=106

700

Random, O3+AVX512

Speed-Up

600
500• MPI Processes per Node400

– 16x2, 24x2, 28x2

300

• N=106

200
100
0
0

100

200

300

400

500

600

700

800

900

Core #– ノード数・コア数が増えると
差異少ない

900
800

28x2

24x2

16x2

Ideal

• N=107– コア数が同じであれば16x2が最良，メモリをより効率的に利用

700

Speed-Up

600
500
400
300
200

Memory
96 GB

N=107

100

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

0
0

100

200

300

400

500

Core #

600

700

800

900

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Socket #0: 0th-27th Cores

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

Socket #1: 28th-55th Cores

118

900

Random, O3+AVX512
Random, O3
Specified, O3+AVX512
Specified, O3
Ideal

800
700

Speed-Up

600

Up to 896 cores(16-nodes)

N=106

500
400
300
200
100
0
0
4500

200

300

400

500

600

700

800

900

Core #
Random, -O3+AVX512
Random, -O3
Specified, -O3+AVX512
Specified, -O3
Ideal

4000
3500
3000

Speed-Up100• 56コア（1ノード）の性能を56.0とする• 複数ノードの場合には，1ノード（56コア）での値を基準とするのがreasonable– L2 cache: 1MB/core– L3: 38.5MB/socket (shared)2500• N=1062000– 通信オーバーヘッド大1500• N=107

1000

N=107

500

– Superlinear強い

0
0

100

200

300

400

500

Core #

600

700

800

900

Intro119

Category

Capacity

X-Way
Set Associative

L1$Data

32 KB/core

8-Way

64B

L1$Instruction

32 KB/core

8-Way

64B

L2

1.00 MB/core

16-Way64B

L338.5 MB/socket

11-Way64B

Memory
96 GB

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

2933 MHz×6ch
140.8 GB/sec

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #0: 0th-27th cores

UPI
UPI
UPI

Intel® Xeon®
Platinum 8280
(Cascade Lake, CLX)
2.7GHz, 28-Cores
2.419 TFLOPS
Soc. #1: 28th-55th cores

Ultra Path Interconnect
10.4 GT/sec ×3= 124.8 GB/sec

Cache
Line

DDR4
DDR4
DDR4
DDR4
DDR4
DDR4

Memory96 GB2933 MHz×6ch
140.8 GB/sec

120

5回の計測における変動
N=106

N=106
100

Random, -O3+AVX512
Random, -O3
Specified, -O3+AVX512
Specified, -O3

80

Max-Min Fluctuation (%)

Max-Min Fluctuation (%)

100

60

40

20

0

Random, -O3+AVX512
Random, -O3
Specified, -O3+AVX512
Specified, -O3

80

60

40

20

0

56

112

224

336

Core #

448

672

896

56

112

224

336

Core #

448

672

896

121

“-O3+AVX512” and “-O3”Best Case of 5 Measurements+: -O3 is better, -: -O3 is worse“-O3” is rather faster in 106 cases with many cores
up to 1 node,56cores

up to 16 nodes,
896 cores

10.0

10.0

Random, 10^6Specified, 10^6
Random, 10^7
Specified, 10^7

6.0
4.0
2.0
0.0
1

2

4

8

16

24

(2.0)
(4.0)
(6.0)
(8.0)
(10.0)

Random, 10^6Specified, 10^6
Random, 10^7
Specified, 10^7

8.0

28

32

48

56

Performance Differnce (%)

Performance Differnce (%)

8.0

6.0
4.0
2.0
0.0
56

112

224

336

(2.0)
(4.0)
(6.0)
(8.0)

Core #(10.0)

Core #

448

672

896

S2-ref122

並列有限要素法：まとめ• 「局所分散データ構造の適切な設計」に尽きる• 問題点– 並列メッシュ生成，並列可視化– 悪条件問題における並列前処理手法– 大規模I/O

123

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

124

もし各プロセスにおいて，各隣接プロ
セスに対応した外点の番号が連続して
いたら
84

81

85

82

83

86

88

87

96

57

58

59

60

61

62

63

64

73

95

49

50

51

52

53

54

55

56

74

94

41

42

43

44

45

46

47

48

80

93

33

34

35

36

37

38

39

40

79

92

25

26

27

28

29

30

31

32

78

91

17

18

19

20

21

22

23

24

77

90

9

10

11

12

13

14

15

16

76

89

1

2

3

4

5

6

7

8

75

65

66

67

68

69

70

71

72

125

[A]{p}= {q} (Original), 1d.f
allocate (stat_send(MPI_STATUS_SIZE,NEIBPETOT))
allocate (stat_recv(MPI_STATUS_SIZE,NEIBPETOT))
allocate (request_send(NEIBPETOT)); allocate (request_recv(NEIBPETOT))
do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= W(kk,P)
enddo
enddo
do neib= 1, NEIBPETOT
is = export_index(neib-1) + 1
len_s= export_index(neib) - export_index(neib-1)
call MPI_Isend (SENDbuf(is), len_s, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD,&
request_send(neib), ierr)
enddo

&
&

do neib= 1, NEIBPETOT
ir = import_index(neib-1) + 1
len_r= import_index(neib) - import_index(neib-1)
call MPI_Irecv (RECVbuf(ir), len_r, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD,&
request_recv(neib), ierr)enddo
call MPI_Waitall (NEIBPETOT, request_recv, stat_recv, ierr)&&do neib= 1, NEIBPETOT
do k= import_index(neib-1)+1, import_index(neib)
kk= import_item(k)W(kk,P)= RECVbuf(k)
enddo
enddo
call MPI_Waitall (NEIBPETOT, request_send, stat_send, ierr)[A]{p}= {q} (Mod.): 受信でコピー不要
1d2.f
allocate (stat_send(MPI_STATUS_SIZE,2*NEIBPETOT))
allocate (request_send(2*NEIBPETOT))
do neib= 1, NEIBPETOT
do k= export_index(neib-1)+1, export_index(neib)
kk= export_item(k)SENDbuf(k)= W(kk,P)
enddo
enddo
do neib= 1, NEIBPETOT
is= export_index(neib-1) + 1
len_s= export_index(neib) - export_index(neib-1)
call MPI_Isend (SENDbuf(is), len_s, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD,&
request_send(neib), ierr)enddo
do neib= 1, NEIBPETOT
ir= import_index(neib-1) + 1
len_r= import_index(neib) - import_index(neib-1)
call MPI_Irecv (W(ir+N,P), len_r, MPI_DOUBLE_PRECISION,&
NEIBPE(neib), 0, MPI_COMM_WORLD,&
request_send(neib+NEIBPETOT), ierr)enddo
call MPI_Waitall (2*NEIBPETOT, request_send, stat_send, ierr)

&
&

&

&

126

127

“Original” and “Modified”
Best Case of 5 Measurements+: “Modified” is better, -: “Modified” is worse“Modified” is faster in 106 cases with many cores
But differences are small
up to 1 node,56cores

up to 16 nodes,
896 cores
10.0

10.0

6.0

Specified, 10^6

8.0

Specified, 10^6Specified, 10^7

6.0

Specified, 10^7

4.0
2.0
0.0
1

2

4

8

16

24

(2.0)
(4.0)
(6.0)

32

48

56

4.0
2.0
0.0
56

112

224

336

(2.0)
(4.0)
(6.0)
(8.0)

(8.0)
(10.0)

28

Performance Differnce (%)

Performance Differnce (%)

8.0

Core #(10.0)

Core #

448

672

896

S2-ref128

更なる検討• NUMA controlの有無による比較

mpiexec.hydra -n ${PJM_MPI_PROC} numactl –l ./a.out
mpiexec.hydra -n ${PJM_MPI_PROC} ./a.out• MPIプロセス数が増えても反復回数が変わらない理由