先週、浮動小数点数を足していくプログラムについて、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();
}
使用したハードウェア
使用したソフトウェア
julia version 1.5.2（内部で使用しているLLVMのバージョンはLLVM version 9.0.1jl）clang version 10.0.0icc (ICC) 19.0.5.281 20190815perfgdb
理論限界
ベンチマークプログラムは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つ必要
- 端数処理が直列化されている
- zmmレジスタを使って並列化可能
C言語の結果（clangでコンパイルした場合）
コンパイルオプションは、-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言語の本気（iccでコンパイルした場合）
C言語の強みは何といっても ハードウェアを作っている会社が専用コンパイラを提供してくれる ことにあります。
Intel社が作ったコンパイラ、iccを使ってコンパイルしてみましょう。
コンパイルオプションは、-axCORE-AVX512 -qopt-zmm-usage=highです。
結果は以下の通りです。
- 62.7G cycles
- 172G insns
- 15.6G branches
- 18秒くらい
理論限界の62.5Gサイクルにかなり近い速度が達成できていることがわかります*5。
最内ループを確認してみると、お互いに依存しない8個のvaddpd命令になっていました。
さすがにこの辺の最適化はちゃんとしていますね。
まとめ
- Juliaの速度は
clangでコンパイルしたC言語と同程度 - 「単純和」プログラムをC言語で書いて
iccでコンパイルした場合、clangでコンパイルした場合やJuliaで書いた場合より3倍くらい速い - コンパイラが自動でやってくれる程度のSIMD化をする手間は、JuliaもC言語も大して変わらないと思う
*1:この手の言語間の速度比較では、不適切なソースコードが使われることにより特定の言語の性能が不当に低く出てしまうということがありがちですが、そういうことにならない有意義な比較ができることを感謝いたします。
*2:-O3や-Ofastだと完全にアンローリングされてしまって見づらかったので-O2にしました（実行にかかるサイクル数は同じでした）。
*3:-funsafe-math-optimizationsは-ffast-mathに含まれるので、-ffast-mathでもよいです。
*4:-funsafe-math-optimizationsは複数のオプションを有効にしますが、実際に役立っているのは-fassociative-math -fno-signed-zerosのようです。
*5:ループ交換が行われているので本来とは全く異なる順序で計算されていますが、加算は1000G回行われています。
