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でオフにできる機能なはずですので,正確な測定環境を求める方はご検討ください.
gccで指定できるアラインメントのサイズには限界がある
C言語におけるデータのアラインメントについての小ネタです.
アラインメント自体については,
データ型のアラインメントとは何か,なぜ必要なのか?
にとても詳しいです.
問題設定
4バイトの変数を16個の配列で確保するとします.
int a[16];
もし各レベルキャッシュのキャッシュラインサイズが64バイトなら,
この配列aを1ラインに収めることができるはずです.
事例1: 配列サイズを定数で指定
__atribute__ ((aligned(64))) を使って,アラインメントを64バイト境界に指定してみましょう.
alignment_stack_const.c
#include <stdio.h> #include <stdint.h> int main() { uintptr_t addr; int __attribute__ ((aligned(64))) a[16]; addr = (uintptr_t)a; printf("a = %p\n", (void *)a); printf("%p %% 64 == %d\n", (void *)a, addr % 64); return 0; }
(ポインタを整数に変換する部分は Cのポインタを整数に変換する - bkブログ を参照)
これをコンパイルして実行すると,aは確かに64バイト境界にアラインされていることが分かります.
$ gcc alignment_stack_const.c $ ./a.out a = 0xbffd3c40 0xbffd3c40 % 64 == 0 $ ./a.out a = 0xbf80ef40 0xbf80ef40 % 64 == 0 $ ./a.out a = 0xbfa08640 0xbfa08640 % 64 == 0 ...
事例2: 配列サイズを変数で指定
比較的新しいCの規約では,配列長を変数で指定することができます.
この機能を使って配列を確保してみます.
alignment_stack_var.c
#include <stdio.h> #include <stdint.h> int main() { uintptr_t addr; int n = 16; int __attribute__ ((aligned(64))) a[n]; addr = (uintptr_t)a; printf("a = %p\n", (void *)a); printf("%p %% 64 == %d\n", (void *)a, addr % 64); return 0; }
先程のソースとは,「16」という配列長を「n」という変数を介して指定している部分以外変わりませんが,
これをコンパイルして実行すると,64バイト境界にアラインメントがとれない場合があります.
*1
*2
$ gcc alignment_stack_var.c $ ./a.out a = 0xbf8e5130 0xbf8e5130 % 64 == 48 $ ./a.out a = 0xbfb98d60 0xbfb98d60 % 64 == 32 $ ./a.out a = 0xbfd47d80 0xbfd47d80 % 64 == 0 ...
何に気をつければよいか
実はアラインメントのサイズには,コンパイラ・リンカが定める最大値があります.
gcc の定義済みマクロ, "__BIGGEST_ALIGNMENT__" を確認することで,
お使いの環境での最大値が分かります.
$ cpp -dD alignment_stack_var.c |grep __BIGGEST_ALIGNMENT__ #define __BIGGEST_ALIGNMENT__ 16
自分の環境では __BIGGEST_ALIGNMENT__ の値が 16 になっていたため,
いつでも64バイトのアラインメントがとれることを期待してはいけないようでした.
一方で, alignment_stack_var.c をコンパイルして実行すると,最低でも16バイト境界にアライン
されていることは確認できるはずです.
perf statでL1,L2(,L3)キャッシュミス測定
この記事は↓に移転しました。
laysakura.github.io
perf stat でパフォーマンスカウンタの値を直接指定し表示
この記事は↓に移転しました。
laysakura.github.io
第六回 カーネル/VM探検隊 の発表資料です
先ほど,第六回 カーネル/VM探検隊 にてLTをさせていただきました.
Google Summer of Codeで取り組むことについてです.
本日発表に使用した資料を公開しますので,コメントなりTwitterなりで質問歓迎です.
アニメーションなしだと判りやすさ(とか)が激減しますので,こちらの OpenOffice Impress で作成したファイルを強く推奨します.
非推奨PDFバージョン(アニメーションなし)
10分のLTということもあり,話しきれなかった部分や説明の不足している部分も多々あります.
近々,やろうとしていることをもう少し詳しく公開する機会(Kernel/VMではなくw)があると思いますので,その時にはこのblogにも公開しますね.