並列有限要素法による一次元定常熱伝導解析プログラム
C言語編
中島 研吾
東京大学情報基盤センター

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 /lustre/gt00/t00XXX/pFEM>$ cp /lustre/gt00/z30088/class_eps/F/s2r-f.tar .
>$ tar xvf s2r-f.tar
Cユーザー>$ cd /lustre/gt00/t00XXX/pFEM>$ cp /lustre/gt00/z30088/class_eps/C/s2r-c.tar .
>$ tar xvf s2r-c.tar
ディレクトリ確認・コンパイル>$ cd mpi/S2-ref>$ mpiifort -O3 -xCORE-AVX2 -align array32byte 1d.f>$ mpicc -O3 -xCORE-AVX2 -align 1d.c
このディレクトリを本講義では <$O-S2r> と呼ぶ。
<$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

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

#0

0

1
0

2

3
2

1

#1

4
3

4

0
3

1
0

2
1

3
2

#2

5
4

4

0
3

1
0

2
1

3
2

S2-ref

15

1D FEM: Numbering of Local ID
0

1
0

N

0
N-1

N

0

…
1

1
0

N-1
N-2

1

1

0
N-1

…

N-1

N-1

#PETot-1:
N+1 nodes
N elements

N-1

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

N-2

…
1

N

N-2

#0:
N+1 nodes
N elements

N+1NS2-ref

プログラム：1d.c（1/11）
諸変数#include <stdio.h>#include <stdlib.h>#include <math.h>#include <assert.h>#include <mpi.h>
int main(int argc, char **argv){

MPIを使用するときの「おまじない」

int NE, N, NP, NPLU, IterMax, NEg, Ng, errno;
double dX, Resid, Eps, Area, QV, COND, QN;
double X1, X2, DL, Ck; double *PHI, *Rhs, *X, *Diag, *AMat;
double *R, *Z, *Q, *P, *DD;
int *Index, *Item, *Icelnod;
double Kmat[2][2], Emat[2][2];int i, j, in1, in2, k, icel, k1, k2, jS;
int iter, nr, neib;
FILE *fp;
double BNorm2, Rho, Rho1=0.0, C1, Alpha, Beta, DNorm2;
int PETot, MyRank, kk, is, ir, len_s, len_r, tag;
int NeibPETot, BufLength, NeibPE[2];int import_index[3], import_item[2];int export_index[3], export_item[2];double SendBuf[2], RecvBuf[2];double BNorm20, Rho0, C10, DNorm20;
double StartTime, EndTime;
int ierr = 1;
MPI_Status *StatSend,*StatRecv;
MPI_Request *RequestSend, *RequestRecv;

16

S2-ref

プログラム：1d.c（2/11）
制御データ読み込み/*// +-------+// | INIT. |// +-------+//=== *//*//-- CONTROL data*/ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &PETot);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &MyRank);

MPI初期化：必須
全プロセス数：PETot
自分のランク番号（0～PETot-1）：MyRankif(MyRank == 0){
fp = fopen("input.dat", "r");assert(fp != NULL);
fscanf(fp, "%d", &NEg);
fscanf(fp, "%lf %lf %lf %lf", &dX, &QV, &Area, &COND);
fscanf(fp, "%d", &IterMax);
fscanf(fp, "%lf", &Eps);
fclose(fp);}ierr = MPI_Bcast(&NEg, 1, MPI_INT, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&IterMax, 1, MPI_INT, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&dX, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&QV, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Area, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&COND, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Eps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

17

S2-ref

プログラム：1d.c（2/11）
制御データ読み込み/*// +-------+// | INIT. |// +-------+//=== *//*//-- CONTROL data*/ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &PETot);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &MyRank);

MPI初期化：必須
全プロセス数：PETot
自分のランク番号（0～PETot-1）：MyRankif(MyRank == 0){
MyRank=0のとき制御データを読み込む
fp = fopen("input.dat", "r");assert(fp != NULL);
fscanf(fp, "%d", &NEg);
NEg：「全」要素数
fscanf(fp, "%lf %lf %lf %lf", &dX, &QV, &Area, &COND);
fscanf(fp, "%d", &IterMax);
fscanf(fp, "%lf", &Eps);
fclose(fp);}ierr = MPI_Bcast(&NEg, 1, MPI_INT, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&IterMax, 1, MPI_INT, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&dX, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&QV, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Area, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&COND, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Eps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

18

S2-ref

プログラム：1d.c（2/11）
制御データ読み込み/*// +-------+// | INIT. |// +-------+//=== *//*//-- CONTROL data*/ierr = MPI_Init(&argc, &argv);
ierr = MPI_Comm_size(MPI_COMM_WORLD, &PETot);
ierr = MPI_Comm_rank(MPI_COMM_WORLD, &MyRank);

MPI初期化：必須
全プロセス数：PETot
自分のランク番号（0～PETot-1）：MyRankif(MyRank == 0){
MyRank=0のとき制御データを読み込む
fp = fopen("input.dat", "r");assert(fp != NULL);
fscanf(fp, "%d", &NEg);
Neg：「全」要素数
fscanf(fp, "%lf %lf %lf %lf", &dX, &QV, &Area, &COND);
fscanf(fp, "%d", &IterMax);
fscanf(fp, "%lf", &Eps);
fclose(fp);}ierr = MPI_Bcast(&NEg, 1, MPI_INT, 0, MPI_COMM_WORLD); 0番プロセスから各プロセスにデータ送信
ierr = MPI_Bcast(&IterMax, 1, MPI_INT, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&dX, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&QV, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Area, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&COND, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
ierr = MPI_Bcast(&Eps, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

19

S2-ref20

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」にメッセージを送信。
• MPI_Bcast (buffer,count,datatype,root,comm)– buffer 任意– count

整数

I/O
I– datatype 整数 I

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

FORTRANMPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION, MPI_CHARACTER etc.

C

MPI_INT, MPI_FLOAT, MPI_DOUBLE, MPI_CHAR etc.

– root– comm

整数
整数

II

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

S2-ref

プログラム：1d.c（3/11）
局所分散メッシュデータ/*//-- LOCAL MESH size*/
Ng= NEg + 1;
N = Ng / PETot;nr = Ng - N*PETot;
if(MyRank < nr) N++;

Ng：総節点数
N ：局所節点数
NgがPETotで割り切れない場合

NE= N - 1 + 2;
NP= N + 2;
if(MyRank == 0) NE= N - 1 + 1;
if(MyRank == 0) NP= N + 1;
if(MyRank == PETot-1) NE= N - 1 + 1;
if(MyRank == PETot-1) NP= N + 1;
if(PETot==1){NE=N-1;}
if(PETot==1){NP=N ;}/*/-- Arrays*/
PHI = calloc(NP, sizeof(double));Diag = calloc(NP, sizeof(double));AMat = calloc(2*NP-2, sizeof(double));Rhs = calloc(NP, sizeof(double));Index= calloc(NP+1, sizeof(int));
Item = calloc(2*NP-2, sizeof(int));
Icelnod= calloc(2*NE, sizeof(int));21

S2-ref

22

プログラム：1d.c（3/11）
局所分散メッシュデータ，各要素→一様/*//-- LOCAL MESH size*/
Ng= NEg + 1;
N = Ng / PETot;nr = Ng - N*PETot;
if(MyRank < nr) N++;

Ng：総節点数
N ：局所節点数（内点）
NgがPETotで割り切れない場合

NE= N - 1 + 2;
局所要素数
NP= N + 2;
内点＋外点（局所総節点数）
if(MyRank == 0) NE= N - 1 + 1;
if(MyRank == 0) NP= N + 1;
if(MyRank == PETot-1) NE= N - 1 + 1;
if(MyRank == PETot-1) NP= N + 1;N0if(PETot==1){NE=N-1;}
if(PETot==1){NP=N ;}/*/-- Arrays*/
PHI = calloc(NP, sizeof(double));Diag = calloc(NP, sizeof(double));AMat = calloc(2*NP-2, sizeof(double));Rhs = calloc(NP, sizeof(double));Index= calloc(NP+1, sizeof(int));
Item = calloc(2*NP-2, sizeof(int));
Icelnod= calloc(2*NE, sizeof(int));

N-1

0

1

…
1

N-1
N-2

N+1
N

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

S2-ref

プログラム：1d.c（3/11）
局所分散メッシュデータ，各要素→一様/*//-- LOCAL MESH size*/
Ng= NEg + 1;
N = Ng / PETot;nr = Ng - N*PETot;
if(MyRank < nr) N++;

Ng：総節点数
N ：局所節点数（内点）
NgがPETotで割り切れない場合

NE= N - 1 + 2;
NP= N + 2;
if(MyRank == 0) NE= N - 1 + 1;
if(MyRank == 0) NP= N + 1;
if(MyRank == PETot-1) NE= N - 1 + 1;
if(MyRank == PETot-1) NP= N + 1;
if(PETot==1){NE=N-1;}
if(PETot==1){NP=N ;}/*/-- Arrays*/
PHI = calloc(NP, sizeof(double));Diag = calloc(NP, sizeof(double));AMat = calloc(2*NP-2, sizeof(double));Rhs = calloc(NP, sizeof(double));Index= calloc(NP+1, sizeof(int));
Item = calloc(2*NP-2, sizeof(int));
Icelnod= calloc(2*NE, sizeof(int));

0

1
0

…
1

N-1
N-2

N
N-1

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

23

S2-ref

プログラム：1d.c（3/11）
局所分散メッシュデータ，各要素→一様/*//-- LOCAL MESH size*/
Ng= NEg + 1;
N = Ng / PETot;nr = Ng - N*PETot;
if(MyRank < nr) N++;

Ng：総節点数
N ：局所節点数（内点）
NgがPETotで割り切れない場合

NE= N - 1 + 2;
NP= N + 2;
if(MyRank == 0) NE= N - 1 + 1;
if(MyRank == 0) NP= N + 1;
if(MyRank == PETot-1) NE= N - 1 + 1;
if(MyRank == PETot-1) NP= N + 1;
if(PETot==1){NE=N-1;}
if(PETot==1){NP=N ;}

N

0
N-1

/*
/-- Arrays*/
PHI = calloc(NP, sizeof(double));Diag = calloc(NP, sizeof(double));AMat = calloc(2*NP-2, sizeof(double));Rhs = calloc(NP, sizeof(double));Index= calloc(NP+1, sizeof(int));
Item = calloc(2*NP-2, sizeof(int));
Icelnod= calloc(2*NE, sizeof(int));

1
0

…
1

N-1
N-2

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

24

S2-ref

プログラム：1d.c（3/11）
局所分散メッシュデータ/*//-- LOCAL MESH size*/
Ng= NEg + 1;
N = Ng / PETot;nr = Ng - N*PETot;
if(MyRank < nr) N++;

Ng：総節点数
N ：局所節点数（内点）
NgがPETotで割り切れない場合

NE= N - 1 + 2;
NP= N + 2;
if(MyRank == 0) NE= N - 1 + 1;
if(MyRank == 0) NP= N + 1;
if(MyRank == PETot-1) NE= N - 1 + 1;
if(MyRank == PETot-1) NP= N + 1;
if(PETot==1){NE=N-1;}
if(PETot==1){NP=N ;}/*/-- Arrays*/
PHI = calloc(NP, sizeof(double));Diag = calloc(NP, sizeof(double));AMat = calloc(2*NP-2, sizeof(double));Rhs = calloc(NP, sizeof(double));Index= calloc(NP+1, sizeof(int));
Item = calloc(2*NP-2, sizeof(int));
Icelnod= calloc(2*NE, sizeof(int));

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

25

S2-ref

プログラム：1d.c（4/11）
配列初期化，要素～節点
for(i=0;i<NP;i++)U[i] = 0.0;
for(i=0;i<NP;i++) Diag[i] = 0.0;
for(i=0;i<NP;i++) Rhs[i] = 0.0;
for(k=0;k<2*NP-2;k++) AMat[k] = 0.0;
for(i=0;i<3;i++) import_index[i]= 0;
for(i=0;i<3;i++) export_index[i]= 0;
for(i=0;i<2;i++) import_item[i]= 0;
for(i=0;i<2;i++) export_item[i]= 0;
for(icel=0;icel<NE;icel++){Icelnod[2*icel ]= icel;
Icelnod[2*icel+1]= icel+1;}
if(PETot>1){
if(MyRank==0){
icel= NE-1;
Icelnod[2*icel ]= N-1;
Icelnod[2*icel+1]= N;}else if(MyRank==PETot-1){
icel= NE-1;
Icelnod[2*icel ]= N;
Icelnod[2*icel+1]= 0;}else{
icel= NE-2;
Icelnod[2*icel ]= N;
Icelnod[2*icel+1]= 0;
icel= NE-1;
Icelnod[2*icel ]= N-1;
Icelnod[2*icel+1]= N+1;
}
}

icel
Icelnod[2*icel]=icel

Icelnod[2*icel+1]=icel+126

S2-ref

プログラム：1d.c（4/11）
配列初期化，要素～節点
for(i=0;i<NP;i++)U[i] = 0.0;
for(i=0;i<NP;i++) Diag[i] = 0.0;
for(i=0;i<NP;i++) Rhs[i] = 0.0;
for(k=0;k<2*NP-2;k++) AMat[k] = 0.0;
for(i=0;i<3;i++) import_index[i]= 0;
for(i=0;i<3;i++) export_index[i]= 0;
for(i=0;i<2;i++) import_item[i]= 0;
for(i=0;i<2;i++) export_item[i]= 0;
for(icel=0;icel<NE;icel++){Icelnod[2*icel ]= icel;
Icelnod[2*icel+1]= icel+1;}
if(PETot>1){
if(MyRank==0){
icel= NE-1;
Icelnod[2*icel ]= N-1;
Icelnod[2*icel+1]= N;}else if(MyRank==PETot-1){
icel= NE-1;
Icelnod[2*icel ]= N;
Icelnod[2*icel+1]= 0;}else{
icel= NE-2;
Icelnod[2*icel ]= N;
Icelnod[2*icel+1]= 0;
icel= NE-1;
Icelnod[2*icel ]= N-1;
Icelnod[2*icel+1]= N+1;}}

「0-1」の要素を「0」とする

0

1
0

N

0
N-1

N

1
0

0
N-1

…

…
1

1
0

N-1
N-2

1

N-1

N-1
N-2

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

N-2

…
1

N
N-1

N+1
N

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

27

S2-ref

プログラム：1d.c（5/11）
Index定義
Kmat[0][0]= +1.0;
Kmat[0][1]= -1.0;
Kmat[1][0]= -1.0;
Kmat[1][1]= +1.0;/*// +--------------+// | CONNECTIVITY |// +--------------+*/
for(i=0;i<N+1;i++)Index[i] = 2;
for(i=N+1;i<NP+1;i++) Index[i] = 1;
Index[0] = 0;
if(MyRank == 0)
Index[1] = 1;
if(MyRank == PETot-1) Index[N] = 1;

0

1
0

…

N-1
N-2

1

N
N-1

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

for(i=0;i<NP;i++){Index[i+1]= Index[i+1] + Index[i];}NPLU= Index[NP];

N

0
N-1

N

1
0

0
N-1

…
1

1
0

N-1

…
1

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

N-2

N-1
N-2

N+1
N

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

28

S2-ref

プログラム：1d.c（6/11）
Item定義for(i=0;i<N;i++){jS = Index[i];if((MyRank==0)&&(i==0)){
Item[jS] = i+1;
}else if((MyRank==PETot-1)&&(i==N-1)){
Item[jS] = i-1;
}else{
Item[jS] = i-1;
Item[jS+1] = i+1;
if(i==0) { Item[jS] = N;}if(i==N-1){ Item[jS+1]= N+1;}
if((MyRank==0)&&(i==N-1)){Item[jS+1]= N;}
}
0
}
i =N;
jS= Index[i];if (MyRank==0) {Item[jS]= N-1;
}else {
Item[jS]= 0;
}

N

1
0

0
N-1

i =N+1;
jS= Index[i];if ((MyRank!=0)&&(MyRank!=PETot-1)) {Item[jS]= N-1;
N
0
}
N-1

…

1
0

…
1

1
0

N-1
N-2

1

N-1

N-1
N-2

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

N-2

…
1

N
N-1

N+1
N

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

29

S2-ref

30

プログラム：1d.c（7/11）
通信テーブル定義/*//-- COMMUNICATION*/
NeibPETot = 2;
if(MyRank == 0)
NeibPETot = 1;
if(MyRank == PETot-1) NeibPETot = 1;
if(PETot == 1)
NeibPETot = 0;
NeibPE[0] = MyRank - 1;
NeibPE[1] = MyRank + 1;0

if(MyRank == 0)
NeibPE[0] = MyRank + 1;
if(MyRank == PETot-1) NeibPE[0] = MyRank - 1;
import_index[1]=1;
import_index[2]=2;
import_item[0]= N;
import_item[1]= N+1;

N

if(MyRank == 0) import_item[0]=N;
if(MyRank == 0) export_item[0]=N-1;0
N-1

export_index[1]=1;
export_index[2]=2;
export_item[0]= 0;
export_item[1]= N-1;

1
0

N

1
0

0
N-1

…

…
1

1
0

BufLength = 1;
StatSend = malloc(sizeof(MPI_Status) * NeibPETot);
StatRecv = malloc(sizeof(MPI_Status) * NeibPETot);
RequestSend = malloc(sizeof(MPI_Request) * NeibPETot);
RequestRecv = malloc(sizeof(MPI_Request) * NeibPETot);

N-1
N-2

1

N-1

N-1
N-2

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

N-2

…
1

N
N-1

N+1
N

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

S2-ref31

MPI_Isend•

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

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

sendbuf 任意
count
整数
datatype 整数
dest
整数
tag
整数

I
I
I
I
I

整数
整数

I
O

– comm– request

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

S2-ref32

MPI_Irecv•

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

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

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

I
I
I
I
I

整数
整数

I
O

– comm– request

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

S2-ref

33

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

整数O

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

S2-ref34

一般化された通信テーブル：送信• 送信相手– NeibPETot，NeibPE[NeibPETot]• それぞれの送信相手に送るメッセージサイズ– export_index[NeibPETot+1]• 「境界点」番号– export_item[export_index[NeibPETot+1]]• それぞれの送信相手に送るメッセージ– SendBuf[export_index[NeibPETot+1]]

S2-ref35

送信（MPI_Isend/Irecv/Waitall）
SendBuf

neib#0neib#1neib#2neib#3

BUFlength_e

BUFlength_e

BUFlength_e

BUFlength_e

export_index[0]export_index[1]export_index[2]export_index[3]export_index[4]export_index[neib]～export_index[neib+1]-1番目のexport_itemがneib番目の隣接領域に送信されるfor (neib=0; neib<NeibPETot;neib++){for (k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= VAL[kk];
送信バッファへの代入}}
for (neib=0; neib<NeibPETot; neib++){tag= 0;iS_e= export_index[neib];iE_e= export_index[neib+1];
BUFlength_e= iE_e - iS_e}

ierr= MPI_Isend(&SendBuf[iS_e], BUFlength_e, MPI_DOUBLE, NeibPE[neib], 0,
MPI_COMM_WORLD, &ReqSend[neib])

MPI_Waitall(NeibPETot, ReqSend, StatSend);

S2-ref36

送信：一次元問題
4

0

1

2

3

5

• 受信相手SendBuf[0]=VAL[0]– NeibPETot，NeibPE[NeibPETot]

SendBuf[1]=VAL[3]• NeibPETot=2, NeibPE[0]= my_rank-1, NeibPE[1]= my_rank+1• それぞれの送信相手に送るメッセージサイズ– export_index[NeibPETot+1]• export_index[0]=0, export_index[1]= 1, export_index[2]= 2• 「境界点」番号– export_item[export_index[NeibPETot+1]]• export_item[0]= 0, export_item[1]= N-1• それぞれの送信相手に送るメッセージ– SendBuf[export_index[NeibPETot+1]]• SendBuf[0]= VAL[0], SendBuf[1]= VAL[N-1]S2-ref37

一般化された通信テーブル：受信• 受信相手– NeibPETot，NeibPE[NeibPETot]• それぞれの受信相手から受け取るメッセージサイズ– export_index[NeibPETot+1]• 「外点」番号– export_item[export_index[NeibPETot+1]]• それぞれの受信相手から受け取るメッセージ– SendBuf[export_index[NeibPETot+1]]

S2-ref38

受信（MPI_Isend/Irecv/Waitall）
for (neib=0; neib<NeibPETot; neib++){tag= 0;iS_i= import_index[neib];iE_i= import_index[neib+1];
BUFlength_i= iE_i - iS_i}

ierr= MPI_Irecv(&RecvBuf[iS_i], BUFlength_i, MPI_DOUBLE, NeibPE[neib], 0,
MPI_COMM_WORLD, &ReqRecv[neib])

MPI_Waitall(NeibPETot, ReqRecv, StatRecv);
for (neib=0; neib<NeibPETot;neib++){for (k=import_index[neib];k<import_index[neib+1];k++){kk= import_item[k];VAL[kk]= RecvBuf[k];受信バッファからの代入}}
import_index[neib]～import_index[neib+1]-1番目のimport_itemがneib番目の隣接領域から受信される

RecvBuf

neib#0neib#1neib#2neib#3

BUFlength_i

BUFlength_i

BUFlength_i

BUFlength_i

import_index[0]import_index[1]import_index[2]import_index[3]import_index[4]S2-ref39

受信:一次元問題
4

0

1

2

3

5

• 受信相手
VAL[4]=RecvBuf[0]– NeibPETot，NeibPE[NeibPETot]

VAL[5]=RecvBuf[1]• NeibPETot=2, NeibPE[0]= my_rank-1, NeibPE[1]= my_rank+1• それぞれの受信相手から受け取るメッセージサイズ– import_index [NeibPETot+1]• import_index[0]=0, import_index[1]= 1, import_index[2]= 2• 「外点」番号– import_item [import_index[NeibPETot+1]]• import_item[0]= N, import_item[1]= N+1• それぞれの受信相手から受け取るメッセージ– import_item [import_index[NeibPETot+1]]• VAL[N]=RecvBuf[0], VAL[N+1]=RecvBuf[1]S2-ref40

一般化された通信テーブル：Fortran5

1

2

SENDbuf(1)=BUF(1)5

1

BUF(5)=RECVbuf(1)

2

3

4

6

SENDbuf(2)=BUF(4)

3

4

6

BUF(6)=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-ref41

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

1

SENDbuf[0]=BUF[0]4

0

BUF[4]=RECVbuf[0]

1

2

3

5

SENDbuf[1]=BUF[3]

2

3

5

BUF[5]=RECVbuf[1]NEIBPETOT= 2
NEIBPE[0]= my_rank - 1
NEIBPE[1]= my_rank + 1
import_index[1]= 0
import_index[2]= 1
import_item [0]= N
import_item [1]= N+1
export_index[1]= 0
export_index[2]= 1
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-ref

42

プログラム：1d.c（8/11）
全体マトリクス生成：1CPUのときと全く同じ：各要素→一様/*// +-----------------+// | MATRIX assemble |// +-----------------+*/
for(icel=0;icel<NE;icel++){in1= Icelnod[2*icel];
in2= Icelnod[2*icel+1];
DL = dX;#0

Ck= Area*COND/DL;
Emat[0][0]= Ck*Kmat[0][0];Emat[0][1]= Ck*Kmat[0][1];Emat[1][0]= Ck*Kmat[1][0];Emat[1][1]= Ck*Kmat[1][1];Diag[in1]= Diag[in1] + Emat[0][0];Diag[in2]= Diag[in2] + Emat[1][1];if ((MyRank==0)&&(icel==0)){
k1=Index[in1];}else {k1=Index[in1]+1;}
k2=Index[in2];
AMat[k1]= AMat[k1] + Emat[0][1];AMat[k2]= AMat[k2] + Emat[1][0];}QN= 0.5*QV*Area*dX;
Rhs[in1]= Rhs[in1] + QN;
Rhs[in2]= Rhs[in2] + QN;

0

1
0

2

3
2

1

#1

4
3

4

0
3

1
0

2
1

3
2

#2

5
4

4

0
3

1
0

2
1

3
2

FEM3D

43

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

NP

N

NP

N

internal

external

FEM3D44

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

NP

N

NP

N

internal

external

pFEM3D-2

45

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

46

全ての要素の計算を実施する
外点を含むオーバーラップ領域の要素の計算も実施
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

FEM3D47

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

NP

N

NP

N

internal

external

FEM3D48

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

NP

N

NP

N

internal

external

S2-ref

49

プログラム：1d.c（9/11）
境界条件：1CPUのときとほとんど同じ/*// +---------------------+// | BOUNDARY conditions |// +---------------------+*/#0/* X=Xmin */if (MyRank==0){i=0;jS= Index[i];AMat[jS]= 0.0;
Diag[i ]= 1.0;
Rhs [i ]= 0.0;

}

0

1
0

for(k=0;k<NPLU;k++){if(Item[k]==0){AMat[k]=0.0;}
}

2

3
2

1

#1

4
3

4

0
3

1
0

2
1

3
2

#2

5
4

4

0
3

1
0

2
1

3
2

S2-ref

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

/*
// +---------------+// | CG iterations |// +---------------+//=== */
R = calloc(NP, sizeof(double));Z = calloc(NP, sizeof(double));P = calloc(NP, sizeof(double));Q = calloc(NP, sizeof(double));DD= calloc(NP, sizeof(double));for(i=0;i<N;i++){DD[i]= 1.0 / Diag[i];}/*//-- {r0}= {b} - [A]{xini} |*/for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= U[kk];}}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|end50

S2-ref51

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

S2-ref

52

前処理，DAXPY/*//-- {z}= [Minv]{r}*/for(i=0;i<N;i++){Z[i] = DD[i] * R[i];}/*//-- {x}= {x} + ALPHA*{p}// {r}= {r} - ALPHA*{q}*/for(i=0;i<N;i++){U[i] += Alpha * P[i];R[i] -= Alpha * Q[i];}S2-ref53

行列ベクトル積（1/2）
通信テーブル使用，{p}の最新値を計算前に取得/*//-- {q}= [A]{p}*/for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= P[kk];}}for(neib=0;neib<NeibPETot;neib++){is= export_index[neib];
len_s= export_index[neib+1] - export_index[neib];
MPI_Isend(&SendBuf[is], len_s, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib]);}for(neib=0;neib<NeibPETot;neib++){ir= import_index[neib];
len_r= import_index[neib+1] - import_index[neib];
MPI_Irecv(&RecvBuf[ir], len_r, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestRecv[neib]);}
MPI_Waitall(NeibPETot, RequestRecv, StatRecv);
for(neib=0;neib<NeibPETot;neib++){for(k=import_index[neib];k<import_index[neib+1];k++){kk= import_item[k];P[kk]=RecvBuf[k];}}S2-ref54

行列ベクトル積（2/2）{q}= [A]{p}MPI_Waitall(NeibPETot, RequestSend, StatSend);
for(i=0;i<N;i++){Q[i] = Diag[i] * P[i];for(j=Index[i];j<Index[i+1];j++){Q[i] += AMat[j]*P[Item[j]];}}S2-ref55

内積各プロセスで計算した値を，MPI_Allreduceで合計/*//-- RHO= {r}{z}*/Rho0= 0.0;
for(i=0;i<N;i++){
Rho0 += R[i] * Z[i];}ierr = MPI_Allreduce(&Rho0, &Rho, 1, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD);S2-ref

56
P#0 A0 B0 C0 D0

MPI_AｌｌreduceP#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• 総和，最大値等を計算して，全プロセスに配信• MPI_Allreduce(sendbuf,recvbuf,count,datatype,op, comm)– sendbuf– recvbuf

任意
任意

I
O

– count– datatype– op– comm

整数
整数
整数
整数

I
I
I
I

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

S2-ref57

CG法（1/5）/*//-- {r0}= {b} - [A]{xini} |*/for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= PHI[kk];}}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)for(neib=0;neib<NeibPETot;neib++){is = export_index[neib];
else
len_s= export_index[neib+1] - export_index[neib];
βi-1= ρi-1/ρi-2
MPI_Isend(&SendBuf[is], len_s, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib]);(i)= z(i-1) + β
(i-1)
p
p
i-1
}
endif
for(neib=0;neib<NeibPETot;neib++){q(i)= [A]p(i)ir = import_index[neib];
len_r= import_index[neib+1] - import_index[neib];αi = ρi-1/p(i)q(i)MPI_Irecv(&RecvBuf[ir], len_r, MPI_DOUBLE, NeibPE[neib],x(i)= x(i-1) + αip(i)0, MPI_COMM_WORLD, &RequestRecv[neib]);}r(i)= r(i-1) - αiq(i)MPI_Waitall(NeibPETot, RequestRecv, StatRecv);
check convergence |r|end
for(neib=0;neib<NeibPETot;neib++){for(k=import_index[neib];k<import_index[neib+1];k++){kk= import_item[k];PHI[kk]=RecvBuf[k];}}MPI_Waitall(NeibPETot, RequestSend, StatSend);

S2-ref58

CG法（2/5）for(i=0;i<N;i++){R[i] = Diag[i]*PHI[i];for(j=Index[i];j<Index[i+1];j++){R[i] += AMat[j]*PHI[Item[j]];}}BNorm20 = 0.0;
for(i=0;i<N;i++){
BNorm20 += Rhs[i] * Rhs[i];R[i] = Rhs[i] - R[i];}ierr = MPI_Allreduce(&BNorm20, &BNorm2, 1, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD);
for(iter=1;iter<=IterMax;iter++){/*//-- {z}= [Minv]{r}*/for(i=0;i<N;i++){Z[i] = DD[i] * R[i];}/*//-- RHO= {r}{z}*/Rho0= 0.0;
for(i=0;i<N;i++){
Rho0 += R[i] * Z[i];}ierr = MPI_Allreduce(&Rho0, &Rho, 1, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD);

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

CG法（3/5）/*//-- {p} = {z} if
ITER=1//
BETA= RHO / RHO1 otherwise*/
if(iter == 1){
for(i=0;i<N;i++){P[i] = Z[i];}}else{
Beta = Rho / Rho1;
for(i=0;i<N;i++){P[i] = Z[i] + Beta*P[i];}}/*//-- {q}= [A]{p}*/for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= P[kk];}}for(neib=0;neib<NeibPETot;neib++){is = export_index[neib];
len_s= export_index[neib+1] - export_index[neib];
MPI_Isend(&SendBuf[is], len_s, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib]);}

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

CG法（4/5）
for(neib=0;neib<NeibPETot;neib++){ir = import_index[neib];
Compute r(0)= b-[A]x(0)len_r= import_index[neib+1] - import_index[neib];
MPI_Irecv(&RecvBuf[ir], len_r, MPI_DOUBLE, NeibPE[neib], for i= 1, 2, …
0, MPI_COMM_WORLD, &RequestRecv[neib]);
solve [M]z(i-1)= r(i-1)}
MPI_Waitall(NeibPETot, RequestRecv, StatRecv);
for(neib=0;neib<NeibPETot;neib++){for(k=import_index[neib];k<import_index[neib+1];k++){kk= import_item[k];P[kk]=RecvBuf[k];}}MPI_Waitall(NeibPETot, RequestSend, StatSend);
for(i=0;i<N;i++){Q[i] = Diag[i] * P[i];for(j=Index[i];j<Index[i+1];j++){Q[i] += AMat[j]*P[Item[j]];}}ρ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|/*//-- ALPHA= RHO / {p}{q}end*/
C10 = 0.0;
for(i=0;i<N;i++){C10 += P[i] * Q[i];}ierr = MPI_Allreduce(&C10, &C1, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
Alpha = Rho / C1;

S2-ref61

CG法（5/5）/*//-- {x}= {x} + ALPHA*{p}//{r}= {r} - ALPHA*{q}*/for(i=0;i<N;i++){
PHI[i] += Alpha * P[i];R[i] -= Alpha * Q[i];}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=1
DNorm20 = 0.0;p(1)= z(0)for(i=0;i<N;i++){
DNorm20 += R[i] * R[i];else}βi-1= ρi-1/ρi-2
ierr = MPI_Allreduce(&DNorm20, &DNorm2, 1, MPI_DOUBLE,(i)= z(i-1) + β
(i-1)
p
p
i-1
MPI_SUM, MPI_COMM_WORLD);
endif
Resid = sqrt(DNorm2/BNorm2);q(i)= [A]p(i)if (MyRank==0)
printf("%8d%s%16.6e¥n", iter, " iters, RESID=", Resid);
αi = ρi-1/p(i)q(i)x(i)= x(i-1) + αip(i)if(Resid <= Eps){
ierr = 0;r(i)= r(i-1) - αiq(i)break;}
check convergence |r|end
Rho1 = Rho;}S2-ref

プログラム：1d.c（11/11）
結果書き出し：各プロセスごとに実施/*//-- OUTPUT*/
printf("¥n%s¥n", "### TEMPERATURE");for(i=0;i<N;i++){printf("%3d%8d%16.6E¥n", MyRank, i+1, PHI[i]);}}ierr = MPI_Finalize();
return ierr;

62

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

63

S2-ref64

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 Cores65

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 Cores66

S2-ref67

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

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

S2-ref69

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

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

S2-ref71

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

Intro72

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

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

74

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75

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

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

76

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”

77

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78

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

79

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80s16.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81s32.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82

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83s56.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

84

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

Radom
Specified #0-#55

GB/sec

200

150

100

50

0
1

10

Core #

100

85

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

Speed-Up

15

10

5

RadomlSpecified #0-#55
Ideal

0
0

10

20

30

Core #

40

50

60

86
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

87

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

250

Radom
Specified #0-#55

GB/sec

200

150

100

1-Socket■ go1.sh, compact

50

0
1

10

Core #

100

88

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

Radom
Specified #0-#55

GB/sec

200

150

100

50

0
1

10

Core #

100

89

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

90

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

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

ideal
super-linear

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

S2-ref

92

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

メモリーコピーも意外に時間かかる（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-ref94

メモリーコピーも意外に時間かかる（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

95

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

96

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=106

500

– Superlinear強い

0
0

100

200

300

400

500

Core #

600

700

800

900

Intro97

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

98

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

99

“-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-ref100

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

101

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

102

もし各プロセスにおいて，各隣接プロ
セスに対応した外点の番号が連続していたら
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

103

[A]{p}= {q} (Original): 1d.cStatSend = malloc(sizeof(MPI_Status) * NeibPETot);
StatRecv = malloc(sizeof(MPI_Status) * NeibPETot);
RequestSend = malloc(sizeof(MPI_Request) * NeibPETot);
RequestRecv = malloc(sizeof(MPI_Request) * NeibPETot);
for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= P[kk];}}for(neib=0;neib<NeibPETot;neib++){is = export_index[neib];
len_s= export_index[neib+1] - export_index[neib];
MPI_Isend(&SendBuf[is], len_s, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib]);}for(neib=0;neib<NeibPETot;neib++){ir = import_index[neib];
len_r= import_index[neib+1] - import_index[neib];
MPI_Irecv(&RecvBuf[ir], len_r, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestRecv[neib]);}
MPI_Waitall(NeibPETot, RequestRecv, StatRecv);
for(neib=0;neib<NeibPETot;neib++){for(k=import_index[neib];k<import_index[neib+1];k++){kk= import_item[k];P[kk]=RecvBuf[k];}}MPI_Waitall(NeibPETot, RequestSend, StatSend);

104[A]{p}= {q} (Mod.): No Copy for RECV:
1d2.c
StatSend = malloc(sizeof(MPI_Status) * 2 * NeibPETot);
RequestSend = malloc(sizeof(MPI_Request) * 2 * NeibPETot);
for(neib=0;neib<NeibPETot;neib++){for(k=export_index[neib];k<export_index[neib+1];k++){kk= export_item[k];SendBuf[k]= P[kk];}}for(neib=0;neib<NeibPETot;neib++){is = export_index[neib];
len_s= export_index[neib+1] - export_index[neib];
MPI_Isend(&SendBuf[is], len_s, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib]);}for(neib=0;neib<NeibPETot;neib++){ir = import_index[neib];
len_r= import_index[neib+1] - import_index[neib];
MPI_Irecv(&P[ir+N], len_r, MPI_DOUBLE, NeibPE[neib],0, MPI_COMM_WORLD, &RequestSend[neib+NeibPETot]);}
MPI_Waitall(2*NeibPETot, RequestSend, StatSend);

105“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