先週、浮動小数点数を足していくプログラムについて、Juliaで書いた場合とC言語で書いた場合とで速度を比較しました(浮動小数点数を足すだけのプログラムがC言語よりJuliaのほうが速い、について - よーる)。
この記事について、黒木玄さん(黒木玄 Gen Kuroki (@genkuroki) | Twitter)から「Juliaでは単純な和には普通にsimdを使うので、JuliaとCの両方でsimdを使って比較した方が良かったと思います。」というコメントをいただきました。
先週の記事で比較に用いたプログラムには並列性がないため、SIMD (single instruction multiple data) 命令を使った高速化は不可能なはずですが、Juliaでsimdを使うと10倍高速化するそうです。
一体どういうことなのだろうと思ってJuliaのマニュアル(Base.SimdLoop.@simd (macro))を見てみると、@simd
というアノテーションは浮動小数点数演算の順序を変える効果があるそうです。
浮動小数点数演算は一般に結合則が成り立ちません。
したがって、浮動小数点演算の順序を変えてしまうと、プログラムの意味が変わってしまいます(一般に計算結果も変わってしまいます)。
先週の記事で比較したプログラムはあくまで「浮動小数点数を前から順番に足していく」プログラムであって、「いい感じに浮動小数点数の和を求める」プログラムではありません。
この辺は、多くの読者が勘違いされているところだと感じました。
とはいえ、「浮動小数点数を前から順番に足していった値を正確に求めたい」というワークロードが常識的ではないのも事実です。
そこで、「浮動小数点数を順番を気にせずに足した値を求める」というワークロードで、JuliaとC言語の速度比較をしてみます。
ちょうど黒木玄さんがJulia側のソースコードを示してくださいました(simd example · GitHub)ので、それを使うことにします*1。
ただし、実行時間が短すぎて正確な評価ができなかったので、外側ループの回数を100倍にしました。
C言語側のソースコードは、私が素直に移植したものになります。
function test_inbounds_simd()
a = fill(1e-10, 1000)
s = 0.0
for j = 1:1000000000
@inbounds @simd for i = 1:1000
s += a[i]
end
end
s
end
test_inbounds_simd()
double test() {
unsigned long long i, j;
double a[1000];
double s = 0.0;
for( i = 0; i < 1000; ++i ) { a[i] = 1e-10; }
for( j = 0; j < 1000000000; ++j ) {
for( i = 0; i < 1000; ++i ) {
s += a[i];
}
}
return s;
}
int main() {
return test();
}
使用したハードウェア
- Intel(R) Xeon(R) Platinum 8280 CPU @ 2.70GHz
- 実際には3.6GHzくらいで動いているっぽい(?)
- AVX-512使用可能(一命令で倍精度浮動小数点数8個を同時に演算可能)
- FMA実行ユニット2基(1cycleに浮動小数点数加算命令が2回実行できる)
- 浮動小数点数加算命令のレイテンシは4cycle
使用したソフトウェア
julia version 1.5.2
(内部で使用しているLLVMのバージョンはLLVM version 9.0.1jl
)
clang version 10.0.0
icc (ICC) 19.0.5.281 20190815
perf
gdb
理論限界
ベンチマークプログラムは1000G個の浮動小数点数を加算するものです。
実験に使用したプロセッサは1サイクルに16個の浮動小数点数を加算することができますので、理論限界は62.5Gサイクルです。
3.6GHzで動いているとすると17秒程度かかることになりますが、周波数が安定しないので実時間は参考程度ということにします。
Juliaの結果
- 199G cycles
- 278G insns
- 43G branches
- 55秒くらい
JITコンパイルされた機械語をgdb
で確認してみると、以下のようになっていました。
0x00007fffc565f0d0: cmp rax,0x3b9aca00
0x00007fffc565f0d6: lea rax,[rax+0x1]
0x00007fffc565f0da: je 0x7fffc565f187
0x00007fffc565f0e0: vmovq xmm0,xmm0
0x00007fffc565f0e4: vxorpd xmm1,xmm1,xmm1
0x00007fffc565f0e8: mov esi,0x18
0x00007fffc565f0ed: vxorpd xmm2,xmm2,xmm2
0x00007fffc565f0f1: vxorpd xmm3,xmm3,xmm3
0x00007fffc565f0f5: nop WORD PTR cs:[rax+rax*1+0x0]
0x00007fffc565f0ff: nop
0x00007fffc565f100: vaddpd zmm0,zmm0,ZMMWORD PTR [rcx+rsi*8-0xc0]
0x00007fffc565f108: vaddpd zmm1,zmm1,ZMMWORD PTR [rcx+rsi*8-0x80]
0x00007fffc565f110: vaddpd zmm2,zmm2,ZMMWORD PTR [rcx+rsi*8-0x40]
0x00007fffc565f118: vaddpd zmm3,zmm3,ZMMWORD PTR [rcx+rsi*8]
0x00007fffc565f11f: add rsi,0x20
0x00007fffc565f123: cmp rsi,0x3f8
0x00007fffc565f12a: jne 0x7fffc565f100
0x00007fffc565f12c: vaddpd zmm0,zmm1,zmm0
0x00007fffc565f132: vaddpd zmm0,zmm2,zmm0
0x00007fffc565f138: vaddpd zmm0,zmm3,zmm0
0x00007fffc565f13e: vextractf64x4 ymm1,zmm0,0x1
0x00007fffc565f145: vaddpd zmm0,zmm0,zmm1
0x00007fffc565f14b: vextractf128 xmm1,ymm0,0x1
0x00007fffc565f151: vaddpd xmm0,xmm0,xmm1
0x00007fffc565f155: vpermilpd xmm1,xmm0,0x1
0x00007fffc565f15b: vaddsd xmm0,xmm0,xmm1
0x00007fffc565f15f: test dl,dl
0x00007fffc565f161: jne 0x7fffc565f0d0
0x00007fffc565f167: xor esi,esi
0x00007fffc565f169: nop DWORD PTR [rax+0x0]
0x00007fffc565f170: vaddsd xmm0,xmm0,QWORD PTR [rcx+rsi*8+0x1f00]
0x00007fffc565f179: inc rsi
0x00007fffc565f17c: cmp rsi,0x8
0x00007fffc565f180: jne 0x7fffc565f170
0x00007fffc565f182: jmp 0x7fffc565f0d0
この機械語コードから予測すると、1000個の浮動小数点数を足すのに
- 192cycle(zmmレジスタを使った主ループが4cycle×31回、マージしていく部分が36cycle、端数処理ループが4cycle×8回)
- 273命令(準備に11命令、主ループが7命令×31回、マージしていく部分が13命令、端数処理ループが4命令×8回)
- 42分岐命令
が必要です。
これはperf
コマンドの結果と整合します。
理論限界と大きく差がついている理由は、以下の通りです。
- 主ループのアンローリングが不足している
- 4cycleレイテンシを隠蔽しつつFMAユニット2基をフル稼働させるためには、依存関係にない命令が最低でも8つ必要
- 端数処理が直列化されている
コンパイルオプションは、-O2 -mavx512f -funsafe-math-optimizations
です*2。
C言語の範囲では「浮動小数点数演算の順序を自由に変更してかまわない」という記述をすることができないので、-funsafe-math-optimizations
をつけます*3*4。
- 194G cycles
- 213G insns
- 32G branches
- 55秒くらい
Juliaで書いた場合と同じくらいの速度になりました。それもそのはずで、JuliaのJITエンジンはLLVM、つまりclang
のバックエンドと同じです。
というわけでコード生成結果や実行速度が大きく変わるはずはないのです。
命令数が少なくなっているのは、JITコンパイルより時間をかけてコンパイルすることできれいなコードを生成できたことによるものと思われます。
.LBB0_5:
vaddpd zmm0, zmm2, zmm0
vaddpd zmm0, zmm3, zmm0
vaddpd zmm0, zmm4, zmm0
vextractf64x4 ymm2, zmm0, 1
vaddpd zmm0, zmm0, zmm2
vextractf128 xmm2, ymm0, 1
vaddpd xmm0, xmm0, xmm2
vpermilpd xmm2, xmm0, 1
vaddsd xmm0, xmm0, xmm2
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, xmm1
add rax, 1
cmp rax, 1000000000
je .LBB0_6
.LBB0_3:
vmovq xmm0, xmm0
vxorpd xmm2, xmm2, xmm2
mov ecx, 56
vxorpd xmm3, xmm3, xmm3
vxorpd xmm4, xmm4, xmm4
.p2align 4, 0x90
.LBB0_4:
vaddpd zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 576]
vaddpd zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 512]
vaddpd zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 448]
vaddpd zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 384]
cmp rcx, 1016
je .LBB0_5
vaddpd zmm0, zmm0, zmmword ptr [rsp + 8*rcx - 320]
vaddpd zmm2, zmm2, zmmword ptr [rsp + 8*rcx - 256]
vaddpd zmm3, zmm3, zmmword ptr [rsp + 8*rcx - 192]
vaddpd zmm4, zmm4, zmmword ptr [rsp + 8*rcx - 128]
add rcx, 64
jmp .LBB0_4
zmmレジスタを使った主ループが二倍アンローリングされたほか、端数ループが八倍(完全)アンローリングされています。これらが命令数削減に役立っているようです。
この機械語コードから予測すると、1000個の浮動小数点数を足すのに
- 192cycle(zmmレジスタを使った主ループが8cycle×15.5回、マージしていく部分が36cycle、端数処理が4cycle×8回)
- 213命令(ループの外が27命令、ループの中が12命令×15.5回)
- 32分岐命令
が必要です。
これはperf
コマンドの結果と整合します。
理論限界と大きく差がついている理由は、Juliaの場合と全く同様です。
どうでもいいですが、配列の最後のほうがループ内定数とみなされてxmmレジスタに乗っていますね……。これは微妙にずるな気がしますが計算時間や命令数に変わりはないのでよしとしましょう。
C言語の強みは何といっても ハードウェアを作っている会社が専用コンパイラを提供してくれる ことにあります。
Intel社が作ったコンパイラ、icc
を使ってコンパイルしてみましょう。
コンパイルオプションは、-axCORE-AVX512 -qopt-zmm-usage=high
です。
結果は以下の通りです。
- 62.7G cycles
- 172G insns
- 15.6G branches
- 18秒くらい
理論限界の62.5Gサイクルにかなり近い速度が達成できていることがわかります*5。
最内ループを確認してみると、お互いに依存しない8個のvaddpd
命令になっていました。
さすがにこの辺の最適化はちゃんとしていますね。
まとめ