SIMD

職業 C++ プログラマの間では割と常識っぽいけど web 上に日本語の情報があんまりない気がする話を初心者向けに解説してみる試み。想定読者はゲーム業界入る前の自分。
今回は SIMD について。

SIMD / SSE

SIMD とは "Single Instruction Multiple Data" の略で、文字通り 1 命令で複数のデータを処理することです。これを使うと、例えば 4 要素のベクトル同士の各種演算を 1 命令で行うことができ、高速化を図れます。PC では PentiumIII 以降の CPU が SIMD の一種 (SSE) を実装しており、近年の CPU では SSE も色々改良されています。
SSE は各種コンパイラの拡張命令 (intrinsics) を使うことで、C++ から利用することができます。今記事では VisualC++/gcc での SSE の使い方とかについて適当に解説してみます。

簡単な使用例

// ssetest.cc
#include <cstdio>
#include <xmmintrin.h>

int main()
{
    __m128 a, b, c;
    a = _mm_setr_ps(1.0f, 2.0f, 3.0f, 4.0f);
    b = _mm_setr_ps(5.0f, 6.0f, 7.0f, 8.0f);
    c = _mm_add_ps(a, b);
    c = _mm_mul_ps(c, a);
    c = _mm_sub_ps(c, b);

    float* r = (float*)&c;
    printf("%.2f, %.2f, %.2f, %.2f\n", r[0], r[1], r[2], r[3]);
}

SSE を使うには xmmintrin.h を include します。__m128 という型が SSE の命令に使われるデータ型 (レジスタに直接収まる) で、16byte にアライメントされた 16byte (128bit) のデータです。あとは見たまんまで、各種 intrinsics を使って 4 要素の float を足したり掛けたり引いたりしています。
コンパイルには、gcc の場合 -msse オプションが必要です。

ベクトル計算への応用

ゲームプログラムの場合、SIMD が最も活躍するのがベクトルや行列などの数学関係だと思われます。SSE を使用して簡単なベクトルクラスを実装してみます。

// ssevec4.cc
#include <cstdio>
#include <xmmintrin.h>

struct vec4_t
{
    union {
        float fv[4];
        struct { float x,y,z,w; };
    };

    vec4_t() { for(int i=0; i<4; ++i){ fv[i]=0.0f; } }
    vec4_t(__m128 v) { *this=v; }
    vec4_t(const float* v) { for(int i=0; i<4; ++i){ fv[i]=v[i]; } }
    vec4_t(float _1, float _2, float _3, float _4) { x=_1; y=_2; z=_3; w=_4; }

    operator __m128() { return reinterpret_cast<__m128&>(*this); }
    operator const __m128() const { return reinterpret_cast<const __m128&>(*this); }
    vec4_t& operator=(__m128 v) { reinterpret_cast<__m128&>(*this)=v; return *this; }
    float& operator[](int i) { return fv[i]; }
    const float& operator[](int i) const { return fv[i]; }

    __m128 operator+(__m128 v) const { return _mm_add_ps(*this, v); }
    __m128 operator-(__m128 v) const { return _mm_sub_ps(*this, v); }
    __m128 operator*(__m128 v) const { return _mm_mul_ps(*this, v); }
    __m128 operator/(__m128 v) const { return _mm_div_ps(*this, v); }
};

// VisualC++ と gcc でアライン指定方法が違うので分岐…
#ifdef _MSC_VER
    #define ALIGN(n) __declspec(align(n))
#else
    #define ALIGN(n) __attribute__((aligned(n)))
#endif
typedef ALIGN(16) vec4_t vec4al_t;

int main()
{
    vec4al_t a(1.0f, 2.0f, 3.0f, 4.0f);
    vec4al_t b(5.0f, 6.0f, 7.0f, 8.0f);
    vec4al_t c = a + b;
    c = c * a;
    c = c - b;
    printf("%.2f, %.2f, %.2f, %.2f\n", c[0], c[1], c[2], c[3]);
}

__m128 をそのままメンバとして持ってしまうと問題が出ることがあります。例えば、std::vector に詰め込もうとすると VisualC++ ではコンパイルエラーが出ます。(アライメントの都合？) この対策にトリッキーなことをやっています。
まず本体は __m128 との互換性を持たせつつアライメント指定なしで定義しておき、後でアライメント指定があるバージョンを typedef します。通常はアライメント指定版を使い、std::vector に詰めたりする場合だけアライメント無指定版を使います。同じ型なので代入は相互にできます。ちなみに、アライメントが 16byte に揃ってないデータで SSE 演算しようとすると遠慮無くクラッシュします。注意が必要です。
また、上記のソースで operator+-*/ の引数と返り値は vec4_t にしたかったのですが、そうすると gcc で遅いコードになってしまったため、仕方なく __m128 でやりとりするようにしています…。

注意点

SSE を使う処理と各要素へのアクセスが交互に頻繁に発生するようなケースでは遅くなります。これは、要素へのアクセスは SSE のレジスタ→メモリ→FPU レジスタ またはその逆のデータの移動が発生するためです。メモリを介する操作は遅いので、SSE のレジスタで完結できる処理はがんばって SSE の命令でデータをこね回してなんとかするのが高速化につながるようです。(要素の交換 (_mm_shuffle_ps() など)、単一要素の計算 (_mm_add_ss() など後ろが "ss" の命令) などを組み合わせてなんとかする。)

プリフェッチ (_mm_prefetch()) でメモリアクセスを最適化する、ストールに備えてできるだけ並列に処理できる構成にしておく (直前の計算結果を次の計算に使わないとか) 、などの注意点もあるようですが、ここらへんはまだ私の理解が浅いので今回は解説を避けます。

コンパイラによる自動ベクトル化など

gcc の 4.2 系以降 と icc は適切な記述をすればコンパイラが自動的に SIMD 化 (ベクトル化) してくれるようです。例えば以下のようなソース。

// vectorize.cc
#include <cstdio>

int main()
{
    float a[4] = {1.0f, 2.0f, 3.0f, 4.0f};
    float b[4] = {5.0f, 6.0f, 7.0f, 8.0f};
    float c[4];

    for(int i=0; i<4; ++i) { c[i] = a[i] + b[i]; }
    for(int i=0; i<4; ++i) { c[i] = c[i] * a[i]; }
    for(int i=0; i<4; ++i) { c[i] = c[i] - b[i]; }

    printf("%.2f, %.2f, %.2f, %.2f\n", c[0], c[1], c[2], c[3]);
}

g++ vectorize.cc -O3 -msse -S

でコンパイルすると、for 3 連打の付近は以下のコードに化けます。(自動ベクトル化のオプションは -ftree-vectorize で、-O3 には同じ効果が含まれます)

    movaps    80(%esp), %xmm0
    addps    64(%esp), %xmm0
    mulps    80(%esp), %xmm0
    subps    64(%esp), %xmm0

レジスタをうまく使ってくれてはいないものの、SSE の命令に置き換えられています。固定長ループがベクトル化されるようです。可変長ループやアンロールしたループの場合はベクトル化されませんでした。また、対象データが構造体のメンバだった場合とかには適切なアライメント指定を付けておいた方がいいと思います。無い場合余計な命令が増えました。

VisualC++ は自動ベクトル化には未対応のようですが、"/arch:SSE /fp:fast" オプションを指定することにより浮動小数点演算を SSE で行うようになり、Load&Store を減らせるようです。

Intel が近い先に出す CPU (Sandy Bridge) では 256bit 幅の SIMD (AVX) が実装され、以降 512bit 幅や 1024bit 幅の SIMD に対応した CPU を出していく予定らしく、たぶんそのへんまで行くと人力で SIMD 化頑張るのは限界なんじゃないかと思われるので、SIMD 化は基本コンパイラに頑張ってもらうようになっていくのかもしれません。

まとめ

現状 SSE 使おうとするとアセンブラみたいな intrinsics を使わないといけなかったり、アライメントやレジスタを意識しないといけなくなったり、急に前時代的な世界に突入します。SIMD 化にしろメニーコア化にしろ、今は過渡期なので、10 年くらい先にはベクトル化もタスク化もコンパイラがそれなりによろしくやってくれる世界になっているかもしれません。しかし現状まだそうなってないので、シビアに速度を追い求める場合は泥に塗れてがんばりましょう。

参考

http://msdn.microsoft.com/ja-jp/library/t467de55.aspx MSDN の SSE のリファレンス
http://d.hatena.ne.jp/gamesyokunin/20041120 アライメント指定版/指定なし版で型を分けるテクニックはこの記事を参考にさせていただきました
http://gcc.gnu.org/projects/tree-ssa/vectorization.html gcc の自動ベクトル化機能について
http://wlog.flatlib.jp/item/941 VisualC++ の /arch:SSE /fp:fast オプションなどについて
http://www.4gamer.net/games/105/G010549/20100831055/ CEDEC2010 の Intel の講演の記事。将来のコンパイラと CPU についてなど