SSE.浮動小数点演算手動最適化は本当に効果的なのか

ちょっと試したので、結果をコピペ。

概要とたてまえ

SSEは、x86のSIMD命令セットに含まれる浮動小数点演算の複数同時実行を行う命令セットです。単精度(float)で4つ、倍精度(double)で2つの演算を1命令で実行できるため、うまく使うと繰り返し同じ演算のパフォーマンスアップを期待できます。
Streaming SIMD Extensions - Wikipedia
画像処理や、確率の計算などでは浮動小数点数の計算を数百万回から回数とかいうレベルではなく3日間くらいの規模で行うことがあるので、

• 少しでも速くなると時間的にとてもうれしい！
• 計算をどこで妥協するかみたいな部分があるため速く計算を終わらせることでやれることの可能性が広がる！

という思いがあります。

CUDAのほうが云々

CUDAは、NVIDIAのGPUが対応しているGPGPUの環境です。GPGPUというのは、グラフィックボードの並列コアを汎用コンピューティングに使う技術で、繰り返し同じ演算を行うような処理で驚きのパフォーマンスアップを期待できます。その実力20〜100倍というウワサです！
CUDA - Wikipedia
これまた夢が広がりますが、

• NVIDIAのグラフィックボードが必要（すごく値段の高いやつがきっといい）
• CUDAドライバのインストールが必要

と実行環境が制限されるため、これがないと動かないソフトは敷居が高くなるため、できればオプション要素にしたいという思いがあります。なので、まずSSE。

コンパイラの最適化で云々

コンパイラが自動でやってくれるとうれしいですが、やってくれるか不安なところがあるのと、はじめから手で書く気があるとループを展開したりSSEが効果的に使えるような構造でプログラムを書くことができます。とりあえずは、手で書いてみたのものと、C言語版(コンパイラの最適化だより)を比較します。

かけ算の結果をひたすら足してみた

float fsum_c(float *input, float *weight, int n)
{
  int i;
  float sum = 0;
  
  for (i = 0; i < n; ++i) {
    sum += input[i] * weight[i];
  }
  
  return sum;
}

2つの配列と要素数を受け取って、両配列の同じ添え字要素の乗算結果の総和を返す関数を試しました。
これをSSE命令セットを使って書くと次のようになります。インラインアセンブラではなく、イントリンシック命令を使っています。イントリンシックは、gccとVC++で共通なので、インラインアセンブリより互換性に優れています。詳しくは、

Compiler Intrinsics (C++)

。

float fsum_sse(float *input, float *weight, int n)
{
  __m128 w, x, u;
  ALIGN16 float mm[4] = {0};
  div_t lp = div(n, 4);
  int i, j, pk_lp = lp.quot, nm_lp = lp.rem;
  float sum = 0;
  
  u = _mm_load_ps(mm);
  
  for (i = 0, j = 0; i < pk_lp; ++i, j += 4) {
    w = _mm_load_ps(&input[j]);  // 4つ分ロード
    x = _mm_load_ps(&weight[j]); // 4つ分ロード
    x = _mm_mul_ps(w, x); // 4つ分乗算
    u = _mm_add_ps(u, x); // 4つ分uに加算
  }
  _mm_store_ps(mm, u); // メモリに書き戻す
  
  sum = mm[0] + mm[1] + mm[2] + mm[3]; // 小計
  
  for (i = 0, j = pk_lp * 4; i < nm_lp; ++i, ++j) {
    sum += input[j] * weight[j]; // 残り分
  }
  
  return sum;
}

単精度では、4つ同時に計算できるため、4の倍数単位でループをまわして、4つずつ計算を行います。余り分は最後に計算して加算します。

結果

上の単精度(float)バージョンと、倍精度(double)バージョンを配列要素数10,003個で、100,000回行う間の時間を計測して結果としました。また、OpenMPによる並列化との組み合わせも行いました。実行環境のCPUは、Intel Core 2 Quad Q9550(2.83GHz x 4コア)です。

// モード: 経過時間(ms), テスト用の値

// CバージョンとSSEバージョン比較
- float
  C: 3594, 2.490360E+003
SSE: 562, 2.490360E+003
6.40 倍高速化!

- double
  C: 1344, 2.473939E+003
SSE: 1078, 2.473939E+003
1.25 倍高速化!

// CバージョンとSSE+OpenMPバージョン比較
- float
  C: 3578, 2.490360E+003
SSE: 156, 2.491235E+003
22.94 倍高速化!

- double
  C: 1344, 2.473939E+003
SSE: 281, 2.473939E+003
4.78 倍高速化!

結果だけ見ると、単精度すごい、倍精度あまり変わらないといった印象です。

kousatsu

単精度では、かなりの効果があるように見えますが、VCの単精度計算は、キャストしているとかで遅いという話があるので、SSEバージョンが速いというよりは、Cバージョンが遅いという疑いがあります（つまりSSE速い！）。doubleのほうは、最適化なしだとかなり遅かったのでコンパイラの最適化でSSEが使われてるようです。OpenMPも含めたバージョンではかなり効果が出ていますが、倍精度はほぼOpenMPによる効果だということと、単精度は演算手順が変わるためか丸め誤差が増幅されてCバージョンとかなり違う値になってしまいました(doubleバージョンはズレないので、並列化がバグっているわけではないと思う）。

コメント

ディスアセンブルしないと正確な状況など分かるわけない、でも／(＾o＾)＼って感じはしている。
丸め誤差はFPUのモードによると思うので、調べる。
コンパイラが最適化しなさそうな、もっと複雑な演算の手動での展開やパックドバイトでの文字列処理をやらせたほうが効果がありそう。

今回試したコード

// cl sum.c vcomp.lib /arch:SSE2 /O2 /openmp

#include <math.h>
#include <emmintrin.h>
#include <xmmintrin.h>
#include <windows.h> // GetTickCount
#include <omp.h>

#ifdef _linux // gcc
#define ALIGN16 __attribute__((aligned(16)))
#else // cl
#define ALIGN16 _declspec(align(16))
#endif

float fsum_c(float *input, float *weight, int n)
{
  int i;
  float sum = 0;
  
  for (i = 0; i < n; ++i) {
    sum += input[i] * weight[i];
  }
  
  return sum;
}

float fsum_sse(float *input, float *weight, int n)
{
  __m128 w, x, u;
  ALIGN16 float mm[4] = {0};
  div_t lp = div(n, 4);
  int i, j, pk_lp = lp.quot, nm_lp = lp.rem;
  float sum = 0;
  
  u = _mm_load_ps(mm);
  
  for (i = 0, j = 0; i < pk_lp; ++i, j += 4) {
    w = _mm_load_ps(&input[j]);
    x = _mm_load_ps(&weight[j]);
    x = _mm_mul_ps(w, x);
    u = _mm_add_ps(u, x);
  }
  _mm_store_ps(mm, u);
  
  sum = mm[0] + mm[1] + mm[2] + mm[3];
  
  for (i = 0, j = pk_lp * 4; i < nm_lp; ++i, ++j) {
    sum += input[j] * weight[j];
  }
  
  return sum;
}

double dsum_c(double *input, double *weight, int n)
{
  int i;
  double sum = 0;
  
  for (i = 0; i < n; ++i) {
    sum += input[i] * weight[i];
  }
  
  return sum;
}

double dsum_sse(double *input, double *weight, int n)
{
  __m128d w, x, u;
  ALIGN16 double mm[2] = {0};
  div_t lp = div(n, 2);
  int i, j, pk_lp = lp.quot, nm_lp = lp.rem;
  double sum = 0;
  
  u = _mm_load_pd(mm);
  
  for (i = 0, j = 0; i < pk_lp; ++i, j += 2) {
    w = _mm_load_pd(&input[j]);
    x = _mm_load_pd(&weight[j]);
    x = _mm_mul_pd(w, x);
    u = _mm_add_pd(u, x);
  }
  _mm_store_pd(mm, u);
  
  sum = mm[0] + mm[1];
  
  for (i = 0, j = pk_lp * 2; i < nm_lp; ++i, ++j) {
    sum += input[j] * weight[j];
  }
  
  return sum;
}

#define TC 100000
#define DC 10003
int main(void)
{
  ALIGN16 float fdata1[DC], fdata2[DC];
  ALIGN16 double ddata1[DC], ddata2[DC];
  float fsum = 0;
  double dsum = 0;
  int i, j;
  long st, t1, t2;

#ifdef _OPENMP
  omp_set_num_threads(omp_get_num_procs() * 8);
#endif
  
  //srand(time());
  
  for (j = 0; j < DC; ++j) {
    fdata1[j] = (float)(rand() % 1000) / 1000;
    fdata2[j] = (float)(rand() % 1000) / 1000;
    ddata1[j] = (double)(rand() % 1000) / 1000;
    ddata2[j] = (double)(rand() % 1000) / 1000;
  }
  
  /////// 単精度
  puts("- float");
  
  // normal C
  fsum = 0.0;
  st = GetTickCount();
  for (i = 0; i < TC; ++i) {
    fdata1[0] = 1.0;
    fsum += fsum_c(fdata1, fdata2, DC) / TC;
  }
  t1 = GetTickCount() - st;
  printf("%3s: %ld, %E\n", "C", t1, fsum);
  
  // SSE
  fsum = 0.0;
  st = GetTickCount();
#ifdef _OPENMP
#pragma omp parallel for reduction(+:fsum) 
#endif
  for (i = 0; i < TC; ++i) {
    fdata1[0] = 1.0;
    fsum += fsum_sse(fdata1, fdata2, DC) / TC;
  }

  t2 = GetTickCount() - st;
  printf("%3s: %ld, %E\n", "SSE", t2, fsum);
  printf("%0.2f 倍高速化!\n\n", (double)t1 / t2);

  
  /////// 倍精度
  puts("- double");
  
  // normal C
  dsum = 0.0;
  st = GetTickCount();
  for (i = 0; i < TC; ++i) {
    ddata1[0] = 1.0;
    dsum += dsum_c(ddata1, ddata2, DC) / TC;
  }
  t1 = GetTickCount() - st;
  printf("%3s: %ld, %E\n", "C", t1, dsum);
  
  // SSE
  dsum = 0.0;
  st = GetTickCount();
#ifdef _OPENMP
#pragma omp parallel for reduction(+:dsum)
#endif
  for (i = 0; i < TC; ++i) {
    ddata1[0] = 1.0;
    dsum += dsum_sse(ddata1, ddata2, DC) / TC;
  }
  t2 = GetTickCount() - st;
  printf("%3s: %ld, %E\n", "SSE", t2, dsum);
  printf("%0.2f 倍高速化!\n\n", (double)t1 / t2);
  
  return 0;
}