«前の日記(2009年02月02日) 最新 次の日記(2009年02月05日)» 編集

ema log


2009年02月03日 [長年日記]

_ [Programming]SIMD演算

体調崩して寝込んでました.そんなこんなで,たまには,プログラムネタを書いておかないとという気分になったので.図が無くて,文字ばっかでごめんなさい.

最近のCPUにはSIMD演算と呼ばれる命令群が搭載されています.たとえば,x86ならSSE,PowerPCならAltiVecなんかがそうです.SIMD*1 というのは,Single Instruction, Multiple Dataの略で,複数データに同じ演算を施すことを意味します.

コンピュータアーキテクチャ 定量的アプローチ 第4版 (IT Architects’ Archive)(John L. Hennessy/ジョン・L・ヘネシー/デイビッド・A・パターソン/David A. Patterson/中條 拓伯/吉瀬 謙二/佐藤 寿倫/天野 英晴)

これらの世界では,データはスカラとベクトルに区別されます. 例えば,SSEやAltiVecでは128bitのレジスタを8〜64bit単位に区切ってベクトルとしています. これらでは,128bitを,32bit整数×4や,32bit浮動小数点×4として扱います.

それの何が嬉しいかというと,一番わかりやすいのは,同時に働く演算器の数が増えることです.例えば,x86 のとある世代の FADD という80bit浮動小数点の加算命令はレイテンシーが6,スループットが1ですが,32bit浮動小数点4つ同士を加算するADDPS(おそらく,ADD Packed Single float の意味)のレイテンシーは5,スループットは2となっています[2].bit数の違いはありますが,浮動小数点の演算器を並列に並べているので,4つの計算を一度に並列に行えるわけです*2

ということで,SIMD演算を利用する際には,(32bitを使うことが多いでしょうから)データを4つずつに区切っての演算が基本になります.特に,SSEではベクトル同士の演算(垂直演算)が中心で,ベクトル内の演算(水平演算)はSSE2ではほとんど存在しません*3.小難しく書くと,要はベクトル同士が独立していることを前提としています.画像で書くと,ある4画素値の計算に,手前の1画素が必要だったりすると面倒なことになるわけです.また,RGBRGBとバイト単位で並んでいるよりは,RRRRGGGGBBBBと各画素が独立している方が扱いやすくなります.演算本体よりもそのデータの用意が面倒ですw

実際の,SSEの利用には,「IA-32 SIMDの扉」というページがわかりやすいです.「ネガティブ」では最も単純な4画素値ずつの反転処理が,「モノクロ変換」では,演算のしやすさを考えて並び替えを行う例が示されています.

ここではインラインアセンブラが使用されていますが,最近のコンパイラ(gcc, icc, Visual C++のいずれでも利用可能)には intrinsic [3] と呼ばれるものが備わっていて,Cの関数呼び出しのようにSSEを利用できます*4intrinsic関数はVCとgccとで共通化され,しかも32bit, 64bit両対応なため汎用性も高いそうです[4]し,どうやらコンパイラの最適化の対象にもなります[5].使い方は「SSE についてのメモ(2) SSE4など」や「SSE2 performance optimization.」などがわかりやすいです.

例えば,SSE2 の利用であれば,emmintrin.h を読み込んで,以下のように計算します.結果は 2, 4, 6, 8 のように出力されます.データの並び順に注意して下さい.

#include <stdio.h>
#include <emmintrin.h>

int main(void) {
  __m128i a,b,c; // 実際には xmm のどれかに割り当てられる.
  int i, dst[4];
  int src[4] = { 1, 2, 3, 4 };

  // データの用意.
  a = _mm_loadu_si128( (__m128i*)src ); // メモリから4つまとめて読み出す
  b = _mm_set_epi32( 4, 3, 2, 1 ); // マクロでレジスタを初期化.↑との並び順の違いに注意.逆アセすると偉いことにw

  // 4つまとめて加算する.4.times {|i| c[i] = a[i] + b[i] } と等価
  c = _mm_add_epi32( a, b );

  // 結果表示用に,メモリに4つまとめて書き戻す
  _mm_storeu_si128( (__m128i*)dst, c );

  for(i=0; i<4; i++) {
    printf("%d, ",dst[i]);
  }
  printf("\n");
}

gcc でコンパイルする場合には,-msee2 が必要です.gcc 4.2.4 で動作確認しました.dst については確実にアラインして,_mm_load_si128, _mm_store_si128 を用いる方が性能面では望ましいです*5

実際に,どのようにコンパイルされるかをみてみると,(ここでは,わかりやすさのために,アラインして,最適化オプションは -O1 -msse2 にして,_mm_set_epi32 を使っていません.O3 にすると,命令が並び変わっていました)

#include <stdio.h>
#include <emmintrin.h>

int main(void) {
  __m128i a,b,c;
  int i;
  int dst[4] __attribute__((aligned(16))) ;
  int src1[4] __attribute__((aligned(16))) = {1,2,3,4};
  int src2[4] __attribute__((aligned(16))) = {2,3,4,5};

  a = _mm_load_si128( (__m128i*)src1 );
  b = _mm_load_si128( (__m128i*)src2 );
  c = _mm_add_epi32(a,b);
  _mm_store_si128( (__m128i*)dst, c );

  for(i=0; i<4; i++) {
  	printf("%d, ",dst[i]);
  }
  printf("\n");
}

これを,コンパイルして objdump -d すると

80483b7:       c7 45 c8 01 00 00 00    movl   $0x1,-0x38(%ebp) // src1 の用意
80483be:       c7 45 cc 02 00 00 00    movl   $0x2,-0x34(%ebp)
80483c5:       c7 45 d0 03 00 00 00    movl   $0x3,-0x30(%ebp)
80483cc:       c7 45 d4 04 00 00 00    movl   $0x4,-0x2c(%ebp)
80483d3:       c7 45 b8 02 00 00 00    movl   $0x2,-0x48(%ebp) // src2 の用意
80483da:       c7 45 bc 03 00 00 00    movl   $0x3,-0x44(%ebp)
80483e1:       c7 45 c0 04 00 00 00    movl   $0x4,-0x40(%ebp)
80483e8:       c7 45 c4 05 00 00 00    movl   $0x5,-0x3c(%ebp)
80483ef:       66 0f 6f 45 c8          movdqa -0x38(%ebp),%xmm0 // src1 の読み出し
80483f4:       66 0f fe 45 b8          paddd  -0x48(%ebp),%xmm0 // src1 と src2 の加算
80483f9:       66 0f 7f 45 d8          movdqa %xmm0,-0x28(%ebp) // dst への保存

のようになります.x86ではメモリをオペランドに指定できるので,movdqa は一つしかないですが,intrinsic と SSE 命令が1対1に対応しながらも,コンパイラによる最適化の恩恵にも与れることが確認できました*6

ほかにも,メモリアクセスなどでも有利だったり,アラインしてないとロードストアがなどなど,メモリの問題がありますが,今回は一番基本的な演算部分だけおおざっぱに書いています.SSE化より,プリフェッチ命令の方が効いたりすることもあるでしょう.

References
[1] コンピュータアーキテクチャ 定量的アプローチ 第4版, p.211
[2] IA-32 インテル(R) アーキテクチャー最適化リファレンス・マニュアル, pp. C-11-C-14 (PDF)
[3] MMX, SSE, and SSE2 Intrinsics (C++)
[4] fast strlen and memchr by SSE2 (mitsunari@cybozu labs)
[5] 本の虫: Compiler Intrinsicの吐くコード
2.4 計算しやすいベクタデータの作成 - PS3 Linux Information Site / Cell/B.E.のパワーを体験しよう
IA-32 SIMDの扉
Compiler Intrinsics
SSE についてのメモ(2) SSE4など
SSE2 performance optimization.
【最新パーツ性能チェック(Vol.23)】いよいよプレスコット登場(PART2)! SSE3の神髄に世界で初めて触れる!

*1 読み方は日本ではシムドが一般的ですが,「へネパタ本」[1] の訳注では「シムディ」と発音して欲しいと書かれています.

*2 (たぶん.というのも,世代によっては整数は64bitの MMX ユニット2回で計算とかしてたような?それでも,普通にやるより速いけど)

*3 例えば,SSE2で,ベクトル内の総和をとることは大変です.いずれにせよ,SIMD演算に合わせてデータを用意しておかなければ大変なのですが・・・

*4 SSE使うような人たちはインラインアセンブラでごりごり書く人だからか,instrinsic の情報が少ない・・・

*5 gcc だと __attribute__((aligned(16))) をつけますが,VC などでは方式が違ったはずなので,ここでは_mm_storeu_si128としています.アラインされていないと,movdquなどでは,境界をまたぐデータの読み出しに複数データを読み出して合成するため遅くなります.

*6 paddd は Packed ADD Double word integer の略だと思われます.歴史的経緯により,x86 で word というと 16bit で,32bit なら double word, 64bitなら quadruple word となります.