昨日ついったーで面白い実験をみかけたので、ちょっと考察してみました。
元記事の方は、JuliaとCで同じ処理を行うプログラム書いて性能を比較されてたんですが、 私としてはJuliaとの比較ではなくHWスペックに対してC言語の実験結果が異様に遅いことに興味をひかれました。
使われていたプログラムは、1000要素の倍精度の配列の値を10000000足し合わせるというもので、 Intel(R) Xeon(R) CPU E5-2603 v3 @ 1.60GHzで測定されたとのことです。
とりあえず、メモリネックなのか演算ネックなのか確認したいので、HWのB/Fを計算してみましょう。
こちらにあるようにE5-2603 V3はHaswell世代の6コアのCPUですが、今回は1コアしか使ってないようなのでHW的には1コア分の性能だけを見ます。
L1Dキャッシュは、32KBなので、1000要素の配列なら全データがオンキャッシュで動作します。初回アクセス時にキャッシュミスは発生しますが、 1千万回も繰り返しているので性能に与える影響はほぼ無視できるはずです。Haswel世代のL1Dキャッシュは、load時は64Bytes/cycleとのことなので、この値をメモリアクセス性能とします*1
演算性能は、1コアで256bit幅のFMA2基ですが今回は加算しか演算が無いのでただの加算器が2つと考えて、8Flop/cycleとします。
したがって、この場合のHWのB/F値は8となります。
一方、プログラム側のB/F値は8Byteロードに対して1演算なので8B/Fとなります。
ということは、このプログラムは1コア内の加算器を限界まで使えるので、理論上の性能の限界は2x4x1.6 = 12.8GFlops となります。 演算数から時間を逆算すると、1000x10000000/(12.8x10003) = 0.78 sec ですね。
こちらのツイートで示されている、Juliaでまっとうな書き方をした時の性能が0.89secということなので、もし同じようなCPUでの測定結果であれば 少なくともこのケースに限ってはJuliaの最適化は優秀と言えそうです。
#Julia言語
— 黒木玄 Gen Kuroki (@genkuroki) 2020年9月27日
Juliaでは単純な和には普通にsimdを使うので、JuliaとCの両方でsimdを使って比較した方が良かったと思います。
Juliaでsimdを使うと10倍以上速くなります。
Juliaのbuiltin sumも同じオーダーで速いです。
ソースコード↓https://t.co/tB5Ti3Hq4A pic.twitter.com/C31z3YPJu7
さて、C言語の方ですが、元記事ではgccで-O3でコンパイルして18.8secという結果だったそうです。
アセンブリも示されているのでまずそちらを見ると、加算命令としてaddsd命令が使われていることがわかります。 つまり、SIMD演算を使っておらず8Byteづつ加算している状態です。 また、twitterで元の記事を書かれた方から指摘をいただいて気付いたんですが、配列の先頭から順に加算するため、 せっかく2ポートで同時に加算命令が発行できるのにこの機能も使えません。 したがって、さきほど計算した理論性能と比較すると、SIMDが使えないことで1/4さらに命令発行が1命令/cycleとなることで1/2になり合わせて1/8の性能しか使えません。
そうすると実行時間で考えるとだいた6.4秒くらいですね。 実測値は、18.8秒なのでさらに3倍近く遅いのですが、おそらくこれは、ADDSD命令のlatencyとThroughputの差と思われます。
これまで1cycle毎に命令が行なわれる*2前提で考えてきたので、1演算あたり1cycleとしていましたが、先行する加算命令の結果を書かれたレジスタから値を読んで後続の加算命令が行なわれているので、ここはlatencyである3cycleが1演算あたり必要な状態と考えるべきです。*3 したがって、さらに3倍遅くて19秒くらい(なんかちょっと多いな・・・)という計算で辻褄が合いそうです。
最適化
さて、ここまで、性能低下の要因を色々とみてきましたが原因は全部同じで、「1要素づつしか加算していない」ということにつきます。*4 じゃ、どうすりゃ良いのかというと
for(i =0; i<1000; i++){
s += a[i];
}
というような最内ループをストリップマイニングして
for(i=0; i<1000; i+=8){
s0 += a[i];
s1 += a[i+1];
s2 += a[i+2];
s3 += a[i+3];
s4 += a[i+4];
s5 += a[i+5];
s6 += a[i+6];
s7 += a[i+7];
}
s=s0+s1+s2+s3+s4+s5+s6+s7
という形に書き換えるようにコンパイラに指示すれば良いわけです。
gccの場合はその名もずばりの-floop-strip-mineというのがあるので、これと(SIMD化のために)-march=haswellを指定してコンパイルすればたぶん最適化してくれるでしょう。
実際に(全体をコンパイルすると読むのが面倒なので)次のような2重ループの関数を作って、-Sでアセンブリを吐いてみます。
#include <stddef.h>
double sum(double* a){
double s=0;
for(size_t j=0; j<10000000;j++){
for(size_t i =0; i<1000; i++){
s += a[i];
}
}
return s;
}
-O3だけの時
> docker run -it --rm -w /usr/src -v ${PWD}:/usr/src gcc:9.1.0 gcc -S -O3 test.c -o test_O3.s
> cat test_O3.s
.file "test.c"
.text
.p2align 4
.globl sum
.type sum, @function
sum:
.LFB0:
.cfi_startproc
movq %rdi, %rsi
movl $10000000, %ecx
pxor %xmm0, %xmm0
leaq 8000(%rdi), %rdx
.L2:
movq %rsi, %rax
.p2align 4,,10
.p2align 3
.L3:
movsd (%rax), %xmm1
addq $16, %rax
addsd %xmm0, %xmm1
movsd -8(%rax), %xmm0
addsd %xmm1, %xmm0
cmpq %rax, %rdx
jne .L3
subq $1, %rcx
jne .L2
ret
.cfi_endproc
.LFE0:
.size sum, .-sum
.ident "GCC: (GNU) 9.1.0"
.section .note.GNU-stack,"",@progbits
-march=haswellと-floop-strip-mineをつけると
> docker run -it --rm -w /usr/src -v ${PWD}:/usr/src gcc:9.1.0 gcc -S -O3 -march=haswell -floop-strip-mine test.c -o test_O3_haswell_stripmine.s
> cat test_O3_haswell_stripmine.s
.file "test.c"
.text
.section .text.unlikely,"ax",@progbits
.LCOLDB1:
.text
.LHOTB1:
.p2align 4
.globl sum
.type sum, @function
sum:
.LFB0:
.cfi_startproc
movq %rdi, %rax
vxorpd %xmm1, %xmm1, %xmm1
leaq 8000(%rdi), %rsi
xorl %ecx, %ecx
.L3:
movq %rax, %rdx
.L2:
vmovupd (%rdx), %xmm2
vmovupd (%rdx), %ymm3
addq $32, %rdx
vaddsd %xmm1, %xmm2, %xmm1
vunpckhpd %xmm2, %xmm2, %xmm2
vextractf128 $0x1, %ymm3, %xmm0
vaddsd %xmm2, %xmm1, %xmm1
vaddsd %xmm0, %xmm1, %xmm1
vunpckhpd %xmm0, %xmm0, %xmm0
vaddsd %xmm0, %xmm1, %xmm1
cmpq %rdx, %rsi
jne .L2
incl %ecx
cmpl $10000000, %ecx
je .L6
jmp .L3
.cfi_endproc
.section .text.unlikely
.cfi_startproc
.type sum.cold, @function
sum.cold:
.LFSB0:
.L6:
vmovapd %xmm1, %xmm0
vzeroupper
ret
.cfi_endproc
.LFE0:
.text
.size sum, .-sum
.section .text.unlikely
.size sum.cold, .-sum.cold
.LCOLDE1:
.text
.LHOTE1:
.ident "GCC: (GNU) 9.1.0"
.section .note.GNU-stack,"",@progbits
(・3・)あるぇ～？
3オペランド方式のvバージョンになっていますが、あいかわらずaddsdが使われているようです。
ちなみに、icc 19でデフォルトオプションだとこんな感じになります。
# mark_description "-std=c11 -S";
.file "test.c"
.text
..TXTST0:
.L_2__routine_start_sum_0:
# -- Begin sum
.text
# mark_begin;
.align 16,0x90
.globl sum
# --- sum(double *)
sum:
# parameter 1: %rdi
..B1.1: # Preds ..B1.0
# Execution count [0.00e+00]
.cfi_startproc
..___tag_value_sum.1:
..L2:
#2.22
xorl %edx, %edx #5.5
pxor %xmm0, %xmm0 #3.11
pxor %xmm1, %xmm1 #3.11
xorl %eax, %eax #4.3
# LOE rdx rbx rbp rdi r12 r13 r14 r15 eax xmm0 xmm1
..B1.2: # Preds ..B1.4 ..B1.1
# Execution count [6.25e+05]
movaps %xmm1, %xmm10 #3.11
movl %eax, %ecx #4.3
movaps %xmm10, %xmm9 #3.11
movaps %xmm9, %xmm8 #3.11
movaps %xmm8, %xmm7 #3.11
movaps %xmm7, %xmm6 #3.11
movaps %xmm6, %xmm5 #3.11
movsd (%rdi,%rdx,8), %xmm2 #6.13
movaps %xmm5, %xmm4 #3.11
unpcklpd %xmm2, %xmm2 #6.13
movaps %xmm4, %xmm3 #3.11
# LOE rdx rbx rbp rdi r12 r13 r14 r15 eax ecx xmm0 xmm1 xmm2 xmm3 xmm4 xmm5 xmm6 xmm7 xmm8 xmm9 xmm10
..B1.3: # Preds ..B1.3 ..B1.2
# Execution count [6.25e+08]
addl $16, %ecx #4.3
addpd %xmm2, %xmm10 #6.8
addpd %xmm2, %xmm9 #6.8
addpd %xmm2, %xmm8 #6.8
addpd %xmm2, %xmm7 #6.8
addpd %xmm2, %xmm6 #6.8
addpd %xmm2, %xmm5 #6.8
addpd %xmm2, %xmm4 #6.8
addpd %xmm2, %xmm3 #6.8
cmpl $10000000, %ecx #4.3
jb ..B1.3 # Prob 99% #4.3
# LOE rdx rbx rbp rdi r12 r13 r14 r15 eax ecx xmm0 xmm1 xmm2 xmm3 xmm4 xmm5 xmm6 xmm7 xmm8 xmm9 xmm10
..B1.4: # Preds ..B1.3
# Execution count [6.25e+05]
addpd %xmm9, %xmm10 #3.11
addpd %xmm7, %xmm8 #3.11
addpd %xmm5, %xmm6 #3.11
addpd %xmm3, %xmm4 #3.11
addpd %xmm8, %xmm10 #3.11
addpd %xmm4, %xmm6 #3.11
addpd %xmm6, %xmm10 #3.11
movaps %xmm10, %xmm2 #3.11
incq %rdx #5.5
unpckhpd %xmm10, %xmm2 #3.11
addsd %xmm2, %xmm10 #3.11
addsd %xmm10, %xmm0 #3.11
cmpq $1000, %rdx #5.5
jb ..B1.2 # Prob 99% #5.5
# LOE rdx rbx rbp rdi r12 r13 r14 r15 eax xmm0 xmm1
..B1.5: # Preds ..B1.4
# Execution count [7.63e-02]
ret #9.10
.align 16,0x90
# LOE
.cfi_endproc
# mark_end;
.type sum,@function
.size sum,.-sum
..LNsum.0:
.data
# -- End sum
.data
.section .note.GNU-stack, ""
# End
なんかめっちゃアンロールしてますね。 -qopt-report=5で出力される最適化レポートを見ると
ループの開始 test.c(5,5)
リマーク #25444: ループの入れ子の交換: ( 1 2 ) --> ( 2 1 )
リマーク #15542: ループ はベクトル化されませんでした: 内部ループがすでにベクトル化されています。 [ test.c(5,5) ]
ループの開始 test.c(4,3)
リマーク #15305: ベクトル化のサポート: ベクトル長 2
リマーク #15399: ベクトル化のサポート: アンロールファクターが 8 に設定されます。
リマーク #15309: ベクトル化のサポート: 正規化されたベクトル化のオーバーヘッド 0.525
リマーク #15301: 変換されたループ がベクトル化されました。
リマーク #15475: --- ベクトルのコストサマリー開始 ---
リマーク #15476: スカラーのコスト: 6
リマーク #15477: ベクトルのコスト: 2.500
リマーク #15478: スピードアップの期待値: 2.390
リマーク #15488: --- ベクトルのコストサマリー終了 ---
リマーク #25015: ループの最大トリップカウントの予測=625000
ループの終了
ループの終了
あ、入れ替えられてる・・・そらそうかorz
