その2で宣言した通り、今回の対象はループだ。その2に続きが発生したので番号的には飛んでしまっているが、続きという点は変わらない。 しかしループとなると絶対にSIMDとループアンロールが絡んで面倒なことになる、と思っているので少し腰が重いが、最適化の勉強というならむしろこれこそやらねばならない。
今回とりあえずこんなコードを書いた。
#include <stdint.h> int32_t accum(const int32_t* begin, const int32_t* end) { int32_t sum = 0; for(; begin != end; ++begin) sum += *begin; return sum; }
begin
からend
が到達可能ならそこまでの範囲を総和する。到達不能なら知らない。
最適化なし
では最適化なしでできたバイナリを逆アセンブルしよう。
$ gcc -std=c99 -c -O0 accumulate.c && objdump -d accumulate.o > accumulate.asm
すると、
accumulate.o: file format elf64-x86-64 Disassembly of section .text: 0000000000000000 <accum>: 0: push %rbp ; いつもの 1: mov %rsp,%rbp ; いつもの 4: mov %rdi,-0x18(%rbp); 第一引数をスタックへ 8: mov %rsi,-0x20(%rbp); 第二引数をスタックへ c: movl $0x0,-0x4(%rbp) ; sumを即値0で初期化 13: jmp 23 <accum+0x23> ; 23へ無条件ジャンプ 15: mov -0x18(%rbp),%rax; beginの持つアドレス値をraxへロード 19: mov (%rax),%eax ; raxの指すメモリの内容をeaxへロード 1b: add %eax,-0x4(%rbp) ; eaxをsumへ足す 1e: addq $0x4,-0x18(%rbp); beginに即値4(sizeof(int32_t))を足す 23: mov -0x18(%rbp),%rax; 現在のポインタbeginをraxレジスタへロード 27: cmp -0x20(%rbp),%rax; endとbeginを比較。同一ならZF==1 2b: jne 15 <accum+0x15> ; ZF==0なら15まで戻る 2d: mov -0x4(%rbp),%eax ; sumを返却値レジスタへ 30: pop %rbp ; いつもの 31: retq
既にかなり短い。そして、ここまでやった内容で十分わかる。で、ここまでの経験からこれを最適化するとスタックメモリを使わないように頑張ると予想され、以下のようなコードが想定される(動作確認とかしてないのでミスや誤解があるかも)。
accumulate.o: file format elf64-x86-64 Disassembly of section .text: 0000000000000000 <accum>: movl $0x0,%eax jmp judge inside_loop: add (%rdi),%eax addq $0x4,%rdi judge: cmp %rsi,%rdi jne inside_loop <accum+0x15> retq
だが、現実は甘くはなかった。
-O3
では最適化しよう。
$ gcc -std=c99 -c -O3 accumulate.c && objdump -d accumulate.o > accumulate_O3.asm
結果、目を疑うようなコードが出てきた。
accumulate.o: file format elf64-x86-64 Disassembly of section .text: 0000000000000000 <accum>: 0: cmp %rsi,%rdi 3: je 128 <accum+0x128> 9: lea 0x4(%rdi),%rcx d: mov %rsi,%r9 10: mov %rdi,%rdx 13: and $0xf,%edx 16: sub %rcx,%r9 19: shr $0x2,%rdx 1d: shr $0x2,%r9 21: neg %rdx 24: lea 0x1(%r9),%r8 28: and $0x3,%edx 2b: cmp %r8,%rdx 2e: cmova %r8,%rdx 32: cmp $0x8,%r8 36: ja 138 <accum+0x138> 3c: mov %r8,%rdx 3f: cmp $0x1,%rdx 43: mov (%rdi),%eax 45: je 9c <accum+0x9c> 47: add 0x4(%rdi),%eax 4a: cmp $0x2,%rdx 4e: lea 0x8(%rdi),%rcx 52: je 9c <accum+0x9c> 54: add 0x8(%rdi),%eax 57: cmp $0x3,%rdx 5b: lea 0xc(%rdi),%rcx 5f: je 9c <accum+0x9c> 61: add 0xc(%rdi),%eax 64: cmp $0x4,%rdx 68: lea 0x10(%rdi),%rcx 6c: je 9c <accum+0x9c> 6e: add 0x10(%rdi),%eax 71: cmp $0x5,%rdx 75: lea 0x14(%rdi),%rcx 79: je 9c <accum+0x9c> 7b: add 0x14(%rdi),%eax 7e: cmp $0x6,%rdx 82: lea 0x18(%rdi),%rcx 86: je 9c <accum+0x9c> 88: add 0x18(%rdi),%eax 8b: cmp $0x8,%rdx 8f: lea 0x1c(%rdi),%rcx 93: jne 9c <accum+0x9c> 95: add 0x1c(%rdi),%eax 98: lea 0x20(%rdi),%rcx 9c: cmp %rdx,%r8 9f: je 130 <accum+0x130> a5: sub %rdx,%r8 a8: sub %rdx,%r9 ab: lea -0x4(%r8),%r10 af: shr $0x2,%r10 b3: add $0x1,%r10 b7: cmp $0x2,%r9 bb: lea 0x0(,%r10,4),%r11 c2: c3: jbe 109 <accum+0x109> c5: pxor %xmm0,%xmm0 c9: lea (%rdi,%rdx,4),%rdi cd: xor %edx,%edx cf: add $0x1,%rdx d3: paddd (%rdi),%xmm0 d7: add $0x10,%rdi db: cmp %rdx,%r10 de: ja cf <accum+0xcf> e0: movdqa %xmm0,%xmm1 e4: lea (%rcx,%r11,4),%rcx e8: psrldq $0x8,%xmm1 ed: paddd %xmm1,%xmm0 f1: movdqa %xmm0,%xmm1 f5: psrldq $0x4,%xmm1 fa: paddd %xmm1,%xmm0 fe: movd %xmm0,%edx 102: add %edx,%eax 104: cmp %r11,%r8 107: je 12a <accum+0x12a> 109: lea 0x4(%rcx),%rdx 10d: add (%rcx),%eax 10f: cmp %rdx,%rsi 112: je 12a <accum+0x12a> 114: lea 0x8(%rcx),%rdx 118: add 0x4(%rcx),%eax 11b: cmp %rdx,%rsi 11e: je 150 <accum+0x150> 120: add 0x8(%rcx),%eax 123: retq 124: nopl 0x0(%rax) 128: xor %eax,%eax 12a: repz retq 12c: nopl 0x0(%rax) 130: repz retq 132: nopw 0x0(%rax,%rax,1) 138: test %rdx,%rdx 13b: jne 3f <accum+0x3f> 141: mov %rdi,%rcx 144: xor %eax,%eax 146: jmpq a5 <accum+0xa5> 14b: nopl 0x0(%rax,%rax,1) 150: repz retq
な ん だ こ れ は
基本的に命令の数が少ないほど速度は向上すると思っていたので、これは以外だ。確かにループを一回一回回るよりも効率のよい(と思われる)方法はたくさんある。冒頭でも言及したループアンローリングやSIMD命令などだ。だがここまでくると本当にこの方が速いのか少し疑わしくなる。
ついでに、くやしいので私が予想したコードを出しそうなオプションを使ってみる。-Os
だ。これはバイナリのサイズを小さくするように働く。
$ gcc -std=c99 -c -Os accumulate.c && objdump -d accumulate.o > accumulate_Os.asm
accumulate.o: file format elf64-x86-64 Disassembly of section .text: 0000000000000000 <accum>: 0: xor %eax,%eax 2: cmp %rsi,%rdi 5: je f <accum+0xf> 7: add (%rdi),%eax 9: add $0x4,%rdi d: jmp 2 <accum+0x2> f: retq
私の書いたコードと異なる点がいくつかある。前回も見たが、ゼロクリアは即値の代入でなくxor
になっている。初っ端からジャンプするのではなく、先にcmp
とje
を発行している。私が書いたコードと比べて、渡された領域のサイズが0だった場合にジャンプが一回少なくなる。だが基本的には同じ発想だ。
また、O2
も試した。
accumulate.o: file format elf64-x86-64 Disassembly of section .text: 0000000000000000 <accum>: 0: xor %eax,%eax 2: cmp %rsi,%rdi 5: je 1d <accum+0x1d> 7: nopw 0x0(%rax,%rax,1) e: 10: add (%rdi),%eax 12: add $0x4,%rdi 16: cmp %rdi,%rsi 19: jne 10 <accum+0x10> 1b: repz retq 1d: repz retq
こちらも同様だ。最低限の命令数で素直にループを回っている。
マイクロベンチマーク
O3
のあのとても長いバイナリで本当に高速になるとは俄には信じられないのでちょっとしたマイクロベンチマークをする。
#include "accumulate.h" #include <stdio.h> #include <stddef.h> int main() { const size_t N = 100000; int32_t a[N]; for(int32_t i=0; i<N; ++i) a[i] = i; for(size_t i = 0; i < N; ++i) { int32_t b = accum(a, a+N); } return 0; }
前回見たとおりLTOなしで別々にコンパイルした後リンクすればインライン化の恐れはない(と思いたかったが、一応確認した)。
ここで、-march=native
も戯れにコンパイルし、逆アセンブルした。余裕があれば後で持ち出すが、SIMDまわりの違いだった。結果は、
$ time ./opts # -Osでコンパイルした 3.23user 0.00system 0:03.23elapsed 99%CPU $ time ./opt3 # -O3でコンパイルした 1.66user 0.00system 0:01.66elapsed 99%CPU $ time ./optnative # -O3 -march=nativeでコンパイルした 1.18user 0.00system 0:01.18elapsed 99%CPU
となった。実際ハヤイ。ちなみに3回くらい試したが、ほぼこの関係は変わらなかった。
解読・前半戦
こうなるとコンパイラの出した答えに納得するしかないので最適化済みのアセンブリに戻ろう。とりあえずもう一度ざっと見て、何が起きているか把握したい。そして、できれば分割して解読したい。長くなると自分の考えを挟みづらいからだ。 しかし長いので、予想を立てる。多分、ここまで何度か言及しているループアンロールとSIMDを使いこなしているに違いない。というかそれ以外では並列化しかこの手のループが早くなる方法を私は知らない。そして多分並列化はされていない。
この予想を念頭に置いてもう一度前からざっと見ていこう。
まず、最初にポインタに対する演算をいくつかしている。その後、似たような命令のパターンが7回ほど続き、それに続く少し複雑そうな命令の中でpaddd
が使われている。確定だ。SIMD命令が使われている。
後ろでSIMD命令を使うことを考え、先に何回SIMD演算をすることになるかの計算と、あまりの分の処理をしておいているのだろう。
大まかな見当がついたので、分割して順に見ていこう。最初の部分は最後尾へのジャンプをいくつか含むので最後尾も一緒に見る。
0000000000000000 <accum>: 0: cmp %rsi,%rdi ; 引数同士を比較 3: je 128 <accum+0x128> ; 同じなら範囲が空なので0を返す ; ここ以降では要素数が1以上と確定する 9: lea 0x4(%rdi),%rcx ; rdi+4をrcxへ代入 d: mov %rsi,%r9 ; endをr9レジスタへコピー 10: mov %rdi,%rdx ; beginをrdxレジスタへコピー 13: and $0xf,%edx ; edxの下4bitだけを残す 16: sub %rcx,%r9 ; r9(end)からrcx(begin+1)を引く 19: shr $0x2,%rdx ; rdxを2回右シフト 1d: shr $0x2,%r9 ; r9を2回右シフト 21: neg %rdx ; rdxの符号を反転 24: lea 0x1(%r9),%r8 ; r9+0x1をr8へ 28: and $0x3,%edx ; edxの下位2bitだけを残す 2b: cmp %r8,%rdx ; r8とrdxを比較 2e: cmova %r8,%rdx ; if(r8 < rdx) rdx = r8 32: cmp $0x8,%r8 ; 8とr8を比較 36: ja 138 <accum+0x138> ; if(8 < r8) goto 138 ;... 128: xor %eax,%eax ; eaxをゼロクリア 12a: repz retq ; return 12c: nopl 0x0(%rax) ; 130: repz retq ; return 132: nopw 0x0(%rax,%rax,1) ; 138: test %rdx,%rdx ; rdxが0ならZF=1 13b: jne 3f <accum+0x3f> ; ZF == 0なら3fへジャンプ 141: mov %rdi,%rcx ; rdiをrcxへ代入 144: xor %eax,%eax ; eaxをゼロクリア 146: jmpq a5 <accum+0xa5> ; a5へジャンプ 14b: nopl 0x0(%rax,%rax,1) ; 150: repz retq
最初の2行はわかりやすい。範囲が空だった場合に即リターンするためのものだ。そして、続く部分では要素が少なくとも1つ以上あるという前提で話を進めている。
先頭では2つのことが同時に行われているようで、意識があっちこっちに振り回されて読みにくい。とりあえず片方を無視してわかりやすい方だけ追っていこう。そちら側だけ抽出するとこうなる。
9: lea 0x4(%rdi),%rcx ; begin+1をrcxへ代入 d: mov %rsi,%r9 ; endをr9レジスタへコピー 16: sub %rcx,%r9 ; r9(end)からrcx(begin+1)を引く 1d: shr $0x2,%r9 ; r9((要素数-1)*4)を2回右シフト 24: lea 0x1(%r9),%r8 ; r9+0x1をr8へ
まず9:
でbegin+1
を用意している(ここでsizeof(int32_t) == 4
なので、0x4(%rdi)
はbegin+1
を意味していることになる)。
続いて、d:
でend
をr9
へコピーし、16:
でそこからbegin+1
を引いている。これで、(要素数-1)*4
の値が得られた(この最後の4倍はsizeof(int32_t)
だ)。
これを1d:
で2回右シフトすることで、要素数-1を得ている。右シフト2回は4での除算になるからだ。そして24:
で要素数そのものがr8
へ代入されている。
得られた値がわかりやすいのでこちらが何をしているかはだいたいわかったと思っていいだろう。反対側は何をしているんだろう。
10: mov %rdi,%rdx ; beginをrdxレジスタへコピー 13: and $0xf,%edx ; edxの下4bitだけを残す 19: shr $0x2,%rdx ; rdxを2回右シフト 21: neg %rdx ; rdxの符号を反転 28: and $0x3,%edx ; edxの下位2bitだけを残す
じっくり何が起きるか考えてみないことには何が得られるかすらわからない。使われ方を追いかけてみたが、色々なところにジャンプするばかりでそもそもこの値が何かわからないと理解の足しにはならなかった。
そもそも、これはbegin
のアドレス値に対してなにかしらを行っているので、要素の数とかそういうことはわからないはずだ(endの値を使っていないので)。だとするとなんだろう。
アドレスそのものの値に関係し、ループするときに、SIMDの前にしておく必要があり、その値によっては先に要素を処理しておいたりしなければならない……。ここで思いついた。メモリアライメントだ。
アライメントだとすれば、下位4bitだけ残した後に2回右シフトしてその後下位2bitだけ残した理由がわかる(符号反転は直感ではわからないが)。16バイト境界にアラインしているか確認しているのだ。
本当にそうだろうか。それは、今後のメモリの使われ方を確認すればわかる。SIMD命令の初っ端を確認すればいい。メモリからロードしているなら、整列していることが要求される。
d3: paddd (%rdi),%xmm0
rdi
のアドレスが直接使われている。paddd
命令は32bit整数4つを一度に加算する命令だ。ソースオペランドがメモリアドレスなら、それは16バイト境界に整列していなければならない。
ここでは要素数の余りの他にも、メモリアライメントの余りも考慮したバイナリが生成されているのだ。後でのSIMD命令を可能な限り効率よくするために!
それによく見てみれば、SIMD演算が終了した後にひっそりと、通常のadd
命令が含まれている。これは、アライメントされた領域を脱した後に要素が残っているならそれを足す、という操作をしているのだ。
109: lea 0x4(%rcx),%rdx ; rcxから32bit進んだアドレスをrdxへ 10d: add (%rcx),%eax ; 返却値へ余った要素を足す 10f: cmp %rdx,%rsi ; 先のrdxとendを比較 112: je 12a <accum+0x12a> ; 一致していればreturn 114: lea 0x8(%rcx),%rdx ; rcxから2つ分進んだアドレスをrdxへ 118: add 0x4(%rcx),%eax ; 返却値へ余った要素を足す 11b: cmp %rdx,%rsi ; rdxとendと比較 11e: je 150 <accum+0x150> ; 一致していればreturn 120: add 0x8(%rcx),%eax ; 最大でも3つしか余らないので 123: retq ; 足してreturn
ちなみに、中盤の同じ命令のパターンが続く部分も、基本的にはこれと同じことをしている。16バイト境界に整列したとみなせる部分まで、手動で足しているのだ。
しかし、それなら最大でも3つのはずで(int32_t
は32bit、つまり4バイトなので、16バイトでアラインするまで4つは超えない)、7つある理由がわからない……と思っていたが、そんなことはない。
SIMDで4つずつ足すとしよう。手元に4つの要素がある。これをSIMDにロードして……それでどうする? 配列の要素は既に全てSIMDレジスタにロードされた。これ以上することはない。
あとは総和を取るだけだが、これは素でやってもそう変わらない。なので、SIMDにして高速化を期待するためには、8個以上の要素数が必要なのだ。だから最初に7個分の処理があるのだろう。
解読・後半戦
最後に、SIMDでループしている部分、先のアセンブリコードの直前の部分を読もう。
c5: pxor %xmm0,%xmm0 ; xmm0をゼロクリア c9: lea (%rdi,%rdx,4),%rdi ; rdi = rdi + rdx * 4 cd: xor %edx,%edx ; edxをゼロクリア cf: add $0x1,%rdx ; rdx += 1(カウンタインクリメント) d3: paddd (%rdi),%xmm0 ; xmm0へrdiから128bit分足す d7: add $0x10,%rdi ; rdiを16byte前進 db: cmp %rdx,%r10 ; rdxをr10(多分4で割った余り)と比較 de: ja cf <accum+0xcf> ; 終了していなければcfへ e0: movdqa %xmm0,%xmm1 ; xmm0をxmm1へコピー e4: lea (%rcx,%r11,4),%rcx ; rcx += r11*4 e8: psrldq $0x8,%xmm1 ; xmm1を8byteシフト ed: paddd %xmm1,%xmm0 ; xmm1をxmm0に足す f1: movdqa %xmm0,%xmm1 ; xmm0をxmm1へコピー f5: psrldq $0x4,%xmm1 ; xmm1を4byteシフト fa: paddd %xmm1,%xmm0 ; xmm1をxmm0に足す fe: movd %xmm0,%edx ; xmm0の下位ビットをedxへ 102: add %edx,%eax ; edxをeaxへ足す 104: cmp %r11,%r8 ; r11とr8を比較 107: je 12a <accum+0x12a> ; 一致すればリターン
ここでは、メモリアラインが済んでいると考えてSIMDループを回した後、総和を取って返却値へ足している。総和はシフトしつつ足すという方法によって行われている。 8byteシフトすると、上2つのint32_tと下2つのint32_tがそれぞれ足されることになり、その後その値をまたシフトして足せば総和が求まる。 これを、余った部分を既に足していた返却値レジスタへ足し、この領域を抜けている。この後に続くのは最後に余った部分の総和だ。
bit単位で完全に理解しながら読んできたわけではないが、コンパイラが異様に賢いことはよくわかった。一応、-march=native
や手動ループアンロールのコードも用意していたのだが、結構長くなったし、それ以上に私が疲れた。
今回はここまでで、今手持ちのストックは日を改めて解読していくことにする。
では。
(追記)
しかし、このコードだと要素数が非常に小さかったときはかえって遅くなるのでは、と考えた私は、7要素(SIMDにできない限界だ)の配列を作り、-O0
、-O3
、-O3 -march=native
で比較してみた。すると、
$ gcc -std=c99 -O0 main.c accumulate.c $ time ./a.out 0.22user 0.00system 0:00.22elapsed 98%CPU $ gcc -std=c99 -O2 main.c accumulate.c $ time ./a.out 0.08user 0.00system 0:00.08elapsed 98%CPU $ gcc -std=c99 -O3 main.c accumulate.c $ time ./a.out 0.10user 0.00system 0:00.10elapsed 99%CPU $ gcc -std=c99 -O3 -march=native main.c accumulate.c $ time ./a.out 0.07user 0.00system 0:00.07elapsed 98%CPU
妥当な結果が出た。O0
だとスタック読み書きの効果が出て遅いが、-O2
だとその分がなくなるので、SIMDの恩恵が受けられない場合はO3
やmarch=native
とも十分張り合える速度になる。
もう少しよく見てみたいと思ったので、ループ回数を10倍にしたところ、
$ gcc -std=c99 -O2 main.c accumulate.c $ time ./a.out 0.46user 0.00system 0:00.47elapsed 99%CPU 0.45user 0.00system 0:00.45elapsed 99%CPU 0.46user 0.00system 0:00.46elapsed 100%CPU $ gcc -std=c99 -O3 main.c accumulate.c $ time ./a.out 0.55user 0.00system 0:00.55elapsed 99%CPU 0.57user 0.00system 0:00.57elapsed 100%CPU 0.57user 0.00system 0:00.57elapsed 99%CPU $ gcc -std=c99 -O3 -march=native main.c accumulate.c $ time ./a.out 0.50user 0.00system 0:00.51elapsed 99%CPU 0.51user 0.00system 0:00.51elapsed 99%CPU 0.48user 0.00system 0:00.49elapsed 99%CPU
となった。若干ではあるが(そして3回ずつしか実行していないので統計的にはあやしいが)、O2
が一番早くなっている。SIMDが使えない条件だと、メモリ配置や要素数の計算はオーバーヘッドでしかない。