先週、浮動小数点数を足していくプログラムについて、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つ必要
端数処理が直列化されている
- 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回行われています。