読者です 読者をやめる 読者になる 読者になる

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でオフにできる機能なはずですので,正確な測定環境を求める方はご検討ください.

*1:この記事と同じようにやれば4flops/clockの倍精度浮動小数点演算も可能だと思います.未確認

*2:このことは, 『インテル 64 アーキテクチャーおよびIA-32アーキテクチャー最適化リファレンス・マニュアル』や,『インテル® アーキテクチャ(IA)浮動小数点ユニット(FPU)、ストリーミングSIMD 拡張命令(SSE)、ストリーミングSIMD 拡張命令2(SSE2)を使用した浮動小数点算術演算』にも明示的に書かれていません