2008年02月23日

はじめに

イントリンシックで MMX や SSE が使えるようになり、手軽に SIMD が利用できるようになった。
でも、いざ使おうとするとどう書けばいいのかわからない場面に出くわすことが多い。
だからと言って、一度ばらして普通に処理すると遅くなってしまう。

ここでは SIMD 化するに当たって直面した問題を解決した方法や他のソースコードから学んだテクニック、検索して見つけたものなどをまとめて行こうと思う。

投稿者 Takenori : 20:48 | トラックバック

比較

SIMD の比較は、単にマスクを生成してくれるだけだ。
比較して分岐するなどは出来ない。
生成されたマスクを使って何とかするしかない。
これはいつも分岐で書いていると、どう書けばいいのかさっぱりわからない問題だ。
でも、これが使えると分岐を消して並列で処理できるようになる。

比較的わかりやすいであろう、theora のフィルタのソースコードの一部を以下に書いた。
MMX を使っている。

まずは条件に一致するものを 0 にするコード

// if( R <= -2L || R >= 2L ) R = 0;
mask_L2 = _mm_cmpgt_pi16( limit_x2_p, R ); // 2L >   R ? FF : 00
R = _mm_and_si64( R, mask_L2 );            // 2L >   R ? R : 0
mask_mL2 = _mm_cmpgt_pi16( R, limit_x2_m );//  R > -2L ? FF : 00
R = _mm_and_si64( R, mask_mL2 );           //  R > -2L ? R : 0

マスクを作って and をとって 0 にしているので、if 文の条件が逆転しているのさえわかれば、特に難しいことはない。

次に条件に一致する時に加算するコード。

// if( R < -L ) R2 = R + 2L;
mask_mL_R = _mm_cmpgt_pi16( limit_m, R );      // -L > R ? FF : 00 ( マスクを作る )
tmp_2L = _mm_and_si64( limit_x2_p, mask_mL_R );// -L > R ? 2L : 0
R2 = _mm_add_pi16( R, tmp_2L );                // -L > R ? R + 2L : R

条件に一致する加算する方をマスクし、加算されないものを 0 にしてから加算する。
こうすることで条件に一致するものだけ加算できる

だいたいこの2つがわかれば、比較は他にも応用できる。

追記:
指摘された内容を以下に加筆。
上の判定は、 R + 32767-2L が 32767- 4L 以下ならば R <= -2L || R >= 2L の条件に当てはまるということになる。
値域を整数の上限にまで移動させてしまえば範囲外の比較が一回で済む。
整数の範囲内の計算で範囲がほぼ固定の場合は 値域を整数の上限か下限にもっていってしまえば 比較一回で範囲外か範囲内かがわかるようになります。
だから、値域を移動させるための加算が一回、比較が一回、マスクが一回の形。
計3命令で実現できるかな。

投稿者 Takenori : 20:49 | トラックバック

並列スワップ

SIMD の比較を使うことで、スワップも一度に行うことが出来る。
複数個同時にスワップが出来、分岐がなくなるので高速化に役立つ。
以下、SSE を使ったコード。

// if( x1 > x2 ) swap( x1, s2 ); を行う
__m128 swapmask = _mm_cmpgt_ps( x1, x2 );
__m128 tmp = x1;
x1 = _mm_add_ps( _mm_andnot_ps( swapmask, x1 ), _mm_and_ps( swapmask, x2 ) );  // x1 = (~mask & x1 ) + (mask & x2 );
x2 = _mm_add_ps( _mm_andnot_ps( swapmask, x2 ), _mm_and_ps( swapmask, tmp ) ); // x2 = (~mask & x2 ) + (mask & tmp );

比較を使ってマスクを作り出し、入れ替え対象のものを 0 にして、代入する値を逆のマスクで交換する値のみ残し、加算している。
MMX や SSE2 などでも同じようにしてできる。
整数の時は加算よりも or を取る方が自然かもしれない。

追記:
上の処理だと、min、max を使えばそれで出来ると指摘されたので追記。
コードは以下。

__m128 tmp = x1;
x1 = _mm_min_ps(x1, x2);
x2 = _mm_max_ps(tmp, x2);

上のソースを実際に使っているところでは、この x1 と x2 条件で uv 値もスワップしているため、min、max では実現できない。
例文のソースコードは比較条件を別のものにした方がよかったかも。

さらに追記:
xor でスワップするテクニックを使えば、さらに速くなると指摘を受けたので追記。
コードは以下。

__m128 tmp = _mm_and_ps(_mm_xor_ps(x1, x2), _mm_cmpgt_ps( x1, x2 ));
x1 = _mm_xor_ps(tmp, x1);
x2 = _mm_xor_ps(tmp, x2);

xor でスワップは最近使っていなかったのと、浮動小数点だったので抜け落ちていた。

投稿者 Takenori : 20:52 | トラックバック

クリッピング

SIMD の比較を使えばクリッピングも並列に出来る。
以下の SSE を使ったコードでは、クリッピング以外に X 値がクリッピングされた分、UV 座標値を移動する処理も同時に行っている。

// if( x1 < 0 ) x1 = 0;
__m128 mx1_mask = _mm_cmplt_ps( x1, zero );// x1 < 0.0
__m128 minus_x1 = _mm_sub_ps( zero, x1 );  // -x1
u1 = _mm_add_ps( u1, _mm_and_ps( _mm_mul_ps( du, minus_x1 ), mx1_mask ) ); // u1 += du * -x1;
v1 = _mm_add_ps( v1, _mm_and_ps( _mm_mul_ps( dv, minus_x1 ), mx1_mask ) ); // v1 += dv * -x1;
x1 = _mm_andnot_ps( mx1_mask, x1 ); // x1 < 0.0 のものを0に、~a0 & b0
// if( x2 > dest_w ) x2 = dest_w;
__m128 mw_mask = _mm_cmpgt_ps( x2, dest_width );
x2 = _mm_add_ps( _mm_andnot_ps( mw_mask, x2 ), _mm_and_ps( dest_width, mw_mask ) ); // if( x2 > dest_w ) x2 = dest_w;

普通。

追記:
x2 側は、以下のように min を使うことで実現できる。
同様に x1 側も x1 だけであれば、max を使うことで実現できる。
ただ、上の例では、uv もその判定結果を使っているので、x1 は max ではなく比較を使った方が良い。

x2 = _mm_min_ps( x2, dest_width ); // dest_width 以下にする

投稿者 Takenori : 20:54 | トラックバック

行列の転置

SIMD の処理は垂直演算を基本としているため、行列の転置が必要になることはよくある。
IDCT などではまず使うだろう。
ただ、行列の転置はそれなりに時間のかかる処理だと思うので、転置後の演算が少なくて、また元に戻さないといけない場合、重くなってしまうかもしれない。
転置後の演算が少ない場合は、他の方法で処理できないか一度考えてみる必要がある。

以下、MMX で 8x8 行列の転置を行うコード。
かなり最適化任せの力業。
最適化に頼らずに実装しようとしたら、いろいろと挟み込んで順番入れ替えるなどする必要があって大変。

// このようにソースを読み込んでいるとする。
__m64 row0 = src[0]; // 00 01 02 03 04 05 06 07
__m64 row1 = src[1]; // 10 11 12 13 14 15 16 17
__m64 row2 = src[2]; // 20 21 22 23 24 25 26 27
__m64 row3 = src[3]; // 30 31 32 33 34 35 36 37
__m64 row4 = src[4]; // 40 41 42 43 44 45 46 47
__m64 row5 = src[5]; // 50 51 52 53 54 55 56 57
__m64 row6 = src[6]; // 60 61 62 63 64 65 66 67
__m64 row7 = src[7]; // 70 71 72 73 74 75 76 77

// 係数行列の転置を行う関数
inline void transpose8_8x8_mmx(
 __m64  row0, __m64  row1, __m64  row2, __m64  row3, __m64  row4, __m64  row5, __m64  row6, __m64 row7,
 __m64& col0, __m64& col1, __m64& col2, __m64& col3, __m64& col4, __m64& col5, __m64& col6, __m64& col7 )
{
  // 係数転置 ( フェーズ 1 )
  __m64 tmp01L = _mm_unpacklo_pi8( row0, row1 ); // 00 10 01 11 02 12 03 13
  __m64 tmp01H = _mm_unpackhi_pi8( row0, row1 ); // 04 14 05 15 06 16 07 17
  __m64 tmp67L = _mm_unpacklo_pi8( row6, row7 ); // 60 70 61 71 62 72 63 73
  __m64 tmp67H = _mm_unpackhi_pi8( row6, row7 ); // 64 74 65 75 66 76 67 77
  __m64 tmp23L = _mm_unpacklo_pi8( row2, row3 ); // 20 30 21 31 22 32 23 33
  __m64 tmp23H = _mm_unpackhi_pi8( row2, row3 ); // 24 34 25 35 26 36 27 37
  __m64 tmp45L = _mm_unpacklo_pi8( row4, row5 ); // 40 50 41 51 42 52 43 53
  __m64 tmp45H = _mm_unpackhi_pi8( row4, row5 ); // 44 54 45 55 46 56 47 57

  // 係数転置 ( フェーズ 2 )
  __m64 tmp0123LL = _mm_unpacklo_pi16( tmp01L, tmp23L ); // 00 10 20 30 01 11 21 31
  __m64 tmp0123LH = _mm_unpackhi_pi16( tmp01L, tmp23L ); // 02 12 22 32 03 13 23 33
  __m64 tmp0123HL = _mm_unpacklo_pi16( tmp01H, tmp23H ); // 04 14 24 34 05 15 25 35
  __m64 tmp0123HH = _mm_unpackhi_pi16( tmp01H, tmp23H ); // 06 16 26 36 07 17 27 37
  __m64 tmp4567LL = _mm_unpacklo_pi16( tmp45L, tmp67L ); // 40 50 60 70 41 51 61 71
  __m64 tmp4567LH = _mm_unpackhi_pi16( tmp45L, tmp67L ); // 42 52 62 72 43 53 63 73
  __m64 tmp4567HL = _mm_unpacklo_pi16( tmp45H, tmp67H ); // 44 54 64 74 45 55 65 75
  __m64 tmp4567HH = _mm_unpackhi_pi16( tmp45H, tmp67H ); // 46 56 66 76 47 57 67 77

  // 係数転置 ( フェーズ 3 )
  col0 = _mm_unpacklo_pi32( tmp0123LL, tmp4567LL ); // 00 10 20 30 40 50 60 70
  col1 = _mm_unpackhi_pi32( tmp0123LL, tmp4567LL ); // 01 11 21 31 41 51 61 71
  col2 = _mm_unpacklo_pi32( tmp0123LH, tmp4567LH ); // 02 12 22 32 42 52 62 72
  col3 = _mm_unpackhi_pi32( tmp0123LH, tmp4567LH ); // 03 13 23 33 43 53 63 73
  col4 = _mm_unpacklo_pi32( tmp0123HL, tmp4567HL ); // 04 14 24 34 44 54 64 74
  col5 = _mm_unpackhi_pi32( tmp0123HL, tmp4567HL ); // 05 15 25 35 45 55 65 75
  col6 = _mm_unpacklo_pi32( tmp0123HH, tmp4567HH ); // 06 16 26 36 46 56 66 76
  col7 = _mm_unpackhi_pi32( tmp0123HH, tmp4567HH ); // 07 17 27 37 47 57 67 77
}

unpack で挟み込みながら変換していく。
MMX で 4x4 の時や、SSE2 で 8x8 の時も基本的に同じ考え方でいける。
ただ、シャッフル系の命令が増えてきているので、そちらを使った方が速くなるかもしれない。

投稿者 Takenori : 20:56 | トラックバック

逆数の高精度化

割り算は遅いので、逆数の掛け算で処理したいのだが、逆数を求める _mm_rcp_ps は 11bit 精度しかなくて、普通の用途には使い辛い ( 11bit 精度だと、10 進数で3桁程度しかない ) 。
でも、割り算は遅い。
と言うことで、11bit 精度で得られた逆数の近似値をニュートン-ラフソン法 ( Newton-Raphson ) を用いて高精度化する。
SSE のソースコードは以下。
精度は 22bit 程度。_mm_div_ps は 23bit 精度なので、あまり変わらない。
とりあえず、よくわからなくても以下のコードコピペして使えばうまく行く。

// 22bit 精度で逆数を求める
// aの逆数 = 2 * a - rcpa * a * a を用いる
inline __m128 m128_rcp_22bit_ps( const __m128& a ) {
  __m128 xm0 = a;
  __m128 xm1 = _mm_rcp_ps(xm0);
  xm0 = _mm_mul_ps( _mm_mul_ps( xm0, xm1 ), xm1 );
  xm1 = _mm_add_ps( xm1, xm1 );
  return _mm_sub_ps( xm1, xm0 );
}
inline __m128 m128_rcp_22bit_ss( const __m128& a ) {
  __m128 xm0 = a;
  __m128 xm1 = _mm_rcp_ss(xm0);
  xm0 = _mm_mul_ss( _mm_mul_ss( xm0, xm1 ), xm1 );
  xm1 = _mm_add_ss( xm1, xm1 );
  return _mm_sub_ss( xm1, xm0 );
}
inline float rcp_sse( float a ) {
  float  ret;
  __m128 xm0 = _mm_set_ss(a);
  _mm_store_ss( &ret, m128_rcp_22bit_ss(xm0) );
  return ret;
}

これで割り算の多い日も安心。

投稿者 Takenori : 20:57 | トラックバック

平方根の逆数の高精度化

平方根の逆数を求める命令として、_mm_rsqrt_ss / _mm_rsqrt_ps がある。
でも、これは 11bit 精度しかない ( 11bit 精度だと、10 進数で3桁程度 ) 。
と言うことで、11bit 精度で得られた平方根の逆数の近似値をニュートン-ラフソン法 ( Newton-Raphson ) を用いて高精度化する。
SSE のソースコードは以下。

// 22bit 精度で平方根の逆数を求める。
// 1/√x = 3/2 * x0 - 1/2 * a * x0^3 を用いる
inline __m128 m128_rsqrt_22bit_ps( const __m128& a ) {
  static const __declspec(align(16)) float valm05[4] = {-0.5f,-0.5f,-0.5f,-0.5f};
  static const __declspec(align(16)) float val15[4] = {1.5f,1.5f,1.5f,1.5f};
  static const __m128 Q_m0p5 = *(__m128*)valm05;
  static const __m128 Q_1p5 = *(__m128*)val15;
  __m128 xm0 = a;
  __m128 xm1 = _mm_rsqrt_ps( xm0 );
  __m128 xm2 = _mm_mul_ps( xm0, xm1 );
  xm2 = _mm_mul_ps( xm2, xm1 );
  xm2 = _mm_mul_ps( xm2, xm1 );
  xm2 = _mm_mul_ps( xm2, Q_m0p5 );
  xm0 = _mm_mul_ps( xm1, Q_1p5 );
  return _mm_add_ps( xm0, xm2 );
}
inline __m128 m128_rsqrt_22bit_ss( const __m128& a ) {
  // 必要ないけど、4つセット
  static const __declspec(align(16)) float valm05[4] = {-0.5f,-0.5f,-0.5f,-0.5f};
  static const __declspec(align(16)) float val15[4] = {1.5f,1.5f,1.5f,1.5f};
  static const __m128 Q_m0p5 = *(__m128*)valm05;
  static const __m128 Q_1p5 = *(__m128*)val15;
  __m128 xm0 = a;
  __m128 xm1 = _mm_rsqrt_ss( xm0 );
  __m128 xm2 = _mm_mul_ss( xm0, xm1 );
  xm2 = _mm_mul_ss( xm2, xm1 );
  xm2 = _mm_mul_ss( xm2, xm1 );
  xm2 = _mm_mul_ss( xm2, Q_m0p5 );
  xm0 = _mm_mul_ss( xm1, Q_1p5 );
  return _mm_add_ss( xm0, xm2 );
}
inline float rsqrt_sse( float a ) {
  float  ret;
  __m128 xm0 = _mm_set_ss(a);
  _mm_store_ss( &ret, m128_rsqrt_22bit_ss( xm0 ) );
  return ret;
}

これと逆数の高精度化を用いれば、平方根が得られるが、平方根は素直に _mm_sqrt_ss / _mm_sqrt_ps を使った方が速い様子。
また、標準関数の実装にもよると思うが、以下のようにして SSE の平方根を使う方が速いこともあるようだ ( VC2005 では以下のほうが速かった ) 。
精度も 23bit なので、困ることはないと思う。

inline float sqrt_sse( float a ) {
  float  ret;
  __m128 xm0 = _mm_set_ss(a);
  _mm_store_ss( &ret, _mm_sqrt_ss( xm0 ) );
  return ret;
}

平方根の逆数は距離で割る時など、初めから逆数になっているのでかけるだけで済む。

追記:
1/√a * a = √a なので、上の平方根の逆数に逆数云々ではなく、単純に平方根の逆数に元の値を掛ければ平方根が求まる。
こうすれば上の高精度化した平方根の逆数に元の数を掛けた方が _mm_sqrt_ss / _mm_sqrt_ps よりも倍近く速い。
2 ~ 1千万までの平方根を求める時間を計測すると、VC 2005 の sqrt、_mm_sqrt_ss 使用版、上の高精度化版*元の数 では、だいたい 3 : 2 : 1 ぐらいの速度比になった。
でも、状況によってはコンパイラの最適化との兼ね合いか、 _mm_sqrt_ss / _mm_sqrt_ps を使う方が速いこともあるよう。詳しくは調べてない。

さらに追記:
( 2 ~ 1千万 ) / 100.0 で平方根を求める時間を計測すると、45 : 33 : 35 ぐらいになった。
小数の時は _mm_sqrt_ss / _mm_sqrt_ps の方が速いようだ。
_mm_rsqrt_ss / _mm_rsqrt_ps はテーブルを引くらしいのでその辺りの兼ね合いか?

投稿者 Takenori : 20:59 | トラックバック

2008年03月03日

SSE からの整数値の取り出し

_mm_cvtt_ss2si を使えば、切り捨てを使用して 32 ビット整数を取得出来るが、取り出せるのは最下位のみ。
SSE2 があれば、4つ取り出せるが、SSE では出来ない。
そのためシャッフルを使うが、これは以下のようにぐるぐる回して取り出したほうが良い様子。

i = _mm_cvtt_ss2si( a );
a = _mm_shuffle_ps( a, a, _MM_SHUFFLE( 0, 3, 2, 1 ) );
...
i = _mm_cvtt_ss2si( a );
a = _mm_shuffle_ps( a, a, _MM_SHUFFLE( 0, 3, 2, 1 ) );
...
i = _mm_cvtt_ss2si( a );
a = _mm_shuffle_ps( a, a, _MM_SHUFFLE( 0, 3, 2, 1 ) );
...
i = _mm_cvtt_ss2si( a );

これを以下のようにすると見た目すっきりするけど、少し遅いみたい。

i = _mm_cvtt_ss2si( a );
...
i = _mm_cvtt_ss2si( _mm_shuffle_ps( a, a, _MM_SHUFFLE( 1, 1, 1, 1 ) ) );
...
i = _mm_cvtt_ss2si( _mm_shuffle_ps( a, a, _MM_SHUFFLE( 2, 2, 2, 2 ) ) );
...
i = _mm_cvtt_ss2si( _mm_shuffle_ps( a, a, _MM_SHUFFLE( 3, 3, 3, 3 ) ) );

たぶん、下の方法だと余計なレジスタコピーが発生するからだと思う。
他でも xmm レジスタを何個か使っているようなソースコードの場合、計測して違いが出るぐらいの速度差になった。

補足:
_mm_cvtt_ps2pi を 使えば、__m64 に2個の 32 ビット整数を一度に取り出せる。
ただ、__m64 を使うと、その後浮動小数点を用いた計算をする場合、_mm_empty をコールしなければならない。
SSE を使っていると、近くで浮動小数点を使っていることが多いと思うが、その時は _mm_empty を忘れずに呼ばないと計算結果がおかしくなる。
_mm_cvtt_ps2pi で値を取り出したあと、MMX である程度処理するのであれば、2個同時に取り出した方が速いかもしれない。
1個ずつ取り出すか、2個同時に取り出すかは、その後の処理次第で、どちらが速いかを計測してから決めた方が良い。
MMX の2個同時程度であれば、MMXを使わずに普通に書いたほうが速い場合も多い。

投稿者 Takenori : 14:58 | トラックバック

2008年03月06日

ハイブリッド データ オーダー

SIMD を使わない場合、以下のような構造体の配列 ( AoS ) を使うことが多いと思う。

struct Vector { float x, y, z; };
Vector vec[100];

でも、この形だと垂直演算が基本の SIMD では扱い辛い。
そこで、SIMD で処理する場合、データ構造は以下のような配列の構造体 ( SoA ) を使うと思う。

struct Vector {
__declspec(align(16)) float x[100], y[100], z[100];
};

ただ、配列サイズが大きくなるとメモリページの再オープンが発生して速度が低下するらしい。
この境界は一般的に 4KB だとか。
そこで AoS と SoA の両方を併せ持ったデータ構造にすることでより効率的に処理できるようになる。

struct Vector {
__declspec(align(16)) float x[8], y[8], z[8];
};
Vector vec[15];

ここで、要素数を 8 としているのは、キャッシュ・ラインのサイズである 32バイトにあわせるため ( キャッシュ・ラインのサイズは CPU によって異なるとか。 ) 。
ただ、横方向に1個ずつずらしながら処理する必要がある場合など、この方法が取り辛い場合もある。
その場合は、キャッシュ・ラインのサイズやメモリページのブロックサイズなどを考えて工夫しないといけない。

詳しくは、インターネット・ストリーミングSIMD 拡張命令を考慮したアプリケーション・チューニング ( PDF ) に載っている。
と言うか、これを読んだ方がわかりやすい。

投稿者 Takenori : 03:34 | トラックバック

2014年03月30日

AVX2 の 集約 (GATHER)

AVX2 に 集約 (GATHER) と言う命令が増えている。
これを使うとメモリ上で飛び飛びのデータを読める。

例えば画像データを縦方向に読みたい場合は、以下のようにすると縦1列で読める。

static const __declspec(align(32)) unsigned long M256_U32_OFFSET[8] = { 0, 1, 2, 3, 4, 5, 6, 7 };
__m256i offset = _mm256_load_si256( (__m256i const *)M256_U32_OFFSET );
offset = _mm256_mullo_epi32( offset, _mm256_set1_epi32(width) );
__m256i color8 = _mm256_i32gather_epi32( (const int*)src, offset, sizeof(int) );

_mm256_i32gather_epi32 の第3引数は定数値である必要があるみたいなので、そこでオフセットを調整することは出来ないようだ。
変数が指定出来れば、stride をここに指定出来てもう少し楽なんだけど。

このようにメモリ上で連続していなくてもデータを読めるので、かなり扱いやすくなった。
ただ、遅い。
メモリ上で連続していないから仕方ないと思うが……
やはり、従来通り出来るだけメモリ上で連続するようにして処理する方がいいのは変わらない。
一応 _mm_prefetch を使うことで少しは改善する。

投稿者 Takenori : 14:10 | トラックバック

2014年04月01日

水平加算での合計計算

SSE3 であれば、以下のように2回水平加算を実行すると、全要素の合計値が全要素に入る。

sum = _mm_hadd_ps( sum, sum );
sum = _mm_hadd_ps( sum, sum );


では、AVX では 3回実行すれば合計値が求まるかと言うとそうはならない。
前後の入れ替えとインターリーブが間に必要になる。
128bit を超えて計算されない。

sum = _mm256_hadd_ps( sum, sum );
sum = _mm256_hadd_ps( sum, sum );
__m256 rsum = _mm256_permute2f128_ps(sum, sum, 0 << 4 | 1 );
sum = _mm256_unpacklo_ps( sum, rsum );
sum = _mm256_hadd_ps( sum, sum );


ちなみに SSE で合計計算しようとするとシャッフルを使って計算することになる。

__m128 tmp = sum;
sum = _mm_shuffle_ps( sum, tmp, _MM_SHUFFLE(1,0,3,2) );
sum = _mm_add_ps( sum, tmp );
tmp = sum;
sum = _mm_shuffle_ps( sum, tmp, _MM_SHUFFLE(2,3,0,1) );
sum = _mm_add_ps( sum, tmp );

投稿者 Takenori : 17:37 | トラックバック

 
Total : Today : Yesterday :