OpenMPによるマルチコア・
メニィコア並列プログラミング入門
C言語編
Part-A2： OpenMP

中島研吾
東京大学情報基盤センター

2• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

3

Hybrid並列プログラミング• スレッド並列＋メッセージパッシング– OpenMP+ MPI– CUDA + MPI, OpenACC + MPI• OpenMPがMPIより簡単ということはない– データ依存性のない計算であれば，機械的にOpenMP指示文を入れれば良い– NUMAになるとより複雑：First Touch Data Placement

4Flat MPI vs. Hybrid
Flat-MPI：Each Core -> Independent

core
core

core
core

core
memorycore

core
memorymemory

core

core

core
core
core

core
core
core
core

memory

core
core
core
core

memorymemoryHybrid：Hierarchal Structure
core
core
core
core

OMP-15

共有メモリ型計算機
MEMORY

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

• SMP
– Symmetric Multi Processors– 複数のCPU（コア）で同じメモリ空間を共有するアーキテクチュアOMP-1

6OpenMPとは（1/2）
http://www.openmp.org• 共有メモリ型並列計算機用のDirectiveの統一規格– MPIやHPFに比べると遅く1996年頃から活動開始– 現在 Ver.4.X• 背景– CrayとSGIの合併（1996）– ASCI計画の開始（1995）• Accelerated Strategic Computing Initiative (ASCI) -> Advanced
Simulation and Computing Program (ASC)• ASCI：核実験のシミュレーションによる代替– 計算機開発，シミュレーションソフトウェア• SMPクラスタにフォーカス– ASCI Red (Intel)，Blue Pacific/White/Purple/BlueGene (IBM)，Blue Mountain（SGI）• SMPクラスタ向けの並列プログラミングの共通API（Application Program
Interface）が必要

OMP-1

7OpenMPとは（2/2）
http://www.openmp.org• 主な計算機ベンダーが集まって OpenMP ARBを結成し，1997年にはもう規格案ができていたそうである– SC98ではすでにOpenMPのチュートリアルがあったし，すでにSGI Origin2000でOpenMP-MPIハイブリッドのシミュレーショ
ンをやっている例もあった。

• OpenMPはFortan版とC/C++版の規格が全く別々に進められてきた。
– Ver.2.5で言語間の仕様を統一• Ver.4.0ではGPU，Intel-MIC等Co-Processor，Accelerator環境での動作も考慮– OpenACCに近づいているOMP-1

8OpenMPの概要• 基本的仕様– プログラムを並列に実行するための動作をユーザーが明示– OpenMP実行環境は，依存関係，衝突，デッドロック，競合条
件，結果としてプログラムが誤った実行につながるような問題に関するチェックは要求されていない。
– プログラムが正しく実行されるよう構成するのはユーザーの責任である。

• 実行モデル– fork-join型並列モデル• 当初はマスタスレッドと呼ばれる単一プログラムとして実行を開始し，「PARALLEL」，「END PARALLEL」ディレクティヴの対で並列構造を構成する。並列構造が現れるとマスタスレッドはスレッドのチームを生成し，そのチームのマスタとなる。

– いわゆる「入れ子構造」も可能であるが，ここでは扱わない

OMP-19

Fork-Join 型並列モデル
MasterMasterthreadthreadthreadthreadthreadthreadMaster

PARALLEL
fork

MasterMasterthreadthreadthreadthreadthreadthreadthreadthreadEND PARALLEL
join

PARALLEL
fork

END PARALLEL
join

OMP-110

スレッド数• 環境変数 OMP_NUM_THREADS– 値の変え方• bash（.bashrc）• csh(.cshrc)

export OMP_NUM_THREADS=8
setenv OMP_NUM_THREADS 8• たとえば，OMP_NUM_THREADS=4とすると，以下のように
i=1~100のループが4分割され，同時に実行される。
do i= 1, 25
do i= 26, 50
do i= 1,100
do i= 51, 75
do i= 76, 100

OMP-1

11

OpenMPに関連する情報• OpenMP Architecture Review Board (ARB)– http://www.openmp.org• 参考文献– Chandra, R. et al.「Parallel Programming in OpenMP」（Morgan Kaufmann）– Quinn, M.J. 「Parallel Programming in C with MPI and
OpenMP」（McGrawHill）– Mattson, T.G. et al. 「Patterns for Parallel Programming」（Addison Wesley）– 牛島「OpenMPによる並列プログラミングと数値計算法」（丸善）– Chapman, B. et al. 「Using OpenMP」（MIT Press）• 富士通他による翻訳：（OpenMP 3.0）必携 !
– http://www.openmp.org/mp-documents/OpenMP30spec-ja.pdf

OMP-1

12OpenMPの特徴• ディレクティヴ（指示行）の形で利用– 挿入直後のループが並列化される– コンパイラがサポートしていなければ，コメントとみなされる

13

OpenMP/Directives
Array Operations
Simple Substitution#pragma omp parallel for private (i)for (i=0; i<N; i++) {X[i] = 0.0;
W[0][i] = 0.0;
W[1][i] = 0.0;
W[2][i] = 0.0;}

DAXPY
#pragma omp parallel for private (i)for (i=0; i<N; i++) {Y[i] = Y[i} + alphat*X[i];}Dot Products
RHO = 0.0;
#pragma omp parallel for private (i)reduction (+:RHO)
for(i=0; i<N; i++) {RHO += W[R][i] * W[Z][i];}OMP-1

14OpenMP/Direceives
Matrix/Vector Products#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]-1];}
W[Q][i] = VAL;}OMP-1

15OpenMPの特徴• ディレクティヴ（指示行）の形で利用– 挿入直後のループが並列化される– コンパイラがサポートしていなければ，コメントとみなされる• 何も指定しなければ，何もしない– 「自動並列化」，「自動ベクトル化」とは異なる。
– 下手なことをするとおかしな結果になる：ベクトル化と同じ– データ分散等（Ordering）は利用者の責任• 共有メモリユニット内のプロセッサ数に応じて，「Thread」が立ち上がる– 「Thread」：MPIでいう「プロセス」に相当する。
– 普通は「Thread数＝共有メモリユニット内プロセッサ数，コア数」であるが最近のアーキテクチャではHyper Threading（HT）がサポートされているものが多い（1コアで2-4スレッド）OMP-116メモリ競合（Memory Contention）
MEMORY

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

• 複雑な処理をしている場合，複数のスレッドがメモリ上の同じアドレスにあるデータを同時に更新する可能性がある。
– 複数のCPUが配列の同じ成分を更新しようとする。
– メモリを複数のコアで共有しているためこのようなことが起こりうる。
– 場合によっては答えが変わるOMP-117メモリ競合（Memory Contention）（続き）
MEMORY

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

C
P
U

• 本演習で扱っている例は，そのようなことが生じないよう，
各スレッドが同時に同じ成分を更新するようなことはな
いようにする（実はそもそも無い） 。
– これはユーザーの責任でやること，である。

• データ依存性（Data Dependency）• コア数（スレッド数）が増えるほど，メモリへの負担が増えて，処理速度は低下する（メモリ飽和）。

OMP-1

18

OpenMPの特徴（続き）• 基本は「#pragma omp parallel for」
• 変数について，グローバル/sharedな変数と，Thread内
でローカルな「private」な変数に分けられる。
– デフォルトは「global/shared」– 内積を求める場合は「reduction」を使う
VAL= 0.0;#pragma omp parallel for private (i,ip)
reduction(+:VAL)
for(ip=0; ip<PEsmpTOT; ip++){for (i=INDEX[ip]; i<INDEX[ip+1]; i++){
VAL= VAL + W[R][i] * W[Z][i];}}W(:,:)，R，Z，PEsmpTOTなどはグローバル変数

OMP-1

19

FORTRANとC
use omp_lib...
!$omp parallel do shared(n,x,y) private(i)do i= 1, nx(i)= x(i) + y(i)enddo!$ omp end parallel do#include <omp.h>...
{#pragma omp parallel for default(none) shared(n,x,y) private(i)for (i=0; i<n; i++)x[i] += y[i];}OMP-120

本講義における方針• OpenMPは多様な機能を持っているが，それらの全てを逐一教えることはしない。
– 講演者も全てを把握，理解しているわけではない。

• 数値解析に必要な最低限の機能のみ学習する。
– 具体的には，講義で扱っているPCG法によるポアソン方程式
ソルバーを動かすために必要な機能のみについて学習する– それ以外の機能については，自習，質問のこと（全てに答えられるとは限らない）。

OMP-1

21

最初にやること• use omp_lib• #include <omp.h>

FORTRANC• 様々な環境変数，インタフェースの定義（OpenMP3.0
以降でサポート）OMP-1

22

OpenMPディレクィヴ（FORTRAN）
sentinel directive_name [clause[[,] clause]…]• 大文字小文字は区別されない。
• sentinel– 接頭辞– FORTRANでは「!$OMP」，「C$OMP」，「*$OMP」，但し自由
ソース形式では「!$OMP」のみ。
– 継続行にはFORTRANと同じルールが適用される。以下はいずれも「!$OMP PARALLEL DO SHARED(A,B,C)」!$OMP PARALLEL DO!$OMP+SHARED (A,B,C)!$OMP PARALLEL DO &!$OMP SHARED (A,B,C)

OMP-1

23

OpenMPディレクィヴ（C）#pragma omp directive_name [clause[[,] clause]…]• 継続行は「＼」• 小文字を使用（変数名以外）#pragma omp parallel for shared (a,b,c)

OMP-1

24

PARALLEL DO/for!$OMP PARALLEL DO[clause[[,] clause] … ](do_loop)!$OMP END PARALLEL DO#pragma omp parallel for [clause[[,] clause] … ](for_loop)• 多重スレッドによって実行される領域を定義し，DOルー
プの並列化を実施する。
• 並び項目（clause）：よく利用するもの– private（list）– shared（list）– default（private|shared|none）– reduction（{operation|intrinsic}: list）

OMP-1

25

REDUCTION
REDUCTION ({operator|intrinsic}: list)
reduction ({operator|intrinsic}: list)• 「MPI_REDUCE」のようなものと思えばよい• Operator– ＋，*，-，.AND.，.OR.，.EQV.，.NEQV.

• Intrinsic– MAX，MIN，IAND，IOR，IEQROMP-126

実例A1：簡単なループ#pragma omp parllel for
for(i=0; i<N; i++){B[i]= (A[i] + B[i]) * 0.50;}• ループの繰り返し変数（ここでは「i」）はデフォルトでprivateなので，明示的に宣言は不要。
• 「END PARALLEL DO」は省略可能。
– C言語ではそもそも存在しないOMP-127

実例A2：REDUCTION#pragma omp parallel default(private) reduction(+:A,B)
for(i=0; i<N; i++){err= work(Alocal, Blocal);
A= A + Alocal;
B= B + Blocal;}

28

OpenMP使用時に呼び出すことのできる
関数群
関数名

内容

int omp_get_num_threads (void)
int omp_get_thread_num (void)
double omp_get_wtime (void)

スレッド総数
自スレッドのID
MPI_Wtimeと同じ

void omp_set_num_threads (int
num_threads)
call omp_set_num_threads (num_threads)

スレッド数設定

29

OpenMP for Dot Products
VAL= 0.0;
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}30OpenMP for Dot Products
VAL= 0.0;
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}VAL= 0.0;#pragma omp parallel for private (i) reduction(+:VAL)
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}Thread #0 W[R][i]*W[Z][i]+Thread #1 W[R][i]*W[Z][i]

+

・
・
・

Thread #2 W[R][i]*W[Z][i]Thread #M W[R][i]*W[Z][i]

+

+

VAL

Directives are just inserted.

31

OpenMP for Dot Products
VAL= 0.0;
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}VAL= 0.0;#pragma omp parallel for private (i) reduction(+:VAL)
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}VAL= 0.0;#pragma omp parallel for private (i,ip)
reduction(+:VAL)
for(ip=0; ip<PEsmpTOT; ip++){for (i=INDEX[ip]; i<INDEX[ip+1]; i++){
VAL= VAL + W[R][i] * W[Z][i];}}OpenMPディレクティヴの挿入
これでも並列計算は可能

多重ループの導入
PEsmpTOT：スレッド数
あらかじめ「INDEX(:)」を用意しておくより確実に並列計算実施（別に効率がよくなるわけでは無い）

32

OpenMP for Dot Products
VAL= 0.0;
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}VAL= 0.0;#pragma omp parallel for private (i) reduction(+:VAL)
for(i=0; i<N; i++){
VAL= VAL + W[R][i] * W[Z][i];}VAL= 0.0;#pragma omp parallel for private (i,ip)
reduction(+:VAL)
for(ip=0; ip<PEsmpTOT; ip++){for (i=INDEX[ip]; i<INDEX[ip+1]; i++){
VAL= VAL + W[R][i] * W[Z][i];}}OpenMPディレクティヴの挿入
これでも並列計算は可能

多重ループの導入
PEsmpTOT：スレッド数
あらかじめ「INDEX(:)」を用意しておくより確実に並列計算実施
PEsmpTOT個のスレッドが立ち上がり，並列に実行

33

OpenMP for Dot Products
VAL= 0.0;#pragma omp parallel for private (i,ip)
reduction(+:VAL)
for(ip=0; ip<PEsmpTOT; ip++){for (i=INDEX[ip]; i<INDEX[ip+1]; i++){
VAL= VAL + W[R][i] * W[Z][i];}}多重ループの導入
PEsmpTOT：スレッド数
あらかじめ「INDEX[:]」を用意しておくより確実に並列計算実施
PEsmpTOT個のスレッドが立ち上がり，並列に実行

各要素が計算されるスレッドを指定できる
e.g.: N=100, PEsmpTOT=4
INDEX[0]= 0
INDEX[1]= 25
INDEX[2]= 50
INDEX[3]= 75
INDEX[4]= 100

34

Matrix-Vector Multiplyfor(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}

#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}

35

Matrix-Vector Multiplyfor(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}

#pragma omp parallel for private (i,ip.VAL,j)for (ip=0; i<PEsmpTOT; ip++){for(i=index[ip]; i<index[ip+1]; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;
}
}

36

omp parallel (do)• omp parallel-omp end
parallelはそのたびにスレッドを生成，消滅させる：fork-join• ループが連続するとオーバーヘッドになる。
• omp parallel + omp
do/omp for
MasterMasterthreadthreadthreadthreadthreadMaster

PARALLEL
fork

threadMasterMasterthreadthreadthreadthreadthreadthreadthreadthreadEND PARALLEL
join

PARALLEL
fork

END PARALLEL
join#pragma omp parallel ...
#pragma omp for {...
#pragma omp for {!$omp parallel ...
!$omp dodo i= 1, N...
!$omp dodo i= 1, N...
!$omp end parallelrequired37• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

38• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

39

ここでの目標• プログラムの並列化• OpenMP指示行を入れるだけでやってみる– “poi_gen.f/c”(poi_gen)，– “solver_PCG.f/c” (solve_PCG)

OMP-140

前処理付共役勾配法
Preconditioned Conjugate Gradient Method （PCG）
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

実際にやるべき計算は：{z} = [M ] {r}−1「近似逆行列」の計算が必要：[M ] ≈ [A] , [M ] ≈ [A]−1−1

究極の前処理：本当の逆行列[M ] = [A] , [M ] = [A]−1−1

対角スケーリング：簡単＝弱い[M ] = [D] , [M ] = [D]−1−1

OMP-141

対角スケーリング，点ヤコビ前処理• 前処理行列として，もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする。
– 対角スケーリング，点ヤコビ（point-Jacobi）前処理
 D1
0

[M ] =  ...

0
 0

0 ...
0
D2
0
...
DN −1
0
0 ...
0

0 

0 
... 

0 
DN 

• solve [M]z(i-1)= r(i-1)という場合に逆行列を簡単に求めることができる。
• 簡単な問題では収束する。

42

A64FXプロセッサ
TofuD コントローラ TofuD コントローラ

0

1

2

3

4

5

6

7

8

9

10 11

0

1

2

3

4

5

6

7

8

9

10 11

CMG#1

リ
ン
グ
パ
ス

0

1

2

3

4

5

6

7

8

9

10 11

0

1

2

3

4

5

6

7

8

9

10 11

CMG#3

HBM2

HBM2

CMG#2
HBM2

HBM2

CMG#0

プロセッサ名

A64FXプロセッサ数
1 (48+アシスタントコア2 or 4)(コア数)
2.2 GHz
周波数3.3792 TFLOPS
理論演算性能32 GiB
メモリ容量
1,024 GB/s
メモリ帯域幅
L1 Cache
64 KiB/core (Inst/Data)
L2 Cache
8 MiB/CMG• 4つのCMG（Core Memory Group），12計算コア/CMG• NUMAアーキテクチャ（Non-Uniform Memory Access）
 メモリは各CMGに搭載されていて独立，異なるCMGのローカ
ルメモリ上のデータをアクセスすることは可能 ローカルメモリ上のデータを使って計算するのが効率的• 大規模並列：各CMGに1-MPIプロセス（12-OpenMPスレッド），プロセッサ内4プロセスのハイブリッド推奨

43

A64FXプロセッサ
CMG#2

1 2
0MPI
#0 3

1 2#230 MPI

5 6 7
124 Threads
8

9

10 11

0

1

2

3

#1 7
5 6
4MPI

リ
ン
グ
パ
ス

5 6 7
124 Threads
8

9

10 11

0

1

2

3

5 6#37
4 MPI

128 Threads
9 10 11

128 Threads
9 10 11

CMG#1

CMG#3

HBM2

HBM2

CMG#0

プロセッサ名HBM2

HBM2

TofuD コントローラ TofuD コントローラ

A64FX

プロセッサ数
1 (48+アシスタントコア2 or 4)(コア数)
2.2 GHz
周波数3.3792 TFLOPS
理論演算性能32 GiB
メモリ容量
1,024 GB/s
メモリ帯域幅
L1 Cache
64 KiB/core (Inst/Data)
L2 Cache
8 MiB/CMG• 4つのCMG（Core Memory Group），12計算コア/CMG• NUMAアーキテクチャ（Non-Uniform Memory Access）
 メモリは各CMGに搭載されていて独立，異なるCMGのローカ
ルメモリ上のデータをアクセスすることは可能 ローカルメモリ上のデータを使って計算するのが効率的• 大規模並列：各CMGに1-MPIプロセス（12-OpenMPスレッド），プロセッサ内4プロセスのハイブリッド推奨

44

A64FX：CMG（Core Memory Group）
TofuD コントローラ TofuD コントローラ
CMG#2

HBM2

1

2

3

0

1

2

3

4

5

6

7

4

5

6

7

8

9

10 11

8

9

10 11

0

1

2

3

0

1

2

3

4

5

6

7

4

5

6

7

8

9

10 11

8

9

10 11

CMG#1

リ
ン
グ
パ
ス

CMG#3

HBM2

0

HBM2

HBM2

CMG#0

45

CMG番号，コア番号，メモリ番号（1/2）
CMG#2

12 13 14 15

36 37 38 39

16 17 18 19

40 41 42 43

20 21 22 23
24 25 26 27
28 29 30 31

リ
ン
グ
パ
ス

44 45 46 47
48 49 50 51
52 53 54 55

32 33 34 35

56 57 58 59

CMG#1

CMG#3

Memory #7

Memory #5

CMG#0

Memory #6

Memory #4

CMG:#0-#3，Core:#12-59，Memory:#4-#7

46

CMG番号，コア番号，メモリ番号（1/2）
CMG#0

CMG#2

12 13 14 15

36 37 38 39

16 17 18 19

40 41 42 43

20 21 22 23

Memory #5

44 45 46 47
48 49 50 51
52 53 54 55

32 33 34 35

56 57 58 59

CMG#1

CMG#3

Memory #7

まずはこの12コア，1-CMG
24 25 26 27だけ使ってやってみる
28 29 30 31

リ
ン
グ
パ
ス

Memory #6

Memory #4

CMG:#0-#3，Core:#12-59，Memory:#4-#7

OMP-147

ファイルコピー：Wisteria/BDEC-01>$ cd /work/gt00/t00XYZ>$ cp /work/gt00/z30088/ompw.tar .
>$ tar xvf ompw.tar>$ cd ompw

これを<$O-ompw>と呼ぶ

以下のディレクトリが出来ていることを確認
run src-c0 src-c1 src-c2 src-c3 src-c3b src-f0 src-f1 src-f2 src-f3 src-f4>$ module load fj ログインしたら必ずこれをタイプ（コンパイラ）>$ make –f makec>$ ls run/solc0
solc0>$ cd run<modify “INPUT.DAT”, “c12.sh”>>$ pjsub c12.sh

48<$O-ompw>/makec
default:

clean:(cd src-c0 ; make )(cd src-c1 ; make )(cd src-c2 ; make )(cd src-c3 ; make )(cd src-c3b ; make )(cd src-c0 ; make clean)(cd src-c1 ; make clean)(cd src-c2 ; make clean)(cd src-c3 ; make clean)(cd src-c3b ; make clean)<$O-ompw>/makec-org
default:

clean:(cd src-c0 ; make –f make-o)(cd src-c1 ; make –f make-o)(cd src-c2 ; make –f make-o)(cd src-c3 ; make –f make-o)(cd src-c0 ; make –f make-o clean)(cd src-c1 ; make –f make-o clean)(cd src-c2 ; make –f make-o clean)(cd src-c3 ; make –f make-o clean)

49<$O-ompw>/src-c0/Makefile, make-oparallel computing by OpenMP
MakefileCC= fccpx
OPTFLAG= -Kfast,openmp -Nclang -msve-vector-bits=512 -ffj-ocl
TARGET= ../run/solc0.SUFFIXES:.SUFFIXES: .o .c.c.o:$(CC) -c $(CFLAGS) $(OPTFLAG) $< -o $@OBJS = input.o pointer_init.o boundary_cell.o cell_metrics.o ¥poi_gen.o solver_PCG.o outucd.o allocate.o main.o
HEADERS = struct.h struct_ext.h pcg.h pcg_ext.h input.h pointer_init.h¥boundary_cell.h cell_metrics.h poi_gen.h solver_PCG.h allocate.h outucd.h
all: $(TARGET)$(TARGET): $(OBJS)$(CC) $(CFLAGS) $(OPTFLAG) -o $@ $(OBJS)$(OBJS): $(HEADERS)
clean:

make-orm -f *.o $(TARGET) *.log *~ *.lst

CC= fccpx
OPTFLAG= -Kfast,openmp
TARGET = ../run/solc0.SUFFIXES:.SUFFIXES: .o .c
.c.o:
(…)

$(CC) -c $(CFLAGS) $(OPTFLAG) $< -o $@OMP-350

Cコンパイラ：2種類のモード

tradモード(-Nnoclang オプション)（デフォルト）• 京およびPRIMEHPC FX100以前のシステム向け富士通コンパイラをベースとする。
• tradモードは，従来の富士通コンパイラとの互換を重視する場合に適している。
• サポートしている仕様は，C89/C99/C11，OpenMP
3.1/OpenMP 4.5(一部)• オプション省略時（デフォルト）は，-Nnoclangオプショ
ン適用• 本ケースの場合遅い

clangモード(-Nclang オプション)• オープンソースソフトウェアであるClang/LLVMコンパ
イラをベースとする。
• clangモードは，最新言語仕様を使用したプログラムや，
オープンソースソフトウェアを翻訳する場合に適している。
• サポートしている仕様は，C89/C99/C11， OpenMP
4.5/OpenMP 5.0(一部)• 本ケースの場合，tradモードより速い，最適化すると差は縮まるが，今回はデフォルトでclangモード使用

51

ジョブ実行
• 実行方法– 基本的にバッチジョブのみ– インタラクティヴの実行は「基本的に」できません• 実行手順– ジョブスクリプトを書きます– ジョブを投入します– ジョブの状態を確認します– 結果を確認します• その他– 実行時には1ノード（48コア）が占有されます– 他のユーザーのジョブに使われることはありません

52

CMG番号，コア番号，メモリ番号（1/2）
CMG#0

CMG#2

12 13 14 15

36 37 38 39

16 17 18 19

40 41 42 43

20 21 22 23

Memory #5

44 45 46 47
48 49 50 51
52 53 54 55

32 33 34 35

56 57 58 59

CMG#1

CMG#3

Memory #7

まずはこの12コア，1-CMG
24 25 26 27だけ使ってやってみる
28 29 30 31

リ
ン
グ
パ
ス

Memory #6

Memory #4

CMG:#0-#3，Core:#12-59，Memory:#4-#7

ジョブスクリプト：c12.sh53• /work/gt00/t00XXX/ompw/run/c12.sh• スケジューラへの指令 ＋ シェルスクリプト#!/bin/sh#PJM -N “c12”#PJM -L rscgrp=tutorial-o#PJM -L node=1#PJM --omp thread=12#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o c12.lst

ジョブ名称（省略可）
実行キュー名（Resource Group）
ノード数（原則＝1）
スレッド数（1-48，しばらくは1-12）
実行時間
グループ名（財布）C
エラー出力ファイル
標準出力ファイル

module load fj
export OMP_NUM_THREADS=12
スレッド数（--omp thread=XXと同じ数）
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactl –l ./solc0
numactl –C 12-23 –m 4 ./solc0

ジョブスクリプト：c12.sh• /work/gt00/t00XXX/ompw/run/c12.sh• numactl：ローカルなリソースを使う

CMG#0

CMG#2

12 13 14 15

36 37 38 39

16 17 18 19
20 21 22 23

Memory #5

40 41 42 43
44 45 46 47
48 49 50 51
52 53 54 55

32 33 34 35

56 57 58 59

CMG#1

CMG#3

Memory #7

まずはこの12コア，1-CMG
24 25 26 27だけ使ってやってみる
28 29 30 31

リ
ン
グ
パ
ス

Memory #6

Memory #4

numactl –l ./solc0
numactl –C 12-23 –m 4 ./solc0

54

55

ジョブスクリプト：c12.shexport XOS_MMM_L_PAGING_POLICY=
demand:demand:demand
環境変数名

指定値（＿はdefault）

説明

XOS_MMM_L_PAGING_POLICY[demand | prepage][demand | prepage][demand | prepage]

各メモリ領域のページング方式（ページの割り当て契機）を選択する設定
 demand：デマンドページング方式
 prepage：プリページング方式
本変数はコロン区切りで 3 つのメモリ領域のページング方式を指定：
 第1指定：静的データの .bss 領域，静的データの .data 領域はページング方式指定の対象外で常に prepage
 第2指定：スタック領域・スレッドスタック領域
 第3 指定は動的メモリ確保領域
指定値以外の値を指定した場合は，「prepage:demand:prepage 」指定とみなす。
複数CMGをまたぐ場合はdemandを推奨！

56

ジョブ投入
>$ cd /work/gt00/t00XYZ>$ cd ompw/run>$ pjsub c12.sh>$ cat c12.lst
INPUT.DAT
128 128 128
1.00e-0 1.00e-00 1.00e-00
1.0e-08

NX NY NZ
DX/DY/DZEPSICCG出力（1/2）[t00XYZ@wisteria01 run]$ cat c12.lst
128
128
128
1
8.958216e+00
101
8.313496e+00
201
2.090443e+00
301
3.811029e-01
401
3.769653e-02
501
9.429978e-04
601
4.940783e-05
701
1.888611e-06
801
2.243179e-08
826
9.818026e-09
8.634058e+00 sec. (solver)
2097152
1.459831e+04
128
128
128
1
8.958216e+00
101
8.313496e+00
201
2.090443e+00
301
3.811029e-01
401
3.769653e-02
501
9.429978e-04
601
4.940783e-05
701
1.888611e-06
801
2.243179e-08
826
9.818026e-09
8.631781e+00 sec. (solver)
2097152
1.459831e+04

57

58

出力（2/2）
5回ずつ実施：ほとんど同じ[t00XYZ@wisteria01 run]$ grep “(sol” c12.lst
5.040907e+00 sec. (solver)
5.040407e+00 sec. (solver)
5.038814e+00 sec. (solver)
5.042296e+00 sec. (solver)
5.033407e+00 sec. (solver)
5.034004e+00 sec. (solver)
5.035210e+00 sec. (solver)
5.047046e+00 sec. (solver)
5.043275e+00 sec. (solver)
5.036129e+00 sec. (solver)

numactl ./solc0

numactl –C 12-23 –m 4 ./solc0

59

利用可能なResource Group（Queue）• 以下の2種類のResource Groupを利用可能• 最大12ノードを使える– lecture-o• 12ノード（576コア），15分，アカウント有効期間中利用可能• 全教育ユーザーで共有– tutorial-o• 12ノード（576コア），15分，講義・演習実施時間帯• Lecture-oよりは多くのジョブを投入可能（混み具合による）60

様々なコマンド• ジョブ実行• ジョブ実行状況• ジョブ停止• ジョブキューの状況• ジョブキューの状況（詳細）• 実行ジョブ情報• ジョブ実行履歴• ジョブ実行制限

pjsub SCRIPT NAME
pjstatpjdel JOB ID
pjstat --rsc
pjstat --rsc –x
pjstat –a
pjstat –H
pjstat --limit

61[t00470@wisteria01 run]$ pjsub f2_48.sh[INFO] PJM 0000 pjsub Job 15713 submitted.
[t00470@wisteria01 run]$ pjsub f3_48.sh[INFO] PJM 0000 pjsub Job 15714 submitted.
[t00470@wisteria01 run]$ pjstat
Wisteria/BDEC-01 scheduled stop time: 2021/05/28(Fri) 09:00:00 (Remain:
JOB_ID
15713
15714

JOB_NAME
f2_48
f3_48

STATUS PROJECTRUNNING gt00
QUEUED gt00

RSCGROUP
lecture-o
lecture-o

START_DATE
05/24 07:34:03--/-- --:--:--[t00470@wisteria01 run]$ pjstat
Wisteria/BDEC-01 scheduled stop time: 2021/05/28(Fri) 09:00:00 (Remain:
JOB_ID
15713
15714

JOB_NAME
f2_48
f3_48

STATUS PROJECTRUNNING gt00
RUNNING gt00

RSCGROUP
lecture-o
lecture-o

4days 1:25:56)
ELAPSE
00:00:02
00:00:00

TOKEN-

NODE GPU
1
1
-

TOKEN
-

NODE GPU
1
1
-

TOKEN
-

NODE GPU1-

4days 1:25:56)

START_DATE
05/24 07:34:03(05/24 07:34)

ELAPSE
00:00:02
00:00:00[t00XYZ@wisteria01 ~]$ pjdel 15714[INFO] PJM 0100 pjdel Accepted Job 15714[t00XYZ@wisteria01 ~]$ pjstat
Wisteria/BDEC-01 scheduled stop time: 2021/05/28(Fri) 09:00:00 (Remain:
JOB_ID
15713

JOB_NAME
f2_48

STATUS PROJECTRUNNING gt00

RSCGROUP
lecture-o

START_DATE
05/24 07:34:03[t00XYZ@wisteria01 ~]$ pjstat
Wisteria/BDEC-01 scheduled stop time: 2021/05/28(Fri) 09:00:00 (Remain:
No unfinished job found.

4days 1:25:56)
ELAPSE
00:00:02

4days 1:21:45)

62[t00XYZ@wisteria01 ~]$ pjstat --rsc
SYSTEM: Odyssey
RSCGRPSTATUS
lecture-o[ENABLE,START]
tutorial-o[DISABLE,STOP]SYSTEM: Aquarius
RSCGRP
lecture-atutorial-a

NODE
96
2x12x16STATUS[ENABLE,START][DISABLE,STOP]NODE GPU
7 56
2 16[t00XYZ@wisteria01 ~]$ pjstat --rsc -x
SYSTEM: Odyssey
RSCGRPSTATUS
lecture-o[ENABLE,START]
tutorial-o[DISABLE,STOP]MIN_NODE MAX_NODE MAX_ELAPSE REMAIN_ELAPSE MEM(GiB) PROJECT
1
12
00:15:00
00:15:00
28 gt00
1
12
00:15:00--:--:-28 gt00

SYSTEM: Aquarius
RSCGRP
lecture-atutorial-a

MIN_NODE MAX_NODE AVAIL_GPU MAX_ELAPSE REMAIN_ELAPSE MEM(GiB) PROJECT
1
1
1,2,4
00:15:00
00:15:00
448 gt00
1
1
1,2,4
00:15:00--:--:-448 gt00

STATUS[ENABLE,START][DISABLE,STOP][t00XYZ@wisteria01 ~]$ pjstat --limit
SYSTEM: Odyssey
PROJECTACCEPT
RUN BULK_ACCEPT
gt00
0/ 1280/ 160/
8BULK_RUN
0/
16

NODE
0/ 2304

SYSTEM: Aquarius
PROJECTACCEPTgt00
0/
4

BULK_RUN
0/
0

GPU
64

RUN
0/
2

BULK_ACCEPT
0/
0

0/

poi_gen (1/2)：主要部
とにかくprivateに注意#pragma omp parallel for private(icel,icN1,icN2,icN3,icN4,icN5,icN6,VOL0,isLU,icou,coef, ii, jj, kk)
for(icel=0; icel<ICELTOT; icel++) {icN1 = NEIBcell[icel][0];icN2 = NEIBcell[icel][1];
icN3 = NEIBcell[icel][2];
icN4 = NEIBcell[icel][3];icN5 = NEIBcell[icel][4];
icN6 = NEIBcell[icel][5];VOL0 = VOLCEL[icel];
isLU=indexLU[icel];
icou= 0;if(icN5 != 0) {coef = RDZ * ZAREA;
D[icel] -= coef;
itemLU[icou+isLU]= icN5-1;
AMAT [icou+isLU]= coef;
icou= icou + 1;}(…)

ii = XYZ[icel][0];jj = XYZ[icel][1];kk = XYZ[icel][2];}BFORCE[icel] = - (double)(ii + jj + kk) * VOLCEL[icel];63poi_gen (2/2)：境界条件部
とにかくprivateに注意#pragma omp parallel for private (ib,icel,coef)for(ib=0; ib<ZmaxCELtot; ib++) {icel = ZmaxCEL[ib];coef = 2.0 * RDZ * ZAREA;
D[icel-1] -= coef;
}

64

65

#include <stdio.h>#include <stdlib.h>#include <string.h>#include <errno.h>#include <math.h>#include <omp.h>

solve_PCG (1/5)
parallel computing by OpenMP
若干インタフェース変更
restrict：ポインタが他の変数から参照されない#include "solver_PCG.h"
extern int
solve_PCG(int N, int *restrict indexLU, int *restrict itemLU,
double *restrict D, double *restrict B, double *restrict X,
double *restrict AMAT, double EPS, int *restrict ITR, int *restrict IER){
double VAL, BNRM2, WVAL, SW, RHO, BETA, RHO1, C1, DNRM2, ALPHA, ERR;
double Stime, Etime;
int i, j, ic, ip, L, ip1, N3;
int R = 0;
int Z = 1;
int Q = 1;
int P = 2;
int DD = 3;
double (*restrict W)[N] = (double (*)[N])malloc(4*sizeof(double[N]));if(W == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}

solve_PCG (2/5)#pragma omp parallel for private (i)for(i=0; i<N; i++) {X[i] = 0.0;
W[1][i] = 0.0;
W[2][i] = 0.0;
W[3][i] = 0.0;}
#pragma omp parallel for private (i)for(i=0; i<N; i++) {W[DD][i]= 1.e0/D[i];}#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * X[i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * X[itemLU[j]];}W[R][i] = B[i] - VAL;}

66

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

BNRM2 = 0.0;
#pragma omp parallel for private (i) reduction (+:BNRM2)
for(i=0; i<N; i++) {BNRM2 += B[i]*B[i];}solve_PCG (3/5)*ITR = N;

67

Compute r(0)= b-[A]x(0)for i= 1, 2, …
Stime = omp_get_wtime();solve [M]z(i-1)= r(i-1)
for(L=0; L<(*ITR); L++) {ρi-1= r(i-1) z(i-1)
if i=1#pragma omp parallel for private (i)p(1)= z(0)for(i=0; i<N; i++) {W[Z][i] = W[R][i]*W[DD][i];else}
βi-1= ρi-1/ρi-2p(i)= z(i-1) + βi-1 p(i-1)
RHO = 0.0;
#pragma omp parallel for private (i) reduction(+:RHO) endif
for(i=0; i<N; i++) {q(i)= [A]p(i)RHO += W[R][i] * W[Z][i];αi = ρi-1/p(i)q(i)}x(i)= x(i-1) + αip(i)if(L == 0) {r(i)= r(i-1) - αiq(i)#pragma omp parallel for private (i)check convergence |r|for(i=0; i<N; i++) {end
W[P][i] = W[Z][i];}} else {BETA = RHO / RHO1;#pragma omp parallel for private (i)for(i=0; i<N; i++) {W[P][i] = W[Z][i] + BETA * W[P][i];}}solve_PCG (4/5)#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}

68

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)elseC1 = 0.0;#pragma omp parallel for private (i) reduction(+:C1)
βi-1= ρi-1/ρi-2
for(i=0; i<N; i++) {p(i)= z(i-1) + βi-1 p(i-1)
C1 += W[P][i] * W[Q][i];endif}
ALPHA = RHO / C1;q(i)= [A]p(i)αi = ρi-1/p(i)q(i)#pragma omp parallel for private (i)for(i=0; i<N; i++) {x(i)= x(i-1) + αip(i)X[i]+= ALPHA * W[P][i];(i)= r(i-1) - α q(i)riW[R][i] -= ALPHA * W[Q][i];check convergence |r|}end
DNRM2 == 0.0;
0.0;
DNRM2
#pragma omp
omp parallel
parallel forfor private
private (i)(i) reduction(+:DNRM2)
reduction(+:DNRM2)#pragmafor(i=0; i<N;
i<N; i++)i++) {{
for(i=0;
DNRM2 +=+= W[R][i]*W[R][i];W[R][i]*W[R][i];DNRM2}}ERR = sqrt(DNRM2/BNRM2);...

solve_PCG (5/5)
Stime = omp_get_wtime();for(L=0; L<(*ITR); L++) {...
if(ERR < EPS) {*IER = 0;
goto N900;} else {RHO1 = RHO;}}*IER = 1;
N900:
Etime = omp_get_wtime();Elapsed Time= Etime - Stime

fprintf(stderr, "%5d%16.6e¥n", L+1, ERR);
fprintf(stderr, "%16.6e sec. (solver)¥n", Etime - Stime);*ITR = L;free(W);
return 0;
}

69

70

c01,c02,c04,c06,c08,c12.sh• /work/gt00/t00XXX/ompw/run/cYZ.sh• スケジューラへの指令 ＋ シェルスクリプト#!/bin/sh#PJM -N “cYZ”#PJM -L rscgrp=tutorial-o#PJM -L node=1#PJM --omp thread=YZ#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o cYZ.lst

ジョブ名称（省略可）
実行キュー名（Resource Group）
ノード数（原則＝1）
スレッド数（1-48，しばらくは1-12）
実行時間
グループ名（財布）
エラー出力ファイル
標準出力ファイル

module load fj
export OMP_NUM_THREADS=YZ
スレッド数（--omp thread=YZと同じ数）
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactl –l ./solc0
numactl –C 12-23 –m 4 ./solc0

71

PCG計算時間: Etime-Stime：Fortran
NX=NY=NZ=128
実行時間：5回測定して最速時間採用
コア数増加による効率低下：メモリ利用効率の低下
Thread #sec

Speed-up

ParallelEfficiency (%)

Large
大

1

50.27

1.00

100.00

2

25.24

1.99

99.60

Granularity：粒度

4

12.98

3.87

96.86

Problem Size/Thread

6

9.24

5.44

90.73

Small
小

8

7.27

6.92

86.50

12

5.09

9.88

82.30

Parallel Efficiency(%)= 100*(Speed-Up)/Thread#72

c04.sh#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

“c04"
rscgrp=tutorial-o
node=1
thread=4
elapse=00:15:00
gt00
err
c04.lstmodule load
fjexport OMP_NUM_THREADS=4
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0-l ./solc0-l ./solc0-l ./solc0-l ./solc0
-C
12-23
-C
12-23
-C
12-23
-C
12-23
-C
12-23

-m
-m
-m
-m
-m

4
4
4
4
4

./solc0./solc0./solc0./solc0./solc0

73

c08.sh#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

“c08"
rscgrp=tutorial-o
node=1
thread=8
elapse=00:15:00
gt00
err
c08.lstmodule load
fjexport OMP_NUM_THREADS=8
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0-l ./solc0-l ./solc0-l ./solc0-l ./solc0
-C
12-23
-C
12-23
-C
12-23
-C
12-23
-C
12-23

-m
-m
-m
-m
-m

4
4
4
4
4

./solc0./solc0./solc0./solc0./solc0

74

CMG数増加PCG計算時間: Etime-Stime：Fortran
NX=NY=NZ=128
実行時間：5回測定して最速時間採用
Large
大Thread #sec

Speed-up

ParallelEfficiency (%)

Granularity：粒度

12

5.09

12.00

100.00Problem Size/Thread

24

2.79

21.88

91.18

Small
小

36

1.99

30.75

85.41

48

1.70

35.97

74.95

Parallel Efficiency(%)= 100*(Speed-Up)/Thread#

75

CMG番号，コア番号，メモリ番号（1/2）
CMG#2

12 13 14 15

36 37 38 39

16 17 18 19

40 41 42 43

20 21 22 23
24 25 26 27
28 29 30 31

リ
ン
グ
パ
ス

44 45 46 47
48 49 50 51
52 53 54 55

32 33 34 35

56 57 58 59

CMG#1

CMG#3

Memory #7

Memory #5

CMG#0

Memory #6

Memory #4

CMG:#0-#3，Core:#12-59，Memory:#4-#7

76

CMG番号，コア番号，メモリ番号（2/2）
CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

42 43
40-C4136-47

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

Memory #7

Memory #5

CMG#0

Memory #6

Memory #4

CMG:#0-#3，Core:#12-59，Memory:#4-#7

77

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21

Memory #5

err
c0_12.lst

CMG#0

25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=12
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0-l ./solc0-l ./solc0-l ./solc0-l ./solc0
-C
12-23
-C
12-23
-C
12-23
-C
12-23
-C
12-23

-m
-m
-m
-m
-m

4
4
4
4
4

./solc0./solc0./solc0./solc0./solc0

Memory #7“c0_12”
rscgrp=tutorial-o
node=1
thread=12
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c0_12.sh78

Memory #4
Memory #5

err
err
c0_24.lst
c0_24.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

-m 34
5 35
32 33

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 58
7 59
56 57

CMG#1
CMG#3
moduleload
fj
module load
fjexport
OMP_NUM_THREADS=24
export OMP_NUM_THREADS=24
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
XOS_MMM_L_PAGING_POLICY=demand:demand:demand
exportnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0./solc0-l-l ./solc0./solc0-l-l ./solc0./solc0-l-l ./solc0./solc0-l-l ./solc0./solc0
-l
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35
-C
12-35

-m
-m
-m
-m
-m
-m
-m
-m
-m
-m

4-5
4-5
4-5
4-5
4-5
4-5
4-5
4-5
4-5
4-5

./solc0./solc0./solc0./solc0./solc0./solc0./solc0./solc0./solc0./solc0

Memory #7“c0_24”“c0_24”
rscgrp=tutorial-o
rscgrp=tutorial-o
node=1
node=1
thread=24
thread=24
elapse=00:15:00
elapse=00:15:00
gt00
gt00

Memory #6#!/bin/sh#!/bin/sh#PJM -N-N#PJM#PJM -L-L#PJM#PJM -L-L#PJM#PJM --omp
--omp
#PJM
#PJM -L-L#PJM#PJM -g-g#PJM#PJM -j-j#PJM#PJM -e-e#PJM#PJM -o-o#PJMc0_24.sh79

Memory #5

err
c0_36.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=36
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0-l ./solc0-l ./solc0-l ./solc0-l ./solc0
-C
12-47
-C
12-47
-C
12-47
-C
12-47
-C
12-47

-m
-m
-m
-m
-m

4-6
4-6
4-6
4-6
4-6

./solc0./solc0./solc0./solc0./solc0

Memory #7“c0_36”
rscgrp=tutorial-o
node=1
thread=36
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c0_36.sh80

Memory #5

err
c0_48.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc0-l ./solc0-l ./solc0-l ./solc0-l ./solc0
-C
12-59
-C
12-59
-C
12-59
-C
12-59
-C
12-59

-m
-m
-m
-m
-m

4-7
4-7
4-7
4-7
4-7

./solc0./solc0./solc0./solc0./solc0

Memory #7“c0_48”
rscgrp=tutorial-o
node=1
thread=48
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c0_48.sh81• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

82

A64FXプロセッサ
TofuD コントローラ TofuD コントローラ

0

1

2

3

4

5

6

7

8

9

10 11

0

1

2

3

4

5

6

7

8

9

10 11

CMG#1

リ
ン
グ
パ
ス

0

1

2

3

4

5

6

7

8

9

10 11

0

1

2

3

4

5

6

7

8

9

10 11

HBM2

HBM2

CMG#2
HBM2

HBM2

CMG#0

プロセッサ名

A64FXプロセッサ数
1 (48+アシスタントコア2 or 4)(コア数)
2.2 GHz
周波数3.3792 TFLOPS
理論演算性能32 GB
メモリ容量1,024 GB/s
メモリ帯域幅

CMG#3• 4つのCMG（Core Memory Group），12計算コア/CMG• NUMAアーキテクチャ（Non-Uniform Memory Access）
 メモリは各CMGに搭載されていて独立，異なるCMGのローカ
ルメモリ上のデータをアクセスすることは可能 ローカルメモリ上のデータを使って計算するのが効率的• 大規模並列：各CMGに1-MPIプロセス（12-OpenMPスレッド），プロセッサ内4プロセスのハイブリッド推奨OMP-383

データをローカルメモリに置く方法（1/2）• NUMAアーキテクチャでは，プログラムにおいて変数や配列を宣言した時点では，物理的メモリ上に記憶領域は確保されず，ある変数を最初にアクセスしたコア（の属するCMG）のローカルメモリ上に，その変数の記憶領域（ページ）が確保される。
• これをFirst Touch Data Placement（First Touch）と呼び，配列の初期化手順により得られる性能が大幅に変化する場合があるため，注意が必要である。
• 例えばある配列を初期化する場合，特に指定しなければ0番のCMGで初期化が行われるため，記憶領域は0番
CMGのローカルメモリ上に確保される。

OMP-384

データをローカルメモリに置く方法（2/2）• したがって，他のCMGでこの配列のデータをアクセスする
場合には，必ず0番CMGのメモリにアクセスする必要があるため，高い性能を得ることは困難である。
• 配列の初期化を，実際の計算の手順にしたがってOpenMPを使って並列に実施すれば，実際に計算を担当
するCMGのメモリにその配列の担当部分の記憶領域が確保され，より効率的に計算を実施することができる。
• 1CMGしか使用しない場合はこのような配慮は不要– OpenMP/MPIハイブリッドで1-CMG当りに1つのMPIプロセスを使用する場合も同様– numactl

OMP-385

First Touch Data Placement“Patterns for Parallel Programming” Mattson, T.G. et al.
To reduce memory traffic in the system, it is important to keep thedata close to the PEs that will work with the data (e.g. NUMA control).
On NUMA computers, this corresponds to making sure the pages ofmemory are allocated and “owned” by the PEs that will be working
with the data contained in the page.
The most common NUMA page-placement algorithm is the“first touch” algorithm, in which the PE first referencing a region ofmemory will have the page holding that memory assigned to it.
A very common technique in OpenMP program is to initialize
data in parallel using the same loop schedule as will be used laterin the computations.

OMP-1

86

First-Touchに対応したバージョン：src-c1
オリジナルは src-c0>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “INPUT.DAT”, “c1_XY.sh”> (XY:12,24,36,48)>$ pjsub c1_XY.shOMP-187[XYZ@wisteria01 run]$ cd ../src-c0[XYZ@wisteria01 src-f0]$ diff poi_gen.c ../src-c1/poi_gen.c
25,29c25,31---

for (i = 0; i <ICELTOT ; i++) {BFORCE[i]=0.0;
D[i] =0.0;
PHI[i]=0.0;
INLU[i] = 0;}

#pragma omp parallel for private (i)for (i = 0; i <ICELTOT ; i++) {BFORCE[i]=0.0;
D[i] =0.0;
PHI[i]=0.0;
INLU[i] = 0;}

src-c0src-c1

OMP-188

for (i = 0; i <=ICELTOT ; i++) {indexLU[i] = 0;}

src-c0--#pragma omp parallel for private (i)for (i = 0; i <=ICELTOT ; i++) {indexLU[i] = 0;}---

src-c1for(i=0; i<ICELTOT; i++) {for(j=indexLU[i]; j<indexLU[i+1]; j++) {itemLU[j]=0;
AMAT[j]=0.0;}}#pragma omp parallel for private (i,j)
for(i=0; i<ICELTOT; i++) {for(j=indexLU[i]; j<indexLU[i+1]; j++) {itemLU[j]=0;
AMAT[j]=0.0;}}

src-c0src-c1

89

Memory #5

err
c1_48.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc1-l ./solc1-l ./solc1-l ./solc1-l ./solc1
-C
12-59
-C
12-59
-C
12-59
-C
12-59
-C
12-59

-m
-m
-m
-m
-m

4-7
4-7
4-7
4-7
4-7

./solc1./solc1./solc1./solc1./solc1

Memory #7“c1_48”
rscgrp=tutorial-o
node=1
thread=48
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c1_48.shCMG数増加
PCG計算時間: Etime-Stime：Fortran
NX=NY=NZ=128
実行時間：5回測定して最速時間採用

src-f0src-f1Thread #sec

Speed-up

ParallelEfficiency (%)

12

5.09

12.00

100.00

24

2.79

21.88

91.18

36

1.99

30.75

85.41

48

1.70

35.97

74.95

12

5.27

-

-

24

2.70

22.61

94.20

36

1.87

32.65

90.69

48

1.52

40.09

83.51

90

91

CMG数増加PCG計算時間: Etime-Stime：C（clang）
NX=NY=NZ=128
実行時間：5回測定して最速時間採用

src-c0src-c1

Parallel
Speed-up
Efficiency (%)

Thread #

sec

12

5.03

12.00

100.00

24

2.69

22.42

93.41

36

1.91

31.57

87.69

48

1.57

38.51

80.23

12

5.22

-

-

24

2.70

22.39

93.29

36

1.85

32.66

90.71

48

1.17

51.70

107.7

Cache is wellutilized,
because theproblem size is
small

92

CMG数増加PCG計算時間: Etime-Stime
NX=NY=NZ=128
実行時間：5回測定して最速時間採用

Original

First-Touch

Language

12

24

36

48

Fortran

5.09

2.79

1.99

1.70

C (clang)

5.03

2.69

1.91

1.57

C (trad)

7.75

4.19

2.90

2.36

Fortran

5.27

2.70

1.87

1.52

C (clang)

5.22

2.70

1.85

C (trad)

7.85

4.05

2.79

1.17
1.72
(112.6 %)93• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

94

疎行列の格納方式
CRS（Compressed Row Storage）
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}ELL (ELLPACK/ITPACK)
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[6*i+j] * W[P][itemLU[6*i+j]];}
W[Q][i] = VAL;}• CRS: Compressed Row Storage– 非ゼロ非対角成分のみ格納⇒メモリ節約できるが，計算効率悪い• ELL：ELLPACK/ITPACK– 非ゼロ非対角成分数固定⇒0のところには0を入れる– 記憶容量，計算量は増えるがメモリアクセス性能向上⇒Prefetch

95

疎行列の格納方式
CRS（Compressed Row Storage）
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}ELL (ELLPACK/ITPACK)
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[6*i+j] * W[P][itemLU[6*i+j]];}
W[Q][i] = VAL;}ELL (ELLPACK/ITPACK) こちらの方がだいぶ遅い
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[i][j] * W[P][itemLU[i][j]];}W[Q][i] = VAL;}

96

poi_gen (1/2)：Privateに注意#pragma omp parallel for private(icel,icN1,icN2,icN3,icN4,icN5,icN6,VOL0,icou,coef,ii,jj,kk)
for(icel=0; icel<ICELTOT; icel++) {icN1 = NEIBcell[icel][0];icN2 = NEIBcell[icel][1];
icN3 = NEIBcell[icel][2];
icN4 = NEIBcell[icel][3];icN5 = NEIBcell[icel][4];
icN6 = NEIBcell[icel][5];(…)

VOL0 = VOLCEL[icel];
icou= 0;if(icN5 != 0) {coef = RDZ * ZAREA;
D[icel] -= coef;
itemLU[6*icel+icou]= icN5-1;
AMAT [6*icel+icou]= coef;
icou= icou + 1;}
if(icN3 != 0) {coef = RDZ * YAREA;
D[icel] -= coef;
itemLU[6*icel+icou]= icN3-1;
AMAT [6*icel+icou]= coef;
icou= icou + 1;
}
}

97

poi_gen (2/2)：Padding
N2= 128 (poi_gen.c)

N

N+1

N+2

N+3

N+4

N+6

N+7

N+8

N=ICELTOT/************ PADDING ************/
icou= 0;for (i = 0; i <ICELTOT ; i++) {for (k = 0; k <6 ; k++) {if (itemLU[6*i+k]==-1) {icou= icou + 1;
itemLU[6*i+k]= ICELTOT-1 + icou;
if (icou==N2){icou=0;}
}
}
}

N+5

N+9 N+10 N+11
N+12 N+13 N+14 N+15
N+16 N+17 N+18 N+19

N
N
N+125 +126 +127
N

N+1

N+2

N+3

N+4

N+5

N+6

N+7

N+8

N+9 N+10

98

solve_PCG (1/2)
W = (double **)malloc(sizeof(double *)*4);
if(W == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}
for(i=0; i<4; i++) {W[i] = (double *)malloc(sizeof(double)*(N+N2));
if(W[i] == NULL) {fprintf(stderr, "Error: %s¥n", strerror(errno));return -1;}}#pragma omp parallel for private (i)for(i=0; i<N; i++) {X[i] = 0.0;
W[1][i] = 0.0;
W[2][i] = 0.0;
W[3][i] = 0.0;}
#pragma omp parallel for private (i)for(i=0; i<N; i++) {W[DD][i] = 1.0 / D[i];}99solve_PCG (2/2)/***************** {q} = [A]{p} *****************/#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[6*i+j] * W[P][itemLU[6*i+j]];}
W[Q][i] = VAL;}OMP-1100

ELLに対応したバージョン
src-c2>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “INPUT.DAT”, “c2_48.sh”>>$ pjsub c2_48.sh101

Memory #5

err
c2_48.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc2-l ./solc2-l ./solc2-l ./solc2-l ./solc2
-C
12-59
-C
12-59
-C
12-59
-C
12-59
-C
12-59

-m
-m
-m
-m
-m

4-7
4-7
4-7
4-7
4-7

./solc2./solc2./solc2./solc2./solc2

Memory #7“c2_48”
rscgrp=tutorial-o
node=1
thread=48
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c2_48.sh出力（1/2）Fortran[t00XYZ@wisteria01 run]$ cat f12.lst
1
8.958216E+00
101
8.313496E+00
201
2.090443E+00
301
3.811029E-01
401
3.769653E-02
501
9.429978E-04
601
4.940783E-05
701
1.888611E-06
801
2.243179E-08
826
9.818026E-09
5.275418E+00 sec. (solver)##ANSWER

2097152

1.459831E+04

18.958216E+00
101
8.313496E+00
201
2.090443E+00
301
3.811029E-01
401
3.769653E-02
501
9.429978E-04
601
4.940783E-05
701
1.888611E-06
801
2.243179E-08
826
9.818026E-09
5.270398E+00 sec. (solver)##ANSWER

2097152

1.459831E+04

102

103

出力（2/2）：Fortran，5回実施：ほぼ同じ[XYZ@wisteria01 run]$ grep “(sol” f1_48.lst
1.480524E+00 sec. (solver)
1.501454E+00 sec. (solver)
1.441297E+00 sec. (solver)
1.483405E+00 sec. (solver)
1.481864E+00 sec. (solver)

numactl –l ./solf1

1.475129E+00 sec. (solver)
1.483695E+00 sec. (solver)
1.485036E+00 sec. (solver)
1.502549E+00 sec. (solver)
1.487192E+00 sec. (solver)

numactl –C 12-59 –m 4-7 ./solf1[XYZ@wisteria01 run]$ grep “(sol” f2_48.lst
7.713702E-01 sec. (solver)
7.568300E-01 sec. (solver)
7.328739E-01 sec. (solver)
7.826090E-01 sec. (solver)
7.884219E-01 sec. (solver)

numactl –l ./solf2

7.546160E-01 sec. (solver)
7.937970E-01 sec. (solver)
7.403760E-01 sec. (solver)
7.745121E-01 sec. (solver)
7.862871E-01 sec. (solver)

numactl –C 12-59 –m 4-7 ./solf2

104• OpenMP• Login to Wisteria/BDEC-01• OpenMPによる並列化（0）（12コアまで）• OpenMPによる並列化（1）（First Touch）• OpenMPによる並列化（2）（＋ELL）• OpenMPによる並列化（3）（＋omp-parallel削減）• OpenMPによる並列化（4）（＋更なる最適化（Fortranのみ））

105

omp parallel (do)• omp parallel-omp end
parallelはそのたびにスレッドを生成，消滅させる：fork-join• ループが連続するとオーバーヘッドになる。
• omp parallel + omp
do/omp for
MasterMasterthreadthreadthreadthreadthreadMaster

PARALLEL
fork

threadMasterMasterthreadthreadthreadthreadthreadthreadthreadthreadEND PARALLEL
join

PARALLEL
fork

END PARALLEL
join#pragma omp parallel ...
#pragma omp for {...
#pragma omp for {!$omp parallel ...
!$omp dodo i= 1, N...
!$omp dodo i= 1, N...
!$omp end parallelrequiredomp parallel for：スレッド生成消滅#pragma omp parallel for private (i,VAL,j)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=indexLU[i]; j<indexLU[i+1]; j++) {VAL += AMAT[j] * W[P][itemLU[j]];}W[Q][i] = VAL;}

106

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)elseC1 = 0.0;#pragma omp parallel for private (i) reduction(+:C1)
βi-1= ρi-1/ρi-2
for(i=0; i<N; i++) {p(i)= z(i-1) + βi-1 p(i-1)
C1 += W[P][i] * W[Q][i];endif}
ALPHA = RHO / C1;q(i)= [A]p(i)αi = ρi-1/p(i)q(i)#pragma omp parallel for private (i)for(i=0; i<N; i++) {x(i)= x(i-1) + αip(i)X[i]+= ALPHA * W[P][i];(i)= r(i-1) - α q(i)riW[R][i] -= ALPHA * W[Q][i];check convergence |r|}end
DNRM2 == 0.0;
0.0;
DNRM2
#pragma omp
omp parallel
parallel forfor private
private (i)(i) reduction(+:DNRM2)
reduction(+:DNRM2)#pragmafor(i=0; i<N;
i<N; i++)i++) {{
for(i=0;
DNRM2 +=+= W[R][i]*W[R][i];W[R][i]*W[R][i];DNRM2}}ERR = sqrt(DNRM2/BNRM2);
...

107

更なる最適化の方針• src-c3, src-f3– 各反復で一回だけomp
parallel を呼ぶ

MasterthreadthreadthreadthreadthreadMaster• src-f4（Fortranのみ）– !$omp parallelは一回のみ呼び出す（反復の最初と最後）– !$omp doもなくす– 内積部分のreductionは並
列化しない

Master

PARALLEL
fork

threadMasterMasterthreadthreadthreadthreadthreadthreadthreadthreadEND PARALLEL
join

PARALLEL
fork

END PARALLEL
join

108

更なる最適化の方針：Fortranの場合
src-f3src-f4!$omp parallel private (…)

do L= 1, ITR!$omp parallel private (i,k,VAL)(…)!$omp end parallel
enddo

(…)

900 continue

900 continue

ITR= Ldo L= 1, ITR

enddo

ITR= L!$omp end parallel

deallocate (W)deallocate (W)returnend

returnend

src_c3 (1/2)*ITR = N;
Stime = omp_get_wtime();for(L=0; L<(*ITR); L++) {#pragma omp parallel private (i,j,VAL){#pragma omp for
for(i=0; i<N; i++) {W[Z][i] = W[R][i]*W[DD][i];}RHO = 0.0;
#pragma omp for reduction(+:RHO)
for(i=0; i<N; i++) {RHO += W[R][i] * W[Z][i];}if(L == 0) {#pragma omp for
for(i=0; i<N; i++) {W[P][i] = W[Z][i];}} else {BETA = RHO / RHO1;#pragma omp for
for(i=0; i<N; i++) {W[P][i] = W[Z][i] + BETA * W[P][i];
}
}

109

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

110

src_c3 (2/2)#pragma omp for
Compute r(0)= b-[A]x(0)for(i=0; i<N; i++) {for i= 1, 2, …
VAL = D[i] * W[P][i];for(j=0; j<6; j++) {solve [M]z(i-1)= r(i-1)
VAL += AMAT[6*i+j] * W[P][itemLU[6*i+j]];(i-1) z(i-1)
ρ
=
r
i-1
}
if i=1
W[Q][i] = VAL;}
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|

C1 = 0.0;#pragma omp for reduction(+:C1)
for(i=0; i<N; i++) {C1 += W[P][i] * W[Q][i];}ALPHA = RHO / C1;#pragma omp for
for(i=0; i<N; i++) {X[i]+= ALPHA * W[P][i];W[R][i] -= ALPHA * W[Q][i];}DNRM2 = 0.0;
#pragma omp for reduction(+:DNRM2)
for(i=0; i<N; i++) {DNRM2 += W[R][i]*W[R][i];}}ERR = sqrt(DNRM2/BNRM2);

end

OMP-1111

src-c3 に対応したバージョン>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “INPUT.DAT”, “c3_48.sh”>>$ pjsub c3_48.sh112

Memory #5

err
c3_48.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solc3-l ./solc3-l ./solc3-l ./solc3-l ./solc3
-C
12-59
-C
12-59
-C
12-59
-C
12-59
-C
12-59

-m
-m
-m
-m
-m

4-7
4-7
4-7
4-7
4-7

./solc3./solc3./solc3./solc3./solc3

Memory #7“c3_48”
rscgrp=tutorial-o
node=1
thread=48
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

c3_48.shOMP-1

113

C言語：tradを使う場合>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “INPUT.DAT”, “c48org.sh”>>$ pjsub c48org.sh
c48org.sh#!/bin/sh#PJM -N "cx48"#PJM -L rscgrp=lecture-o#PJM -L node=1#PJM --omp thread=48#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o c48org_160.lst
module load fj
export OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactl –l ./solc0org
numactl –l ./solc0org
numactl –l ./solc0org
numactl –l ./solc0org
numactl –l ./solc0org…

114

PCG法：計算時間：N=1283，48スレッド
N= 2,097,152
FortranC (clang)C (trad)src0（初期設定）

1.671

1.564

2.354

src1（First Touch）

1.480

1.122

1.720

src2（+ELL）

0.747

0.809

1.127

src3
（+omp-parallel削減）

0.707

0.834

0.854

115

PCG法：計算時間：N=1603，48スレッド
N= 4,096,000
FortranC (clang)C (trad)src0（初期設定）

3.610

3.484

4.067

src1（First Touch）

2.993

2.228

3.425

src2（+ELL）

1.534

1.690

2.340

src3
（+omp-parallel削減）

1.556

1.693

1.742

116

PCG法：計算時間：N=2003，48スレッド
N= 8,000,000
FortranC (clang)C (trad)src0（初期設定）

7.666

8.321

9.397

src1（First Touch）

6.952

5.102

8.008

src2（+ELL）

3.421

3.910

5.381

src3
（+omp-parallel削減）

3.440

3.920

3.824

117

PCG法：計算時間：N=2563，48スレッド
N= 16,777,216
FortranC (clang)C (trad)src0（初期設定）

34.30824.77225.547src1（First Touch）

32.20222.17223.814src2（+ELL）

8.916

10.761

14.566src3（+omp-parallel削減）

8.915

10.76410.415

118

src_f4 (1/5)
parallel computing by OpenMP
module solver_PCG
contains!C!C*** solve_PCG!C
subroutine solve_PCG&&( N, NPLU, indexLU, itemLU, D, B, X, AMAT, EPS, ITR, IER)
use omp_lib
implicit REAL*8 (A-H,O-Z)real(kind=8), dimension(N):: D
real(kind=8), dimension(N):: B
real(kind=8), dimension(N):: X
real(kind=8), dimension(NPLU) :: AMAT
integer, dimension(0:N)
integer, dimension(NPLU):: indexLU:: itemLU

real(kind=8), dimension(:,:), allocatable :: W
integer(kind=4), dimension(: ), allocatable :: SMPindex
integer, parameter :: R= 1
integer, parameter :: Z= 2
integer, parameter :: Q= 2
integer, parameter :: P= 3
integer, parameter :: DD= 4
real(kind=8), dimension(:) , allocatable :: W_RHO, W_C1, W_DNRM2

allocate (W(N+N2,4))!$omp parallel do private(i)do i= 1, N
X(i) = 0.d0
W(i,2)= 0.0D0
W(i,3)= 0.0D0W(i,DD)= 1.d0/D(i)enddo!$omp parallel do private(i)do i= N+1, N+N2
X(i) = 0.d0
W(i,2)= 0.0D0
W(i,3)= 0.0D0W(i,DD)= 1.d0/D(i)enddo!$omp parallel
PEsmpTOT= omp_get_num_threads()!$omp end parallel119

src_f4 (2/5)

PEsmpTOT:総スレッド数取得

allocate (SMPindex(0:PEsmpTOT))
SMPindex(0:PEsmpTOT):各スレッド受持要素番号
SMPindex(0)= 0
m = N/PEsmpTOT
nr= N - PEsmpTOT*m
do ip= 1, PEsmpTOT
SMPindex(ip)= m
if (ip.le.nr) SMPindex(ip)= m+1
enddo
do ip= 1, PEsmpTOT
SMPindex(ip)= SMPindex(ip)+SMPindex(ip-1)
enddo
allocate (W_RHO(PEsmpTOT), W_C1(PEsmpTOT), W_DNRM2(PEsmpTOT))

内積用

120

src_f4 (3/5)!$omp parallel do private(i,VAL,k)do i= 1, N
VAL= D(i)*X(i)do k= 1, 6
VAL= VAL + AMAT(k,i)*X(itemLU(k,i))
enddo
W(i,R)= B(i) - VAL
enddo
BNRM2= 0.0D0!$omp parallel do private(i) reduction(+:BNRM2)
do i= 1, N
BNRM2 = BNRM2 + B(i) **2
enddo

121

ITR= N
Stime= omp_get_wtime()!$omp parallel private(L,ip,ip1,ip2,i,k,VAL)!$omp&
private(RHO,BETA,RHO1,C1,ALPHA,DNRM2)
do L= 1, ITR
ip = omp_get_thread_num()+1
ip1= SMPindex(ip-1)+1
ip2= SMPindex(ip)!$omp simd
do i= ip1, ip2
W(i,Z)= W(i,R)*W(i,DD)
enddo
W_RHO(ip)= 0.0d0!$omp simd
do i= ip1, ip2
W_RHO(ip)= W_RHO(ip) + W(i,R)*W(i,Z)
enddo!$omp barrier
RHO= 0.d0!$omp simd
do i = 1, PEsmpTOT
RHO= RHO + W_RHO(i)enddo
if ( L.eq.1 ) then!$omp simd
do i= ip1, ip2
W(i,P)= W(i,Z)
enddo
else
BETA= RHO / RHO1!$omp simd
do i= ip1, ip2
W(i,P)= W(i,Z) + BETA*W(i,P)
enddo
endif
!$omp barrier

src_f4 (4/5)
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 i= ip1, ip2
VAL= D(i)*W(i,P)
do k= 1, 6
VAL= VAL + AMAT(k,i)*W(itemLU(k,i),P)
enddo
W(i,Q)= VAL
enddo
W_C1(ip)= 0.0d0!$omp simd
do i= ip1, ip2
W_C1(ip)= W_C1(ip) + W(i,P)*W(i,Q)
enddo!$omp barrier
C1= 0.d0!$omp simd
do i = 1, PEsmpTOT
C1= C1 + W_C1(i)enddo
ALPHA= RHO / C1!$omp simd
do i= ip1, ip2
X(i) = X(i)+ ALPHA * W(i,P)
W(i,R)= W(i,R) - ALPHA * W(i,Q)
enddo
W_DNRM2(ip)= 0.0d0!$omp simd
do i= ip1, ip2
W_DNRM2(ip)= W_DNRM2(ip) + W(i,R)**2
enddo!$omp barrier
DNRM2= 0.d0!$omp simd
do i = 1, PEsmpTOTDNRM2= DNRM2 + W_DNRM2(i)enddo
ERR = dsqrt(DNRM2/BNRM2)...

122

src_f4 (5/5)
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OMP-1123

src4 に対応したバージョン>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “INPUT.DAT”, “f4_48.sh”>>$ pjsub f4_48.sh124

Memory #5

err
f4_48.lst

CMG#0

CMG#2

13 14 15
12numactl

37 38 39
36numactl

18 19
16-C1712-23

-m 22
4 23
20 21
25 26 27
24numactl
30 31
28-C2924-35

リ
ン
グ
パ
ス

42 43
40-C4136-47

-m 46
6 47
44 45
49 50 51
48numactl
54 55
52-C5348-59

-m 34
5 35
32 33

-m 58
7 59
56 57

CMG#1

CMG#3

module load
fjexport OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand
numactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactlnumactl-l ./solf4-l ./solf4-l ./solf4-l ./solf4-l ./solf4
-C
12-59
-C
12-59
-C
12-59
-C
12-59
-C
12-59

-m
-m
-m
-m
-m

4-7
4-7
4-7
4-7
4-7

./solf4./solf4./solf4./solf4./solf4

Memory #7"f4_48”
rscgrp=tutorial-o
node=1
thread=48
elapse=00:15:00
gt00

Memory #6#!/bin/sh
#PJM
-N
#PJM
-L
#PJM
-L
#PJM
--omp
#PJM
-L
#PJM
-g
#PJM
-j
#PJM
-e
#PJM
-o

Memory #4

ｆ4_48.sh125

PCG法：計算時間：N=1283，48スレッド
N= 2,097,152
FortranC (clang)C (trad)src0（初期設定）

1.671

1.564

2.354

src1（First Touch）

1.480

1.122

1.720

src2（+ELL）

0.747

0.809

1.127

0.707

0.834

0.854

-

0.764

-

0.676

-

-

src3
（+omp-parallel削減）
src3b（+ clang loop unroll_count）
src4(+omp-do無し，reduction無し)

126

src_c3b (5/5)
src-c3#pragma omp for
for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[6*i+j]*W[P][itemLU[6*i+j]];}
W[Q][i] = VAL;}
src-c3b: clangのみ#pragma omp for#pragma clang loop unroll_count(8)for(i=0; i<N; i++) {VAL = D[i] * W[P][i];for(j=0; j<6; j++) {VAL += AMAT[6*i+j]*W[P][itemLU[6*i+j]];}
W[Q][i] = VAL;}

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

127

PCG法：計算時間：N=1603，48スレッド
N= 4,096,000
FortranC (clang)C (trad)src0（初期設定）

3.610

3.484

4.067

src1（First Touch）

2.993

2.228

3.425

src2（+ELL）

1.534

1.690

2.340

1.556

1.693

1.742

-

1.586

-

1.435

-

-

src3
（+omp-parallel削減）
src3b（+ clang loop unroll_count）
src4(+omp-do無し，reduction無し)

128

PCG法：計算時間：N=2003，48スレッド
N= 8,000,000
FortranC (clang)C (trad)src0（初期設定）

7.666

8.321

9.397

src1（First Touch）

6.952

5.102

8.008

src2（+ELL）

3.421

3.910

5.381

3.440

3.920

3.824

-

3.624

-

3.276

-

-

src3
（+omp-parallel削減）
src3b（+ clang loop unroll_count）
src4(+omp-do無し，reduction無し)

129

PCG法：計算時間：N=2563，48スレッド
N= 16,777,216
FortranC (clang)C (trad)src0（初期設定）

34.30824.77225.547src1（First Touch）

32.20222.17223.814src2（+ELL）

8.916

10.761

14.5668.915

10.76410.415-10.003

-

8.620

-

-

src3
（+omp-parallel削減）
src3b（+ clang loop unroll_count）
src4(+omp-do無し，reduction無し)

OMP-3130

まとめ• 有限体積法によるポアソン方程式ソルバー，PCG法による連立一次方程式求解⇒OpenMP並列化• Wisteria/BDEC-01（Odyssey）（A64FX）による実習– 複数CMGを使う– Cコンパイラの挙動⇒2021年11月時点より改善• 様々な最適化– First Touch Data Placement– ELL，OpenMPオーバーヘッド削減• 一日速習に丁度よい分量，であると考えている– アンケートで忌憚のないご意見をお願いいたします• データ依存性を含む場合（ICCG）は，3日間，2022年秋以降に開催予定

OMP-3131

性能評価：詳細プロファイラ• 測定部をプログラム内で指定（複数設定可）• 17回実行• 計算性能，消費電力等詳細なデータを得られる• Excel（マクロ）– https://www.dropbox.com/s/kat9ny5aoxp7cqm/cpu_pa_rep
ort.xlsm?dl=0

OMP-3

詳細プロファイラ利用（1/4）

測定部分をプログラム内で指定，特殊なコンパイルオプ
ション不要：solver_PCG.c/f#include "fj_tool/fapp.h"call fapp_start ("CG",1,0)
Stime = omp_get_wtime()

fapp_start ("CG",1,0);
Stime = omp_get_wtime();for(L=0; L<(*ITR); L++) {...
if(ERR < EPS) {*IER = 0;
goto N900;} else {RHO1 = RHO;}}*IER = 1;
N900:
Etime = omp_get_wtime();
fapp_stop ("CG",1,0);
return 0;}

do L= 1, ITR...
if (ERR .lt. EPS) then
IER = 0goto 900
else
RHO1 = RHO
endif
enddo
IER = 1900 continue
Etime= omp_get_wtime()call fapp_stop ("CG",1,0)
returnend

132

OMP-1

133

詳細プロファイラ利用（2/4）
Wisteria/BDEC-01上での実行>$ cd /work/gt00/t00XYZ/ompw>$ cd run<modify “fapp.sh”, “data.sh”>>$ pjsub fapp.sh(after finishing)>$ pjsub data.sh>$ ls pa*.csv
pa1.csv … pa17.csv

OMP-1

詳細プロファイラ利用（3/4）

134

Wisteria/BDEC-01上での実行スクリプト
各ディレクトリが既に生成されている場合は中身を空にしておくこと（上書きはしない）
fapp.shdata.sh#!/bin/sh#PJM -N "fapp"#PJM -L rscgrp=tutorial-o#PJM –L node=1#PJM –omp thread=48#PJM –L elapse=00:15:00#PJM –g gt00
#PJM –j#PJM –e err#PJM –o test.lst#!/bin/sh#PJM -N "data"#PJM -L rscgrp=lecture-o#PJM -L node=1#PJM --mpi proc=1#PJM -L elapse=00:15:00#PJM -g gt00
#PJM -j#PJM -e err#PJM -o data.lst

module load fj
export OMP_NUM_THREADS=48
export XOS_MMM_L_PAGING_POLICY=demand:demand:demand

module load fj
module load fjmpi

fapp -C -d ./repo01 -Hevent=pa1 ./solf2
fapp -C -d ./repo02 -Hevent=pa2 ./solf2
fapp -C -d ./repo03 -Hevent=pa3 ./solf2
fapp -C -d ./repo04 -Hevent=pa4 ./solf2
fapp -C -d ./repo05 -Hevent=pa5 ./solf2
fapp -C -d ./repo06 -Hevent=pa6 ./solf2
fapp -C -d ./repo07 -Hevent=pa7 ./solf2
fapp -C -d ./repo08 -Hevent=pa8 ./solf2
fapp -C -d ./repo09 -Hevent=pa9 ./solf2
fapp -C -d ./repo10 -Hevent=pa10 ./solf2
fapp -C -d ./repo11 -Hevent=pa11 ./solf2
fapp -C -d ./repo12 -Hevent=pa12 ./solf2
fapp -C -d ./repo13 -Hevent=pa13 ./solf2
fapp -C -d ./repo14 -Hevent=pa14 ./solf2
fapp -C -d ./repo15 -Hevent=pa15 ./solf2
fapp -C -d ./repo16 -Hevent=pa16 ./solf2
fapp -C -d ./repo17 -Hevent=pa17 ./solf2

fapp -A -d./repo01 -Icpupa,mpi -tcsv -o pa1.csv
fapp -A -d./repo02 -Icpupa,mpi -tcsv -o pa2.csv
fapp -A -d./repo03 -Icpupa,mpi -tcsv -o pa3.csv
fapp -A -d./repo04 -Icpupa,mpi -tcsv -o pa4.csv
fapp -A -d./repo05 -Icpupa,mpi -tcsv -o pa5.csv
fapp -A -d./repo06 -Icpupa,mpi -tcsv -o pa6.csv
fapp -A -d./repo07 -Icpupa,mpi -tcsv -o pa7.csv
fapp -A -d./repo08 -Icpupa,mpi -tcsv -o pa8.csv
fapp -A -d./repo09 -Icpupa,mpi -tcsv -o pa9.csv
fapp -A -d./repo10 -Icpupa,mpi -tcsv -o pa10.csv
fapp -A -d./repo11 -Icpupa,mpi -tcsv -o pa11.csv
fapp -A -d./repo12 -Icpupa,mpi -tcsv -o pa12.csv
fapp -A -d./repo13 -Icpupa,mpi -tcsv -o pa13.csv
fapp -A -d./repo14 -Icpupa,mpi -tcsv -o pa14.csv
fapp -A -d./repo15 -Icpupa,mpi -tcsv -o pa15.csv
fapp -A -d./repo16 -Icpupa,mpi -tcsv -o pa16.csv
fapp -A -d./repo17 -Icpupa,mpi -tcsv -o pa17.csv

OMP-1

135

詳細プロファイラ利用（4/4）
PC上での実行• pa*.csvをPC上へコピー>$ scp t00XYZ@wisteria.cc.u-tokyo.ac.jp:/work/gt00/t00XYZ/ompw/run/pa*.csv .

• 上記pa*.csvとExcelマクロファイルを同じディレクトリに置く• Excelマクロファイルをダブルクリック– 指示に従う– AllではなくCGを選択

136

PCG法：計算時間：N=1603，48スレッド
N= 4,096,000，ノード当たり換算，Fortran
最適化⇒メモリビジー⇒消費電力増加，Joule値は低下
計算時間(sec.)

対ピー
ク性能
比
(%)

SIMD
化率
(%)

MemoryThroughput
(%)

solf0

3.69

1.59

20.0

30.3

solf1
(First Touch)

2.98

1.97

28.8

37.5

solf2
(+ ELL)

1.58

3.73

49.7

70.0

solf3
(+ reduced“ompparallel”)
solf4（+ further
optimization)

1.58

1.45

3.72

4.06

48.4

55.1

69.8

75.8

Instruction

Power (W)Core
L2
Effective
Node
Memory81.6
3.39×1011 8.27×1010 10.9
112.
20.2
92.1
2.35×1011 5.33×1010 10.8
136.
33.3
104.
1.19×1011 4.17×1010 15.0
170.
51.7
Load/
Store

10

101.
14.9
51.0

167.

1.10×1011 3.78×1010

102.
15.8
56.5

174.

11

1.22×10

4.10×10

137

PCG法：性能評価：N=1603，48スレッド
Fortran，各CMG
src0:Initialsrc1:First Touch

Cycle Accounting execution time(s)src2: +ELL

Cycle Accounting execution time(s)4.0E+00

Other instruction commit

Cycle Accounting execution time(s)4.0E+00

4.0E+00

4 instruction commit
3 instruction commit

3.5E+00

3.5E+00

3.5E+00

3.0E+00

3.0E+00

3.0E+00

2.5E+00

2.5E+00

2.5E+00

2.0E+00

2.0E+00

2.0E+00

2 instruction commit
1 instruction commit
Barrier synchronization wait
Instruction fetch wait
Store port busy wait
Other wait
Branch instruction wait
Floating-point operation wait
Integer operation wait
Floating-point load L1D cache access wait (*)

1.5E+00

1.5E+00

1.5E+00

1.0E+00

1.0E+00

1.0E+00

Floating-point load L2 cache access wait
Integer load L1D cache access wait
Integer load L2 cache access wait
Floating-point load memory access wait
Integer load memory access wait

5.0E-01

5.0E-01

5.0E-01

0.0E+00

0.0E+00

0.0E+00

Prefetch port busy wait by software prefetch
Prefetch port busy wait by hardware prefetch

Process

Process

Process000■ Prefetch Port Busy Wait (H/W)■ Integer Load L1D Cache Access Wait
■ Floating-Point Operation Wait
■ ■ ■ Instruction Commit

Thread 11

Thread 9

Thread 10

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0

Thread 11

Thread 10

Thread 9

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0

Thread 11

Thread 10

Thread 9

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0---------------------------------------Floating-point busy rate execution time
Integer busy rate execution time
L1 busy rate execution time
L2 busy rate execution time
Memory busy rate execution time(*)Include wait time for integer L1D cache access

138

PCG法：性能評価：N=1603，48スレッド
Fortran，各CMG
src2:ELL

src3:omp-parallel
削減

Cycle Accounting execution time(s)src4: omp-do無し，reduction無し

Cycle Accounting execution time(s)Other instruction commit

Cycle Accounting execution time(s)1.8E+00

1.8E+00

1.8E+00

1.6E+00

1.6E+00

1.6E+00

4 instruction commit
3 instruction commit
2 instruction commit
1 instruction commit

1.4E+00

1.4E+00

1.4E+00

1.2E+00

1.2E+00

1.2E+00

1.0E+00

1.0E+00

1.0E+00

8.0E-01

8.0E-01

8.0E-01

6.0E-01

6.0E-01

6.0E-01

4.0E-01

4.0E-01

4.0E-01

2.0E-01

2.0E-01

2.0E-01

0.0E+00

0.0E+00

0.0E+00

Barrier synchronization wait
Instruction fetch wait
Store port busy wait
Other wait
Branch instruction wait
Floating-point operation wait
Integer operation wait
Floating-point load L1D cache access wait (*)
Floating-point load L2 cache access wait
Integer load L1D cache access wait
Integer load L2 cache access wait
Floating-point load memory access wait
Integer load memory access wait
Prefetch port busy wait by software prefetch
Prefetch port busy wait by hardware prefetch

Process0

Process0

Process0■ Prefetch Port Busy Wait (H/W)■ Floating-Point Load Memory Access Wait
■ Floating-Point Load L2D Cache Access Wait
■ Floating-Point Load L1D Cache Access Wait
■ ■ ■ Instruction Commit

Thread 11

Thread 9

Thread 10

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0

Thread 11

Thread 10

Thread 9

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0

Thread 11

Thread 10

Thread 9

Thread 8

Thread 7

Thread 6

Thread 5

Thread 4

Thread 3

Thread 2

Thread 1

Thread 0---------------------------------------Floating-point busy rate execution time
Integer busy rate execution time
L1 busy rate execution time
L2 busy rate execution time
Memory busy rate execution time(*)Include wait time for integer L1D cache access139

PCG法：計算時間：N=1603，48スレッド
N= 4,096,000，ノード当たり換算，C（clang）
最適化⇒メモリビジー⇒消費電力増加，Joule値は低下
Time(sec.)

Peak
Perf.
Ratio
(%)

SIMD
Ratio
(%)

MemoryThroughput
(%)

solc0

3.53

0.90

2.95

31.7

solc1
(First Touch)

2.42

1.34

4.82

46.2

solc2
(+ ELL)

1.73

1.88

8.21

63.8

solc3
(+ reduced“ompparallel”)
solc3b（+ clang loop
unroll_count）

1.75

1.62

1.86

2.00

8.51

10.4

62.9

67.8

Instruction

Power (W)Core
L2
Effective
Node
Memory103.
4.83×1011 1.65×1011 12.9
134.
17.9
107.
2.97×1011 1.10×1011 11.9
153.
33.8
104.
1.74×1011 8.92×1010 13.9
163.
44.8
Load/
Store

10

108.
14.2
45.0

167.

1.37×1011 8.22×1010

109.
14.6
44.4

168.

11

1.68×10

8.88×10

140

PCG法：計算時間：N=1603，48スレッド
N= 4,096,000，ノード当たり換算，C（trad）
最適化⇒メモリビジー⇒消費電力増加，Joule値は低下
Time(sec.)

Peak
Perf.
Ratio
(%)

SIMD
Ratio
(%)

MemoryThroughput(%)

solc0org

4.14

2.66

27.3

27.1

solc1org(First Touch)

3.52

3.13

37.4

31.8

solc2org(+ ELL)

2.34

4.70

58.1

47.1

solc3org(+ reduced“ompparallel”)

1.70

3.46

47.1

64.8

Instruction

Power (W)Core
L2
Effective
Node
Memory91.3
3.78×1011 7.40×1010 11.3
121.
18.4
90.
2.75×1011 6.14×1010 9.86
128.
28.4
95.1
1.63×1011 4.77×1010 11.7
144.
37.6
Load/
Store

11

1.17×10

10

4.14×10

96.7
13.9
48.7

159.