スカラチューニングとOpenMPによるコードの高速化
松本洋介
千葉大学理学研究科
謝辞
C言語への対応：簑島敬（JAMSTEC）
宇宙磁気流体・プラズマシミュレーションサマースクール
2015年8月4日 千葉大学統合情報センター

内容
イントロダクション● スカラチューニング
● OpenMPによる並列化● 最近のHPC分野の動向● まとめ●

イントロダクション

命令実行の流れ

http://japan.intel.com/contents/museum/mpuworks/index.html

実行にかかる時間は主に、1.命令による実行（②、③）2.データの入出力（①、④）による。

メモリの階層構造
CPU
レジスタ
L1キャッシュL2キャッシュL3キャッシュDDR3, DDR4, HMCで接続

メインメモリregister : サイクル
● L1$: 32 kB, 数サイクル● L2$: 6MB, >180 GB/s
● Memory: 16GB, 64GB/s●「京」の場合同一チップ上

私のスパコン利用暦
2000
vpp800(Fujitsu)

2003 20042005

vpp5000(Fujitsu)

2009

2012

FX1(Fujitsu)
FX10(Fujitsu)
京(Fujitsu)
sx6(NEC)

sx9(NEC)
XT4(Cray)
XC30(Cray)2015私のスパコン利用暦
ベクトル計算機時代
2000
vpp800(Fujitsu)

2003 20042005

vpp5000(Fujitsu)

2009

2012

FX1(Fujitsu)
FX10(Fujitsu)
京(Fujitsu)
sx6(NEC)

sx9(NEC)
XT4(Cray)
XC30(Cray)2015私のスパコン利用暦
ベクトル計算機時代
2000
vpp800(Fujitsu)

2003 20042005

vpp5000(Fujitsu)

スカラ計算機時代
2009

2012

FX1(Fujitsu)
FX10(Fujitsu)
京(Fujitsu)
sx6(NEC)

sx9(NEC)
XT4(Cray)
XC30(Cray)2015スカラ／ベクトル？
●

スカラ
同じ車を何台も作る作業に例えると…

●ベクトル（データのパイプライン処理）
サイクル

２台目

３台目

do k=1,100
do j=1,100
do i=1,100a(i,j,k) = c*b(i,j,k)+d(i,j,k)
enddo
enddo
enddo

近年の演算処理の高速化のしくみ
サイクル時間（1/周波数）の短縮（→周波数～3GHzで頭打ち）● ベクトル化– パイプライン処理– SIMD（複数演算器）● メモリ構造の階層化● 並列化（マルチコア/MIMD）●

近年の演算処理の高速化のしくみ
サイクル時間（1/周波数）の短縮（→周波数～3GHzで頭打ち）● ベクトル化– パイプライン処理– SIMD（複数演算器）● メモリ構造の階層化● 並列化（マルチコア/MIMD）●

近年の演算処理の高速化のしくみ
サイクル時間（1/周波数）の短縮（→周波数～3GHzで頭打ち）● ベクトル化– パイプライン処理
ユーザから見たら同じ– SIMD（複数演算器）● メモリ構造の階層化● 並列化（マルチコア/MIMD）●

近年の演算処理の高速化のしくみ
サイクル時間（1/周波数）の短縮（→周波数～3GHzで頭打ち）● ベクトル化– パイプライン処理
ユーザから見たら同じ– SIMD（複数演算器）● メモリ構造の階層化● 並列化（マルチコア/MIMD）●

近年の計算機では、SIMD化、キャッシュヒット率の向上、マルチコアによる並列化が高速化のポイント

スカラチューニング

対象
宇宙磁気流体プラズマシミュレーションにかかわること● すなわち、– 差分法：磁気流体（MHD）・ブラソフシミュレーション– 粒子法：電磁粒子（PIC）シミュレーション● 行列の演算（例：LU分解など）は対象外● Fortran, C●注意

一般に、チューニングすると可読性が損なわれます。まずは読みやすいコードを書き、充分テストしてバグを除いてからチューニングを行いましょう。

チューニングが必要？
無理にしなくて良いです。（好きでもしんどい）● 最先端（大規模）シミュレーション研究では必須。なぜなら、、、● 1ランで数週間→2倍の速度向上で10日単位の短縮● 「京」などの大規模計算申請書類では、実行効率・並列化率などの情報が求められる。
● 実行効率15%以上あれば、計算機資源の獲得において、他分野との競争力になる。
●

チューニングの手順
編集
コンパイル

コンパイルリスト

最適化されているか？
Yes

実行プロファイラ

一番重いルーチンは？

Noコンパイルリスト●●

最適化情報の詳細を出力– インライン展開等の最適化– SIMD化– 並列化
コンパイルオプション– gcc/gfortran: N/A– icc/ifort: -opt-report, -vec-report, -parreportプロファイラの利用
各サブルーチンの経過時間を計測● ホットスポット（一番処理が重いサブルーチン）から最適化● 商用コンパイラ（intel, PGI, スパコン等）では、詳細情報（キャッシュミス率、FLOPS）が得られる● GNUでは、gprof● gprofの使い方●––––gfortran (gcc) -pg test.f90
ifort (icc) -p test.f90./a.out
gprof ./a.out gmon.out> output.txt

output.txt

スカラチューニングのポイント
コンパイラ（＝人）にやさしいプログラム構造– ループ内で分岐は使わない（if文の代わりにmin, max,
sign で、 goto文は不可）– ループ内処理を単純にする（SIMD化促進）– 外部関数のインライン展開● データの局所化を高める– 繰り返し使用するデータはなるべくひとまとめにして、キャッシュに乗るようにする。
– 一時変数の再利用– 連続アクセス– ポインタは使わない（Fortran）●

基本的なtips
割り算を掛け算に– a(i) = b(i)/c → c = 1.0/c ; a(i) = b(i)*c● べき乗表記はなるべく使わない– a(i) = b(i)**2 → a(i) = b(i)*b(i)– a(i) = b(i)**0.5 → a(i) = sqrt(b(i))● 因数分解をして演算数を削減– y=a*x*x*x*x+b*x*x*x+c*x*x+d*x → y=x*(d+x*(c+x*(b+x*(a))))– 演算回数13 → 7● 一時変数は出来る限り再利用（レジスタの節約）●

分岐処理の回避1
Fortran:do i=1,nx
if(a /= 0.0)thenb(i) = c(i)/aelseb(i) = c(i)
endif
enddo

if(a == 0.0) a=1.0
a = 1.0/a
do i=1,nxb(i) = c(i)*aenddo

分岐処理の回避2（minmod関数）
Fortran:do i=1,nx
if(a(i)*b(i) < 0.0)thenc(i) = 0elsec(i) = sign(1.0,a(i)) &*min(abs(a(i)),abs(b(i)))
endif
enddo

do i=1,nxc(i) = sign(1.0,a(i))&*max(0.0,&min(abs(a(i)),&sign(1.0,a(i))*b(i)) &)enddo

SIMD: Single Instruction Multiple Data
ユーザレベルではベクトル化と同じ。ただし、ベクトル長は2～4と、ベクトル計算機のそれ（256）に比べてずっと短い。
● 最内側ループに対してベクトル化● 最近ではループ内にIF文が入っていてもSIMD化してくれる（マスク付きSIMD化）。真率が高ければ効果的。
● コンパイルオプションで最適化– gcc/gfortran: -m{avx, sse4}– icc/ifort: -x{avx,sse4}●

from K computer manualfrom K computer manual

SIMD化の阻害例
SIMD化されない

書き方の工夫SIMD化される

例１: ループ番号間に依存性がある場合
a(1) = dx
do i=2,nxa(i) = a(i-1)+dx
enddo

do i=1,nxa(i) = i*dx
enddo

例２: ループ番号によって処理が異なる場合
do i=1,nx
if(a(i) < 0)thenb(i-1) = c*a(i)elseb(i+1) = c*a(i)
endif
enddo

do i=1,nx
w1 = 0.5*(1.0-sign(1.0,a(i)))
w2 = 0.5*(1.0+sign(1.0,a(i)))
b(i-1) = c*w1*a(i)b(i+1) = c*w2*a(i)enddo配列の宣言とメモリ空間1 (Fortran)
dimension a(nx,ny)

配列aのメモリ空間上での配置は、
a(i-1,j-1) a(i,j-1) a(i+1,j-1)…...

a(i-1,j)nxa(i,j)

a(i+1,j)…...

a(i-1,j+1) a(i,j+1) a(i+1,j+1)nx

do i=1,nx
do j=1,nya(i,j) = i+j
enddo
enddo

do j=1,ny
do i=1,nxa(i,j) = i+j
enddo
enddo

間隔nxで飛び飛びにアドレスにアクセスすることになるので、メモリへの書き込みが非常に遅いアドレスに連続アクセスするので、メモリへの書き込みが速い配列の宣言とメモリ空間1 (C/C++)
float a[ny][nx];配列aのメモリ空間上での配置は、
a[j-1][i-1]a[j-1][i] a[j-1][i+1]…...

a[j][i-1]nx

a[j][i]a[j][i+1]…...

a[j+1][i-1] a[j+1][i] a[j+1][i+1]nxfor(i=0;i<nx;i++){for(j=0;j<ny;j++){a[j][i] = i+j;}}

for(j=0;j<ny;j++){for(i=0;i<nx;i++){a[j][i] = i+j;}}

間隔nxで飛び飛びにアドレスにアクセスすることになるので、メモリへの書き込みが非常に遅い

アドレスに連続アクセスするので、メモリへの書き込みが速い配列の宣言とメモリ空間2連続の式：rn+1 = f(rn,vxn,vyn)
次のステップに進むためには、自分自身（r）の他に速度場が必要
dimension rho(nx,ny), vx(nx,ny), vy(nx,ny), ...と変数を個別に用意する代わりに、dimension f(8,nx,ny) ! 1:rho, 2:p, 3-5:v, 6-8:Bのように、一つの変数にまとめて配列を用意する。このようにすると、必要となる各物理量がメモリ空間上の近い位置に配置される（キャッシュラインにのりやすい）。→システム方程式を解くための工夫

配列の宣言とメモリ空間2（続き）
TypeA: f(nx,nx,nz,8)
TypeB: f(8,nx,nx,nz)

近年のキャッシュ重視型のスパコンにおいて効果的
Fukazawa et al., IEEE Trans. Plasma Sci., 2010.

配列の宣言とメモリ空間3（C言語）
C言語で静的に配列を宣言する場合は、
float a[ny][nx];とするが、領域分割の並列計算では動的に（mallocで）配列を確保する場合が多く、上記の宣言では難しい。２次元配列を１次元配列として宣言する方が、メモリ空間上で連続的に領域を確保できる。
double *a;a=(double*)malloc(sizeof(double)*nx*ny);for (j=0;j<ny;j++){for (i=0;i<nx;i++){a[nx*j+i] = i+j;}}キャッシュと配列ブロック
２次キャッシュ容量～数MB– 倍精度で100x100x100グリッド分程度– MHD計算の1ノードあたりの配列数としては、まだちょっと足りない– PIC計算では、セル内の粒子が必要なグリッド上の場のデータがキャッシュに乗らない● 配列ブロック– 必要なデータだけをキャッシュに収まる程度の別の小さな配列に事前に格納– PIC法で(i,j,k)セルに属する粒子が必要な場の情報をあらかじめパック●

tmp(1:6,-1:1,-1:1,-1:1) = f(1:6,i-1:i+1,j-1:j+1,k-1:k+1)

インライン展開●

外部（ユーザー定義）関数はプログラムの可読性向上に一役。しかし、、
do i=1,nxa(i) = myfunc(b(i))
enddoのように、ループ内で繰り返し呼び出す場合、呼び出しのオーバーヘッドが大きい。関数内の手続きが短い場合は、内容をその場所に展開する→インライン展開● コンパイル時に指定（同一ファイル内に定義される関数）– gcc/gfortran: -O3 もしくは -finline-functions– icc: -O{2,3}, ifort: -finline● コンパイル時に指定（別ファイル内に定義される関数）– icc/ifort: -fast もしくは -ipo

OpenMPによるコードの並列化

アムダールの法則
全処理=1
並列数=1

逐次処理

並列化可能部分 (p)

(1-p)

並列数=2

p/2

...

…

並列数=n

性能向上率 =p/n1p(1− p)+n

少なくも並列化率p>0.95である必要ありスレッド並列とプロセス並列
スレッド並列プロセス並列

プロセス
スレッド1

スレッド2

スレッド3

プロセス1プロセス2

プロセス3

メモリ空間

メモリ空間

メモリ空間

プロセス4

プロセス5

プロセス6

メモリ空間

メモリ空間

メモリ空間

スレッドn

プライベート変 プライベート変 プライベート変 プライベート変数数数数

グローバル変数
メモリ空間

ハイブリッド並列
プロセス1-3

プロセス4-6

プロセス7
thrd1 thrd2 thrd3 thrd4

プロセス7-9

メモリ空間

この例では全72並列● プロセス間はMPIによる通信● 各プロセスに4スレッド● スレッド数分プロセス数を削減● MPIによる通信／同期待ちのオーバヘッドを軽減● 出力ファイル数の削減●

プロセス10-12

プロセス13-15

プロセス16-18

スレッド並列計算を行うためのAPI● コンパイルオプションで有効– gcc/gfortran: -fopenmp– icc/ifort: -openmp● プログラムに指示行を挿入（オプション無効時はコメント行と見なされる（C言語は警告される場合も）● 自動並列化に比べて柔軟に最適化が可能● 標準規格なため、マシン／コンパイラに依らずポータブル● 2015年8月現在、OpenMP 4.1。SIMD化の指示行、アクセラレータ（後述）への対応● http://www.openmp.org●

スレッド数の設定
基本的にはシェルの環境変数 $OMP_NUM_THREADS
でスレッド数を指定する– setenv OMP_NUM_THREADS 8– 指定しなければ、システムの全コア数● プログラム内部で関数で設定（omp_lib/omp.hをインクルードする必要あり）●

Fortran:!$use omp_lib
integer, parameter :: nthrd =8
call
omp_set_num_threads(nthrd)

C:#include <omp.h>
int nthrd=8;
omp_set_num_threads(nthrd);

全体の流れ：fork-join モデル
Fortran:
逐次処理
fork

並列処理

C:
#include <stdio.h>#include <omp.h>
int main(void){puts(“serial region”);#pragma omp parallel{…
…
…
puts(“parallel region”);…
…
…
}
puts(“serial region”);

join
逐次処理

return 0;}

ループの並列化*$OMP_SCHEDULE／SCHEDULE句で分担
方法変更可

i=1-100を各スレッドが!$OMP PARALLEL DO
均等に分担#pragma omp parallel for
do i=1,100
for (i=0;i<100;i++){b(i) = c*a(i)b[i]=c*a[i];enddo!$OMP END PARALLEL DOスレッドの立ち上げは}
なるべくまとめてcall mysub(b)mysub(b);#pragma omp parallel{#pragma omp forfor (i=0;i<100;i++){d[i]=c*b[i];}!$OMP PARALLEL!$OMP DO
do i=1,100d(i) = c*b(i)enddo!$OMP END DO
pragma omp for の#pragma omp for!$OMP DO
直後のforループが並列
for (i=0;i<100;i++){do i=1,100e[i]=c*d[i];処理される。間に”{”をe(i) = c*d(i)}
入れてはならない
enddo}!$OMP END DO!$OMP END PARALLEL

多重ループの並列化
do j=1,100!$OMP PARALLEL DO
do i=1,100b(i,j) = c*a(i,j)
enddo!$OMP END PARALLEL DO
enddo!$OMP PARALLEL DO &!$OMP PRIVATE(i)do j=1,100
do i=1,100b(i,j) = c*a(i,j)
enddo
enddo!$OMP END PARALLEL DO

スレッドの立ち上げが100回も行われ、オーバーヘッドが大きい

for (j=0;j<100;j++){#pragma omp parallel forfor (i=0;i<100;i++){b[j][i]=c*a[j][i];}}最外ループを並列化#pragma omp parallel
内側ループのカウンタ{
変数 i はプライベート#pragma omp for private(i)宣言が必要。
for (j=0;j<100;j++){for (i=0;i<100;i++){b[j][i]=c*a[j][i];}}}多重ループの並列化（続き）●●

多重ループでは最外ループを並列化するのが基本。ループの内側に指示行を入れると、外側ループの回転数分スレッドのfork/joinが行われ、スレッド立ち上げのオーバーヘッドが大きくなる。
内側にあるループのカウンタ変数（i, j, ..）はスレッド固有の変数とする必要があるため、PRIVATE宣言をする。そうしな
いと、スレッド間で上書きしてしまう。

グローバル／プライベート変数!$OMP PARALLEL DO
do i=1,100
tmp = myfunc(i)a(i) = tmp
enddo!$OMP END PARALLEL DO

スレッド間でtmpを上書きしまうので正しい結果が得られな
い#pragma omp parallel forfor (i=0;i<100;i++){tmp=myfunc(i);a[i]=tmp;}

Cの場合はループ内で変数宣言すれば問題なし。

!$OMP PARALLEL DO &!$OMP PRIVATE(tmp)
do i=1,100
tmp = myfunc(i)a(i) = tmp
enddo!$OMP END PARALLEL DO#pragma omp parallel{#pragma omp for private(tmp)for (i=0;i<100;i++){tmp=myfunc(i);a[i]=tmp;}}#pragma omp parallel forfor (i=0;i<100;i++){double tmp;
tmp=myfunc(i);a[i]=tmp;}ループ内変数の演算 (REDUCTION)
sum = 0.0!$OMP PARALLEL DO &!$OMP REDUCTION(+:sum)
do i=1,10
sum = sum+i
enddo!$OMP END PARALLEL DO

sum=1.0;#pragma omp parallel for reduction(+:sum)for (i=0;i<10;i++){sum+=i;}

総和（+）以外には、最大（max）、最小（min）が実用上使われる。

単スレッド処理 (SINGLE)!$OMP PARALLEL!$OMP DOdo i=1,100b(i) = c*a(i)enddo!$OMP END DO!$OMP SINGLE
call output(b)!$OMP END SINGLE!$OMP DOdo i=1,100d(i) = c*b(i)enddo!$OMP END DO!$OMP END PARALLEL

スレッドの立ち上げ #pragma omp parallelを最初に一回だけ {#pragma omp forfor (i=0;i<100;i++){b[i]=c*a[i];}#pragma omp single
途中で逐次処理が入る {場合はSINGLEで対処output(b);}#pragma omp forfor (i=0;i<100;i++){d[i]=c*b[i];}}スレッドの立ち上げ回数はなるべく少なく。データ入出力など、途中で逐次処理が必要な場合に使う。

バリア同期の回避 (NOWAIT)!$OMP PARALLEL!$OMP DOdo i=1,100b(i) = c*a(i)enddo!$OMP END DO NOWAIT!$OMP DOdo i=1,100d(i) = c*b(i)enddo!$OMP END DO!$OMP DO
do i=1,200e(i) = c*d(i)enddo!$OMP END DO NOWAIT!$OMP END PARALLEL

ループの終わりで暗黙に行われるスレッド間の同期待ちをNOWAITで回避#pragma omp parallel{#pragma omp for nowaitfor (i=0;i<100;i++){b[i]=c*a[i];}#pragma omp forfor (i=0;i<100;i++){d[i]=c*b[i];次のループではスレッドに対する変数 d の割り当 }て範囲が変わるので、同期が必要（注意）#pragma omp for nowaitfor (i=0;i<100;i+=2){
e[i]=c*d[i];}}

スレッド数が大きい場合に高速化に寄与する

OpenMP実装上の注意点
ユーザが並列処理箇所を明示するため、並列計算に伴う問題発生はプログラマが責任を負う（自動並列化との違い）。
● 並列処理してはいけない箇所でも、明示したら並列化されてしまう● スレッド内でグローバル/プライベート変数を間違えると結果が不定● NOWAITで必要な同期を忘れると結果が不定● 同じプログラムを数回は実行して、結果が変わらないことの確認が必要● 実装は簡単だけど、デバッグに注意が必要●

最近のHPC分野の動向

TOP500 (2015年6月現在）“ペタFLOPS・メガW時代”

http://www.top500.org/lists/2015/06/

TOP500 (2015年6月現在）“ペタFLOPS・メガW時代”

http://www.top500.org/lists/2015/06/

GREEN500（性能／消費電力）
PEZY（日本のベンチャー企業）が１－３位独占！

http://www.green500.org/lists/green201506
先週土曜日発表があったばかり

GPGPU vs. MIC vs. PEZY-SC
NVIDIA TESLA– ゲーム用途のGPUをHPCに応用（GPGPU)– CUDA/OpenACCによるプログラミング– SIMD– PGI Fortranでも可能（NVIDIAが買収）• Intel Xeon Phi– x86互換のコプロセッサ（~ 60 core）– 既存のコードから容易に拡張可能– SIMD/MIMD– 最新のランキングでは、TOP20に残らず• PEZY-SC– 日本のベンチャー企業が設計– メニーコアプロセッサ（1024個）– MIMD– C/C++しかコンパイラがないようです●エクサFLOPS・メガW時代へ
電力消費量はこれ以上増やせないので、これから専用CPUと組み合わせたスパコンが国内でも増えてくる● 汎用／専用CPU構成のヘテロジニアスなシステムへ● →ユーザのプログラム負担が増える可能性● シミュレーション研究者の宿命だが、5-10年くらいの周期でス
パコンシステムのトレンドに振り回される– ベクトル vs. スーパースカラ– MPI vs. HPF (High Performance Fortran)– 私は2009年に手持ちのコード（MHD/PIC）をスクラッチから再コーディング● スパコン情勢に注意しつつ、研究を進めましょう●まとめ●スカラチューニング
高速化のためのCPUの機能（SIMD）をいかに使い倒すか– キャッシュチューニング
–●OpenMPによるスレッド並列化
指示行を最外ループの手前にいれるだけ（簡単！）– スレッド並列化によりプロセス数を減らし、プロセス間通信のオーバーヘッドを軽減：ハイブリッド並列化
–

●

今後の展望
次世代のスパコンでは電力消費量問題が顕在化– 汎用／専用CPUで構成されるヘテロジニアスシステムに– →ハイブリッド並列化はますます必須–参考資料
プロセッサを支える技術、Hisa Ando著、技術評論社● 各スパコンマニュアル●

http://www.nag-j.co.jp/openMP/index.htm● http://accc.riken.jp/hpc/training/●