コードのチューニング
八木学（理化学研究所 計算科学研究機構）

謝辞
松本洋介氏（千葉大学）
KOBE HPC Spring School 20172017年3月13日 神戸大学計算科学教育センター

自己紹介
・専門は宇宙プラズマの磁気圏領域※太陽起源の超音速プラズマ流と惑星磁場の相互作用により生じる領域
・磁気流体（MHD）シミュレーションや、テスト
粒子計算を用いた水星磁気圏のモデリング・観測データのデータベース整備、及びデータ解析ソフトウェアの開発など
・Fortran, IDL, Python等を計算・解析に使用

磁気流体方程式
𝜕𝜌=− 𝐯∙𝛁 𝜌−𝜌 𝛁∙𝐯
𝜕𝑡𝜕𝐯
11= − 𝐯 ∙ 𝛁 𝐯 − 𝛁𝑃 + 𝐉 × 𝐁
𝜕𝑡
𝜌
𝜌
𝜕𝑃
= − 𝐯 ∙ 𝛁 𝑃 − 𝛾𝑃 𝛁 ∙ 𝐯𝜕𝑡𝜕𝐁= 𝛁 × (𝐯 × 𝐁)
𝜕𝑡
基本形は一般的な流体方程式とほぼ同じ
違いは磁場に関する項（ローレンツ力）が入る程度）

磁気流体方程式
𝜕𝜌=− 𝐯∙𝛁 𝜌−𝜌 𝛁∙𝐯
𝜕𝑡𝜕𝐯
11= − 𝐯 ∙ 𝛁 𝐯 − 𝛁𝑃 + 𝐉 × 𝐁𝜕𝑡𝜌𝜌𝜕𝑃= − 𝐯 ∙ 𝛁 𝑃 − 𝛾𝑃 𝛁 ∙ 𝐯𝜕𝑡
𝜕𝐁
Induction Equation= 𝛁 × (𝐯 × 𝐁)
𝜕𝑡

Lorentz force

基本形は一般的な流体方程式とほぼ同じ
違いは磁場に関する項（ローレンツ力）が入る程度）

磁気流体方程式
𝜕𝜌=− 𝐯∙𝛁 𝜌−𝜌 𝛁∙𝐯
𝜕𝑡𝜕𝐯
11= − 𝐯 ∙ 𝛁 𝐯 − 𝛁𝑃 + 𝐉 × 𝐁𝜕𝑡𝜌𝜌𝜕𝑃= − 𝐯 ∙ 𝛁 𝑃 − 𝛾𝑃 𝛁 ∙ 𝐯𝜕𝑡
𝜕𝐁
Induction Equation= 𝛁 × (𝐯 × 𝐁)
𝜕𝑡

Lorentz force

基本形は一般的な流体方程式とほぼ同じ
違いは磁場に関する項（ローレンツ力）が入る程度→今回は省略（つまり圧縮性の中性流体とほぼ同じ）磁気流体方程式
移流項

𝜕𝜌=− 𝐯∙𝛁 𝜌−𝜌 𝛁∙𝐯
𝜕𝑡𝜕𝐯
11= − 𝐯 ∙ 𝛁 𝐯 − 𝛁𝑃 + 𝐉 × 𝐁𝜕𝑡𝜌𝜌𝜕𝑃= − 𝐯 ∙ 𝛁 𝑃 − 𝛾𝑃 𝛁 ∙ 𝐯𝜕𝑡𝜕𝐁Induction Equation= 𝛁 × (𝐯 × 𝐁)
𝜕𝑡

Lorentz force

基本形は一般的な流体方程式とほぼ同じ
違いは磁場に関する項（ローレンツ力）が入る程度→今回は省略（つまり圧縮性の中性流体とほぼ同じ）

差分法
・もっともシンプル（原始的）かつ汎用性のあるスキーム
・故に大抵の時間発展方程式に適用可能𝜕𝑓
𝜕𝑓 𝜕2 𝑓
= 𝐹(𝑓, , 2 , … )𝜕𝑡𝜕𝑥 𝜕𝑥
空間方向：中心差分二次精度：

𝜕𝑓𝑓 −𝑓= 𝑖+1 𝑖−1
𝜕𝑥2∆𝑥𝜕𝑓−𝑓 +8𝑓 −8𝑓 +𝑓
四次精度： = 𝑖+2 𝑖+1 𝑖−1 𝑖−2𝜕𝑥12∆𝑥

時間積分：オイラー法：𝑓 𝑛+1 = 𝑓 𝑛 + 𝐹(𝑥)∆𝑡
ルンゲクッタ（４次）：
𝑘1 = 𝑆(𝑓 𝑛 , 𝑡𝑛 )𝑘1
1𝑛𝑘2 = 𝑆 𝑓 + ∆𝑡, 𝑡𝑛 + ∆𝑡
22𝑘2
1𝑛𝑘3 = 𝑆 𝑓 + ∆𝑡, 𝑡𝑛 + ∆𝑡
22𝑘4 = 𝑆 𝑓 𝑛 + 𝑘3 ∆𝑡, 𝑡𝑛 + ∆𝑡
1𝑓 𝑛+1 = 𝑓 𝑛 + 𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 ∆𝑡
6

セミラグランジュ法（移流）𝜕𝑓𝜕𝑓= −𝑣𝜕𝑡
𝜕𝑥∴ 𝑓 = 𝑒 −𝑖(𝑥−𝑣𝑡)

i-1

３次精度（３次関数による補間）F x = ax 3 + bx 2 + cx + d𝑓 𝑛+1 = 𝑎(−𝑣∆𝑡)3 +𝑏(−𝑣∆𝑡)2 +𝑐 −𝑣∆𝑡 + 𝑑𝑣 > 0の場合

ii-2

セミラグランジュ法：
・方程式から移流項のみを抽出して解く
・それ以外の部分は、普通の差分で処理・物理量をグリッド間で補間（風上型）
・速度の符号から風上を判定

𝑣∆𝑡i+1𝑣 < 0の場合

i-1ii+2𝑣∆𝑡i+1

チューニングの必要性について
・大規模シミュレーションを行う上で重要。
→数カ月間ジョブを流すようなケースでは研究
効率を大きく左右する→競争の激しい業界では、素早く計算結果を出す必要があるとはいえ・・・
・計算を１日早くするためのチューニングに１年かかってしまうようでは意味がない。
・バランスが重要チューニングの方針
・最近のコンパイラの最適化能力は高い→コンパイラが最適化しやすいコードを書く
・処理時間の最も重い箇所（ホットスポット）を重点的にチューニングする→プロファイラを使ってホットスポットを特定
・質の良い高速なライブラリは積極的に利用する→利用頻度の高い典型問題（フーリエ変換、固有値問題など）に関しては、開発も活発

チューニングの方針
・チューニングを行うと、コードの可読性が犠牲になる場合もあり、大規模なアルゴリズムの改良は難しくなる→アルゴリズムはチューニングする前に検討せよ→バグのチェックも同様
・１つの方針として、まずは可読性の高いコードを書き、アルゴリズムの検討やバグ出しを終えたところで最後にチューニングを行う。

コンパイラによる最適化
・コンパイル時にオプションを指定することで、コンパイラが自動的に最適化してくれる。

・インライン展開
・SIMD化・自動並列化・自動並列化は計算順序の入れ替えを行うため、稀に計算結果が変わってしまうことがある。
・明示的に最適化レベルを０にする（-O0）ことで、結果が変わっていないか確認できる。

コンパイラによる最適化
・ コンパイラによっては、強力なコンパイルオプションをまとめたオプションを用意していることもある--- Intel fortran (ifort) : -fast--- FUJITSU (京, π-Computer) : -Kfast
・ デフォルトで最適化オプションが指定されている場合も多い
・詳しくは各コンパイラのマニュアル参照（スパコンを使う場合、スパコンマニュアルを参照するとよい）

SIMD: Single Instruction Multiple DataＩｎｔｅｌ系：SSE, AVX・１つの命令を同時に複数のデータに対して行う演算、すなわち一種のベクトル計算である。
・ベクトル長は2~4程度と所謂ベクトル計算機に比べるとずっと短いが、それでもSIMD化できるのとできないのでは数倍速度が変わる。
・コンパイルオプションで有効化

gfortran -mavx test.f90
frtpx -Ksimd=1

プロファイラの使用
各サブルーチンの経過時間等を計測
Linuxでは gprof コマンドを利用可能
gfortran, Intel fortranの場合
gfortran -pg test.f90
or ifort -p test.f90./a.out
gprof ./a.out gmon.out

富士通システム（π-computer）
fipp -C -d <dir> ./a.out# 非MPI
fipp -C -d <dir> mpiexec ./a.out # MPIジョブ
fipppx -A -d <dir>

表示例：
Each sample counts as 0.01 seconds.
% cumulative self
self total
time seconds seconds calls us/call us/call name
各ルーチンの実行時間（割合）

呼び出し
回数

ルーチン名

メモリの階層構造
fast & small

CPU内部

・上に行くほど高速・低容量

レジスタ・階層により、アクセス速度が桁で変わる（が、容量も桁で落ちる）

L1

L2

RAM

slow & large

Disk

アクセス速度：L2キャッシュ > 100GB/sec
メインメモリ 10~ GB/sec
HDD 100MB~1GB/sec

割り算・べき乗
除算は加算や乗算に比べて計算コストが重いた
め、乗算の形で表記した方が速い
べき乗表記も可能であれば置き換える
do i=1,nx
f(i) = x(i)/cenddo

d = 1.0/c
do i=1,nx
f(i) = x(i) * denddo

do i=1,nx
f(i) = x(i)**2g(i) = x(i)**0.5
enddo

do i=1,nx
f(i) = x(i)*x(i)g(i) = sqrt((i))
enddo多項式の計算多項式は因数分解したほうが計算量が減る。
丸目誤差の出方が変わってしまうため、計算結果に注意する必要がある。
do i=1,nx
f(i) = a*x(i)*x(i)*x(i) + b*x(i)*x(i) + c*x(i) + d
enddo

do i=1,nx
f(i) = ((a*x(i) + b)*x(i) + c)*x(i) + d
enddo

ループ内の依存関係
ループ番号間に依存関係があるとSIMD化の妨げになり、実行速度が大きく損なわれる。工夫次第で依存を解消できるのであれば、書き換えたい。
●x(1) = dx
do i=2,nxx(i) = x(i-1) + dxenddo

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

実際には、こんな単純に行くことは少ないが・・・ループ内IF文の除去
計算ループ内のIF文は、コンパイラによる自動最適化（SIMD化）の妨げになるため、できる限り別の記述で置き換えるほうが好ましい。
●

min, max, sign等で代用できるケースあり●

do i=1,nx
if (x(i) > 0.0) thenf(i) = c * x(i)elsef(i) = -c * x(i)
endif
enddo

do i=1,nx
f(i) = sign(1.0,x(i))*c*x(i)enddo

ループ内IF文の除去
計算ループ内のIF文は、コンパイラによる自動最適化（SIMD化）の妨げになるため、できる限り別の記述で置き換えるほうが好ましい。
●

min, max, sign等で代用できるケースあり●

do i=1,nx
if (x(i) > 0.0) thenf(i) = c * x(i)elsef(i) = d * x(i)
endif
enddo

do i=1,nx
f(i) = 0.5 * x(i) &*(c+d+(c–d)*sign(1.0, x(i)))
enddo演習問題1-1
演習：1d_v101.f90 を改善せよ
1d_v101.f90 には簡単に修正できるポイントがいくつか存在している。
（ヒント）・ループ内のif文はsignを使えば消せる
・多項式は因数分解可能
・割り算は、掛け算の形式に変換
解答例：1d_v102.f90

補足：πコンピュータ上の操作
ファイルのコピー
ログイン（cygwinの場合）
ssh -X xxxx@pi.ircpi.kobe-u.ac.jp
演習用のディレクトリの作成[xxxx@pi ~]$ cd[xxxx@pi ~]$ mkdir 2017spring[xxxx@pi ~]$ cd 2017spring
演習用のファイルを自分のhomeディレクトリにコピー[xxxx@pi ~]$ cp /tmp/2017spring/code.tar ./[xxxx@pi ~]# tar xvf code.tar

補足：πコンピュータ上の操作"code"のディレクトリ構成
basic：
今から使うプログラム（非並列）openmp：
OpenMP用のプログラム
hybrid：ハイブリッド並列用のプログラム
python：
結果確認のためのプロット用[xxxx@pi 2017spring]$ cd code/basic[xxxx@pi basic]$ ls
1d_v101.f90 1d_v102.f90 1d_v201.f90.......

補足：πコンピュータ上の操作
コンパイル[xxxx@pi basic]$ frtpx -Kfast 1d_v101.f90
frtpx : Fortran90のコンパイルコマンド-Kfast : 自動最適化オプション詰め合わせ。
＊詳しくはπコンピュータマニュアル参照*C言語の場合は fccpx補足：πコンピュータ上の操作
ジョブスクリプトの投入[xxxx@pi basic]$ pjsub run.shrun.sh：スケジューラへの指令＋シェルスクリプト#!/bin/sh#PJM -L "node=1"#PJM -L "elapse=00:01:00"#PJM =L "rscgrp=small"#PJM -j#PJM -o "output.lst"./a.out

ノード数
実行時間
実行キュー
標準出力ファイル名
実行ファイル補足：πコンピュータ上の操作プログラム一覧：
1d_v101.f90 ： １次元移流方程式
1d_v102.f90 ： １次元移流方程式（改善案）
1d_v201.f90 ： １次元流体（オイラー法、空間２次）
1d_v202.f90 ： １次元流体（オイラー法、空間２次）
1d_v203.f90 ： １次元流体（オイラー法、空間４次）
1d_v204.f90 ： １次元流体（ルンゲクッタ、空間４次）
2d_v101.f90 ： ２次元流体（ルンゲクッタ、空間４次）

補足：πコンピュータ上の操作プログラム中のパラメータ（移流）：1d_v101.f90, 1d_v102.f90
nx : 空間グリッド数
delt : 時間刻み幅（Δt）
velc : 速度（v）
プログラム中のパラメータ（流体） ：1d_v20x.f90, 2d_v102.f90
ff(... , ieq), gf(... , ieq)
ieq=1: 密度（ρ）
ieq=2: 圧力（P）
ieq=3: Vx
ieq=4: Vy（２次元のみ）補足：πコンピュータ上の操作
データ出力先：1d_v101.f90, 1d_v102.f90

dat_1d_v11d_v201.f90, 1d_v202.f90
1d_v203.f90, 1d_v204.f90

dat_1d_v22d_v101.f90

dat_2d補足：πコンピュータ上の操作
データ出力の確認（Pythonプログラム）プログラムのあるディレクトリにパスを通す[xxxx@pi]$ emacs ~/.bashrc# .bashrc# Source global definitions.......
export PATH=$HOME/2017spring/code/python:$PATH#ファイルの最後にこの１行を追加

補足：πコンピュータ上の操作
データ出力の確認（Pythonプログラム）

１次元[xxxx@pi]$ cd ~/2017spring/code/basic/dat_1d_v1[xxxx@pi dat_1d_v1]$ plot1d.py 010000.dat
ファイル名

プロット結果補足：πコンピュータ上の操作
データ出力の確認（Pythonプログラム）

２次元[xxxx@pi]$ cd ~/2017spring/code/basic/dat_2d[xxxx@pi dat_2d]$ plot2d.py pre_001000.dat
ファイル名

プロット結果インライン展開
外部関数やサブルーチンはプログラムの可読性向上に役立つが、計算ループ内で繰り返し呼び出す場合は、呼び出しのオーバヘッドが大きい→内容をその場に展開する
do i=1,nx
f(i) = a + b * c(i)call myfunc(f(i))
enddosubroutine myfunc(g)g = g + 1.0
end subroutine

do i=1,nx
f(i) = a + b * c(i)f(i) = f(i) + 1.0
enddoインライン展開
・最近のコンパイラは、オプションを指定することでインライン展開もやってくれる。
コンパイラによるインライン展開の指定
gfortran -O3
ifort -ipo
frtpx -Kilfunc
・・・ただし、うまく機能しない場合もあるメモリの連続性
real(8) :: a(nx,ny)a(i-1,j-1)

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

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

a(2,j)

a[i][0]a[i][1]

内側の配列がメモリ空間上で連続する。

配列のメモリ空間上での配置(C言語)

double a[nx][ny]a[i-1][j-1]a[i-1][j]a[i-1][j+1]....

a[i-1][ny-1]

外側の配列がメモリ空間上で連続する。

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

for (i=0; i<nx; i++) {for (j=0; j<ny; j++) {a[i][j] = i + j}}nx 間隔でメモリアクセスが飛び飛びであり、メモリへの書き込みが非常に遅い

アドレスに連続アクセスするように修正

C言語の場合は、こちらが連続Byte/Flopの概念
１回の演算に対して必要なメモリアクセス量をByte/Flopで定義する
do i = 1, nx
a(i) = a(i) + b * c(i)enddo・この例では１回のループにつき、２回の浮動小数点演算に対して３回のデータロードが必要。
→倍精度とすると 3*8Byte / 2Flop = 12Byte/Flop・Byte/Flopは小さいほうが良いループアンローリング
Do ループを展開し、ループ内の実行文を増やす- ループオーバーヘッドの減少（ループ毎に発生するループの終了条件テストを削減）- レジスタブロッキングを行う
do i = 1, nx
a(i) = a(i) + b * c(i)enddo

do i = 1, nx, 4a(i) = a(i) + b * c(i)a(i+1) = a(i+1) + b * c(i+1)
a(i+2) = a(i+2) + b * c(i+2)
a(i+3) = a(i+3) + b * c(i+3)
enddoループアンローリング
注意点：- アンロールの段数に対して十分なレジスタが必要- ループ回数はアンロール段数で割り切れなければならない- コンパイラの最適化で自動的に行われることもある
do i = 1, nx
a(i) = a(i) + b * c(i)enddo

do i = 1, nx, 4a(i) = a(i) + b * c(i)a(i+1) = a(i+1) + b * c(i+1)
a(i+2) = a(i+2) + b * c(i+2)
a(i+3) = a(i+3) + b * c(i+3)
enddoループ分割・融合
２つのループを１つにまとめることで、オーバーヘッドを減らす（ループ融合）
do i = 1, nx
a(i) = b(i) * c(i)enddo
do i = 1, nx
d(i) = e(i) * f(i)enddo

do i = 1, nx
a(i) = b(i) * c(i)d(i) = e(i) * f(i)enddo

ループ内の依存関係によりSIMD化されない場合がある。
両者を分けることでSIMD化する（ループ分配）参考・スカラーチューニングとOpenMPによるコードの高速化 （松本洋介）
http://www.astro.phys.s.chiba-u.ac.jp/hpci/ss2013/presenymatsumoto.pdf
・チューニング技法入門 （青山幸也）
http://accc.riken.jp/wpcontent/uploads/2015/06/secure_4467_tuning-technique_main.pdf
・πコンピュータマニュアル
http://www.eccse.kobe-u.ac.jp/pi-computer/