SSEを使って8flops/clockを実現する
この記事の記述には誤りがあることが分かっています．代わりにこちらの記事をご覧ください．
SSE(Streaming SIMD Extensions)では，1clockに4つの単精度浮動小数点演算が
可能で，SSE2を使うと1clockに2つの倍精度浮動小数点演算も可能です．
しかし，一部のプロセッサでは，乗算と加算を同時に行なうことができ，
これにより 8flops (floating-point operations) / clock の単精度浮動小数点演算が実現できます． *1
使用したCPU
この記事の内容は，
で調査・実験しています．
それ以外の環境では，記述が当てはまらないことも考えられます．
もしお気づきの点がございましたらコメント頂けたら嬉しいです．
乗算と加算を並列に行なうコツ
SSEでは，128-bitのxmmレジスタを使って計算を進めます．
8flops/clockを実現するためにはこのレジスタの使い方が非常に重要です．
例えば，次のように各命令に依存関係が発生しないようにレジスタを使い回しても
4flops/clockがせいぜいです．
(この記事でのアセンブリの表記は，GCCのインラインアセンブラと同じです．
"Op Reg1, Reg2" だと，結果は Reg2 に格納されると考えてください．)
mulps %xmm0, %xmm1 addps %xmm2, %xmm3 mulps %xmm4, %xmm5 addps %xmm6, %xmm7 mulps %xmm8, %xmm9 addps %xmm10, %xmm11 mulps %xmm12, %xmm13 addps %xmm14, %xmm15
速度を出すためには，「mulpsの結果を直後のaddpsに渡す(或いはその逆)」ということが
重要なようです．*2
直感的には，
1. mulpsの結果が出る (1clock目)
2. mulpsの結果を結果をフォワーディングして，addpsの入力とし，addpsを行なう (2clock目)
3. mulpsを行なう (2clock目)
といった形でフォワーディングが行われているのだと推察できます．
具体的にはこういうコードを書けば良いです．
mulps %xmm0, %xmm1 addps %xmm1, %xmm15 mulps %xmm0, %xmm2 addps %xmm2, %xmm14 mulps %xmm0, %xmm3 addps %xmm3, %xmm13 mulps %xmm0, %xmm4 addps %xmm4, %xmm12
8flops/clock出るコード
(この記事の終わりに，コンパイルしてすぐに実行できるフルのコードも掲載します)
実際に，フォワーディングを使って 8flops/clock を実現したコードのメイン部分を記載します．
ざっくり読む上での注意点も述べておきます．
- print_GFLOPS, print_throughput は計測結果を出力する関数
- rdtsc() は，CPUのRDTSC命令を用いて，(ほぼ)CPUクロックベースで時間を計測するための関数
- mulps SRC, DEST で，ソースにメモリ上のデータを使う場合，そのデータは16バイト境界にアラインされている必要がある
コード
void sse_mulps_addps_forwarding() { const int flops_per_instruction = 4; const int instructions_per_loop = 8; const double flops = flops_per_instruction * instructions_per_loop * LOOP; uint64_t clk0, clk1, cycles; int i; float __attribute__ ((aligned(16))) a[4] = {0.0, 0.0, 0.0, 0.0}; printf("-- sse_mulps_addps_forwarding --\n"); zero_all_xmm(); clk0 = rdtsc(); for (i = 0; i < LOOP; ++i) { asm volatile ( "mulps %0, %%xmm0\n\t" "addps %%xmm0, %%xmm4\n\t" "mulps %0, %%xmm1\n\t" "addps %%xmm1, %%xmm5\n\t" "mulps %0, %%xmm2\n\t" "addps %%xmm2, %%xmm6\n\t" "mulps %0, %%xmm3\n\t" "addps %%xmm3, %%xmm7\n\t" : :"m"(a[0]) ); } clk1 = rdtsc(); cycles = clk1 - clk0; print_GFLOPS(flops, cycles); print_throughput(instructions_per_loop * LOOP, cycles); }
このコードで演算器をフルに使った 8flops/clock が出るということは，(aの乗った)L1キャッシュからデータをロードするコストは，xmmレジスタからデータをロードするコストに並んでいると考えられますね．
4flops/clockがせいぜいなコード
mulpsとaddpsの依存関係を排除したコードを記載します．
直感的にはこれが速そうですが，4flops/clock程度しか出ないことを確認しました．
void sse_mulps_addps_no_dependency() { const int flops_per_instruction = 4; const int instructions_per_loop = 8; const double flops = flops_per_instruction * instructions_per_loop * LOOP; uint64_t clk0, clk1, cycles; int i; printf("-- sse_mulps_addps_no_dependency --\n"); zero_all_xmm(); clk0 = rdtsc(); for (i = 0; i < LOOP; ++i) { asm volatile ("mulps %xmm0, %xmm1\n\t" "addps %xmm2, %xmm3\n\t" "mulps %xmm4, %xmm5\n\t" "addps %xmm6, %xmm7\n\t" "mulps %xmm0, %xmm1\n\t" "addps %xmm2, %xmm3\n\t" "mulps %xmm4, %xmm5\n\t" "addps %xmm6, %xmm7\n\t" ); } clk1 = rdtsc(); cycles = clk1 - clk0; print_GFLOPS(flops, cycles); print_throughput(instructions_per_loop * LOOP, cycles); }
すぐに実験できるコードと実験方法
以下のコードをコピペし， "sse_8flops_per_clock.c" など適当なファイルを作ってください．
その際，"#define GHz 2.00" の行を，お使いのCPUの動作周波数で書き換える必要があります．
コンパイル方法は，
$ gcc -O3 -funroll-loops sse_8flops_per_clock.c -o sse_8flops_per_clock
辺りだと性能が出るはずです．
特にループ・アンローリングのための "-funroll-loops" は性能に直結します．
実行すると， sse_mulps_addps_forwarding() では 8flops/clock 程度， sse_mulps_addps_no_dependency() では 4flops/clock 程度出ていることが確認できるかと思います．
Intel Xeon CPU E7540 @ 2.00GHz (後述のDVFS,TurboBoost共に無効) での実験結果
$ ./sse_8flops_per_clock -- sse_mulps_addps_forwarding -- GFLOPS @ 2.00GHz: 7.986 [flops/clock] = 15.972 [GFLOPS] (67108864 flops in 8403324 clock = 0.004202 sec) Throughput: 1.996 [instructions/clock] (16777216 instrucions in 8403324 clock) -- sse_mulps_addps_no_dependency -- GFLOPS @ 2.00GHz: 3.993 [flops/clock] = 7.985 [GFLOPS] (67108864 flops in 16808016 clock = 0.008404 sec) Throughput: 0.998 [instructions/clock] (16777216 instrucions in 16808016 clock)
Intel Core2 Duo CPU P9400 @ 2.40GHz (後述のDVFS無効,TurboBoost有効) での実験結果
(TurboBoostのせいで見かけ上の効率が理論値よりも大きくなっています)
$ ./sse_8flops_per_clock -- sse_mulps_addps_no_dependency -- GFLOPS @ 2.40GHz: 4.127 [flops/clock] = 9.905 [GFLOPS] (67108864 flops in 16261029 clock = 0.006775 sec) Throughput: 1.032 [instructions/clock] (16777216 instrucions in 16261029 clock) -- sse_mulps_addps_forwarding -- GFLOPS @ 2.40GHz: 7.923 [flops/clock] = 19.015 [GFLOPS] (67108864 flops in 8470251 clock = 0.003529 sec) Throughput: 1.981 [instructions/clock] (16777216 instrucions in 8470251 clock)
#include <time.h> #include <sys/time.h> #include <stdint.h> #include <inttypes.h> #include <stdio.h> #include <stdlib.h> #define GHz 2.00 static inline uint64_t rdtsc() { uint64_t ret; #if defined _LP64 asm volatile ( "rdtsc\n\t" "mov $32, %%rdx\n\t" "orq %%rdx, %%rax\n\t" "mov %%rax, %0\n\t" :"=m"(ret) : :"%rax", "%rdx" ); #else asm volatile ( "rdtsc\n\t" "mov %%eax, %0\n\t" "mov %%edx, %1\n\t" :"=m"(((uint32_t*)&ret)[0]), "=m"(((uint32_t*)&ret)[1]) : :"%eax", "%edx" ); #endif return ret; } void print_GFLOPS(double flops, uint64_t cycles) { double GFLOPS = flops * GHz / cycles; double sec = (double)cycles * 1e-9 / GHz; printf("GFLOPS @ %.2fGHz:\n %.3f [flops/clock] = %.3f [GFLOPS] (%.0f flops in %"PRIu64" clock = %f sec)\n", GHz, flops / (double)cycles, GFLOPS, flops, cycles, sec); } void print_throughput(uint64_t instructions, uint64_t cycles) { printf("Throughput:\n %.3f [instructions/clock] (%"PRIu64" instrucions in %"PRIu64" clock)\n", (double)instructions / (double)cycles, instructions, cycles); } #define zero_all_xmm() \ do { \ asm volatile \ ("xorps %xmm0, %xmm0\n\t" \ "xorps %xmm1, %xmm1\n\t" \ "xorps %xmm2, %xmm2\n\t" \ "xorps %xmm3, %xmm3\n\t" \ "xorps %xmm4, %xmm4\n\t" \ "xorps %xmm5, %xmm5\n\t" \ "xorps %xmm6, %xmm6\n\t" \ "xorps %xmm7, %xmm7\n\t" \ ); \ } while (0) #define LOOP (1 << 21) void sse_mulps_addps_forwarding() { const int flops_per_instruction = 4; const int instructions_per_loop = 8; const double flops = flops_per_instruction * instructions_per_loop * LOOP; uint64_t clk0, clk1, cycles; int i; float __attribute__ ((aligned(16))) a[4] = {0.0, 0.0, 0.0, 0.0}; printf("-- sse_mulps_addps_forwarding --\n"); zero_all_xmm(); clk0 = rdtsc(); for (i = 0; i < LOOP; ++i) { asm volatile ( "mulps %0, %%xmm0\n\t" "addps %%xmm0, %%xmm4\n\t" "mulps %0, %%xmm1\n\t" "addps %%xmm1, %%xmm5\n\t" "mulps %0, %%xmm2\n\t" "addps %%xmm2, %%xmm6\n\t" "mulps %0, %%xmm3\n\t" "addps %%xmm3, %%xmm7\n\t" : :"m"(a[0]) ); } clk1 = rdtsc(); cycles = clk1 - clk0; print_GFLOPS(flops, cycles); print_throughput(instructions_per_loop * LOOP, cycles); } void sse_mulps_addps_no_dependency() { const int flops_per_instruction = 4; const int instructions_per_loop = 8; const double flops = flops_per_instruction * instructions_per_loop * LOOP; uint64_t clk0, clk1, cycles; int i; printf("-- sse_mulps_addps_no_dependency --\n"); zero_all_xmm(); clk0 = rdtsc(); for (i = 0; i < LOOP; ++i) { asm volatile ("mulps %xmm0, %xmm1\n\t" "addps %xmm2, %xmm3\n\t" "mulps %xmm4, %xmm5\n\t" "addps %xmm6, %xmm7\n\t" "mulps %xmm0, %xmm1\n\t" "addps %xmm2, %xmm3\n\t" "mulps %xmm4, %xmm5\n\t" "addps %xmm6, %xmm7\n\t" ); } clk1 = rdtsc(); cycles = clk1 - clk0; print_GFLOPS(flops, cycles); print_throughput(instructions_per_loop * LOOP, cycles); } int main() { sse_mulps_addps_no_dependency(); sse_mulps_addps_forwarding(); return 0; }
xmmレジスタが8本な環境に合わせて書きました．
しかし，16本使用可能な環境でも，使用レジスタを増やせば速くなるという訳ではなさそうです．
測定環境についての注意
測定は，RDTSC命令を用いて行っています．
この命令で読み取るカウンタは，CPUの通常動作時の速度で更新されます(自分はそう理解しています)．
従って，CPUの動作周波数を動的に変化させる DVFS (Dynamic Voltage and Frequency Scaling) や TurboBoost が有効になっているマシンだと，正確な測定はできません(8flops/clockを超えるように見えたりする)．
どちらもBIOSでオフにできる機能なはずですので，正確な測定環境を求める方はご検討ください．
