interleave を最適化するお話。SIMD 型を束ねるクラスを作ったときの事だけど、一般化して残しておく。
まず次のような極小の静的配列がある。boost::array とかでもいい。
template<typename T, int N>
struct array
{
T v[N];
};
N は 1 から 4 くらいを想定している。このインスタンスがいくつかある時、これらをインターリーブしたい。使用例としては複素数とか XYZ 座標とか RGBA とか。
たとえば 4 つの変数をインターリーブする場合は次のような関数で実現出来る。
//cl /c /O1 /FA interleave.cpp
#include <string.h>
template<typename T, int N>
struct array
{
T v[N];
};
template<typename T, int N>
inline void interleave(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d)
{
enum { nTotal = N * 4 };
T tmp[nTotal];
memcpy(&tmp[N*0], a.v, sizeof(a.v));
memcpy(&tmp[N*1], b.v, sizeof(b.v));
memcpy(&tmp[N*2], c.v, sizeof(c.v));
memcpy(&tmp[N*3], d.v, sizeof(d.v));
for ( int i = 0; i < N; ++i )
{
int iA = (i + N*0) * N;
int iB = (i + N*1) * N;
int iC = (i + N*2) * N;
int iD = (i + N*3) * N;
a.v[i] = tmp[iA % nTotal + iA / nTotal];
b.v[i] = tmp[iB % nTotal + iB / nTotal];
c.v[i] = tmp[iC % nTotal + iC / nTotal];
d.v[i] = tmp[iD % nTotal + iD / nTotal];
}
}
template
inline void interleave<int, 2>(
array<int, 2> & a,
array<int, 2> & b,
array<int, 2> & c,
array<int, 2> & d);
memcpy() は説明用に簡略化しているだけなのでここでは深く考えない。
上のコードでは試しに要素数が 2 個の int 配列として実体化している。これを VC++ 2010 SP1 (32-bit) でコンパイルすると次のようになる。
; for より前を省略
$LL3@interleave:
; Line 31
lea eax, DWORD PTR [esi-8] ; 添字の計算
cdq
mov ebx, edi
idiv ebx
add edx, eax
mov eax, DWORD PTR _tmp$[ebp+edx*4] ; 配列のコピー
mov edx, DWORD PTR tv425[ebp]
mov DWORD PTR [edx+ecx], eax
; Line 32
lea eax, DWORD PTR [esi-4] ; 添字の計算
cdq
idiv ebx
add edx, eax
mov eax, DWORD PTR _tmp$[ebp+edx*4] ; 配列のコピー
mov DWORD PTR [ecx], eax
; Line 33
mov eax, esi ; 添字の計算
cdq
idiv ebx
add edx, eax
mov eax, DWORD PTR _tmp$[ebp+edx*4] ; 配列のコピー
mov edx, DWORD PTR tv429[ebp]
mov DWORD PTR [edx+ecx], eax
; Line 34
lea eax, DWORD PTR [esi+4] ; 添字の計算
cdq
idiv ebx
add esi, 2
add edx, eax
mov eax, DWORD PTR _tmp$[ebp+edx*4] ; 配列のコピー
mov edx, DWORD PTR tv433[ebp]
mov DWORD PTR [edx+ecx], eax
add ecx, 4
cmp esi, 12 ; 0000000cH
jl SHORT $LL3@interleave
; for より後を省略
添字の計算に除算が用いられている。プロセッサによるが、除算には 10 から 50 クロックほどかかる。今回 N は 2 であり配列はキャッシュメモリに乗っていると仮定すると、処理の殆どは除算ということになる。
しかし N が 2 であるときの nTotal は 8 と決まっているわけで、iA % 8 は iA & 7 に出来るはず。ビット演算は大抵 1 クロック以下である。なのにそうしないのは VC++ の最適化が弱く、iA が正の値しか取らないという事実を知らないから。iA が負の場合は iA & 7 では余りが求められないので素直に除算にしたのだろう。
そこで iA や他の添字を符号なしにするとこうなる。
//cl /c /O1 /FA interleave2.cpp
#include <string.h>
template<typename T, int N>
struct array
{
T v[N];
};
template<typename T, int N>
inline void interleave(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d)
{
enum { nTotal = N * 4 };
T tmp[nTotal];
memcpy(&tmp[N*0], a.v, sizeof(a.v));
memcpy(&tmp[N*1], b.v, sizeof(b.v));
memcpy(&tmp[N*2], c.v, sizeof(c.v));
memcpy(&tmp[N*3], d.v, sizeof(d.v));
for ( unsigned i = 0; i < N; ++i )
{
unsigned iA = (i + N*0) * N;
unsigned iB = (i + N*1) * N;
unsigned iC = (i + N*2) * N;
unsigned iD = (i + N*3) * N;
a.v[i] = tmp[iA % nTotal + iA / nTotal];
b.v[i] = tmp[iB % nTotal + iB / nTotal];
c.v[i] = tmp[iC % nTotal + iC / nTotal];
d.v[i] = tmp[iD % nTotal + iD / nTotal];
}
}
template
inline void interleave<int, 2>(
array<int, 2> & a,
array<int, 2> & b,
array<int, 2> & c,
array<int, 2> & d);
; for より前を省略
$LL3@interleave:
; Line 31
lea edx, DWORD PTR [eax-8]
shr edx, 3
mov esi, eax
and esi, 7
add edx, esi
mov edx, DWORD PTR _tmp$[ebp+edx*4]
mov DWORD PTR [edi+ecx], edx
; Line 32
lea edx, DWORD PTR [eax-4]
shr edx, 3
lea esi, DWORD PTR [eax-4]
and esi, 7
add edx, esi
mov edx, DWORD PTR _tmp$[ebp+edx*4]
mov DWORD PTR [ecx], edx
; Line 33
mov edx, eax
shr edx, 3
mov esi, eax
and esi, 7
add edx, esi
mov edx, DWORD PTR _tmp$[ebp+edx*4]
mov DWORD PTR [ebx+ecx], edx
; Line 34
lea edx, DWORD PTR [eax+4]
mov esi, edx
shr esi, 3
and edx, 7
add esi, edx
mov edx, DWORD PTR _tmp$[ebp+esi*4]
mov esi, DWORD PTR tv467[ebp]
mov DWORD PTR [esi+ecx], edx
add eax, 2
add ecx, 4
cmp eax, 12 ; 0000000cH
jb SHORT $LL3@interleave
; for より後を省略
除算が消えた。しかしまだループが残っている。ループをアンロールすれば、配列の添字は計算しなくても決まるはずだ。しかし N が template パラメータなので、ハードコーディングでアンロールする事は出来ない。そこで、template の特殊化を使ってアンロールする。
//cl /c /O1 /FA interleave3.cpp
#include <string.h>
#if defined(_MSC_VER)
# define FORCEINLINE __forceinline
#elif defined(__GNUC__)
# define FORCEINLINE __attribute__((always_inline))
#endif
template<typename T, int N>
struct array
{
T v[N];
};
template<typename T, int N, int i>
FORCEINLINE void interleave_base(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d,
T tmp[])
{
enum { nTotal = N * 4 };
unsigned iA = (i + N*0) * N;
unsigned iB = (i + N*1) * N;
unsigned iC = (i + N*2) * N;
unsigned iD = (i + N*3) * N;
a.v[i] = tmp[iA % nTotal + iA / nTotal];
b.v[i] = tmp[iB % nTotal + iB / nTotal];
c.v[i] = tmp[iC % nTotal + iC / nTotal];
d.v[i] = tmp[iD % nTotal + iD / nTotal];
}
template<typename T, int N, int i=N-1>
struct interleave_helper
{
static FORCEINLINE void interleave(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d,
T tmp[])
{
interleave_base<T, N, i>(a, b, c, d, tmp);
interleave_helper<T, N, i-1>::interleave(a, b, c, d, tmp);
}
};
template<typename T, int N>
struct interleave_helper<T, N, 0>
{
static FORCEINLINE void interleave(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d,
T tmp[])
{
interleave_base<T, N, 0>(a, b, c, d, tmp);
}
};
template<typename T, int N>
FORCEINLINE void interleave(
array<T, N> & a,
array<T, N> & b,
array<T, N> & c,
array<T, N> & d)
{
enum { nTotal = N * 4 };
T tmp[nTotal];
memcpy(&tmp[N*0], a.v, sizeof(a.v));
memcpy(&tmp[N*1], b.v, sizeof(b.v));
memcpy(&tmp[N*2], c.v, sizeof(c.v));
memcpy(&tmp[N*3], d.v, sizeof(d.v));
interleave_helper<T, N>::interleave(a, b, c, d, tmp);
}
template
FORCEINLINE void interleave<int, 2>(
array<int, 2> & a,
array<int, 2> & b,
array<int, 2> & c,
array<int, 2> & d);
; for より前を省略
mov eax, DWORD PTR _tmp$[ebp+8]
mov ecx, DWORD PTR _tmp$[ebp+28]
mov DWORD PTR [esi+4], eax
mov eax, DWORD PTR _tmp$[ebp+24]
mov DWORD PTR [edi+4], eax
mov eax, DWORD PTR _tmp$[ebp+12]
mov DWORD PTR [ebx+4], eax
mov eax, DWORD PTR _d$GSCopy$[ebp]
mov DWORD PTR [eax+4], ecx
mov ecx, DWORD PTR _tmp$[ebp]
mov DWORD PTR [esi], ecx
mov ecx, DWORD PTR _tmp$[ebp+16]
mov DWORD PTR [edi], ecx
mov ecx, DWORD PTR _tmp$[ebp+4]
mov DWORD PTR [ebx], ecx
mov ecx, DWORD PTR _tmp$[ebp+20]
add esp, 48 ; 00000030H
mov DWORD PTR [eax], ecx
; for より後を省略
全く演算のない、ただのメモリのコピーになった。
実は、g++ 4.5 では -O1 でも一番最初のコードでメモリのコピーにまで最適化される。clang++ 2.9 だと一番最初のコードで 2 番目のコードのような & とビットシフトにまで最適化される。
