ストリップマイニング
超高性能プログラミング技術のメモ（７）
高性能プログラミングの技術を忘れないようにノートを書いています。今回は、ストリップマイニングについて書いておこうと思います。
ストリップマイニング
ストリップマイニングは、高性能プログラミングでよく用いられる改変方法です。ループを固定長の長さ毎に分けて実行するようにします。その方法は、ベクトル型コンピュータのベクトライザー（ベクトル化処理器）が行なっている改変と同じで、下図のように、ループを二重ループに分割します。
x86_64アーキテクチャのスカラー型コンピュータにはベクトル計算器は装備されておらず、そのためコンパイラにベクトライザーは搭載されていません。しかし、x86アーキテクチャの拡張命令セットMMX, SSE, AVXなどのSIMD（Single Instruction Multiple Data; 単一命令複数データ）処理系統は、１つの命令で複数のデータを同時に処理するまさにベクトル計算を行います。ベクトル長（１命令で処理できるデータ数）は、MMXで64bit、SSEで128bit、AVXで256bit、最新のAVX-512で512bitとなっています。
x86_64向けのコンパイラには、ベクトライザーは搭載されていませんが、最適化オプションでSIMD化を促進することができます。ただし、SIMD化可能かどうかの判断はコンパイラの実装に依存しており、コンパイラ毎に都度調べてみるしかありません。手元のコンパイラでは、上図のように、ループ内に複数の処理が含まれるとSIMD化を見送り、単一の処理にするとSIMD化されました。
ストリップマイニングは、このようにループを分割し、SIMD化を促進する方法として使うことができます。（上図の場合は、元のループを２つに分けるだけでも構いませんが。）
SIMDに最適化できない場合
具体的に、下記のコードをコンパイルしてみましょう。倍精度実数の配列の加算を２つ同時に実行する関数です。
#include <stdlib.h>
void dxpy2( size_t m, const double* A, const double* B, double* C
, const double* D, const double* E, double* F )
{
size_t i;
for( i=0; i<m; i++ ){
C[i] = A[i] + B[i];
F[i] = D[i] + E[i];
}
}これを、gcc -O3 -Sでコンパイルすると最適化済みのアセンブラコードが出力されます。加算を行なっているループの最深部は、下記のコードになります。
LBB0_4: ## =>This Inner Loop Header: Depth=1
movsd (%rsi,%rax,8), %xmm0 ## xmm0 = mem[0],zero
addsd (%rdx,%rax,8), %xmm0
movsd %xmm0, (%rcx,%rax,8)
movsd (%r8,%rax,8), %xmm0 ## xmm0 = mem[0],zero
addsd (%r9,%rax,8), %xmm0
movsd %xmm0, (%r10,%rax,8)
movsd 8(%rsi,%rax,8), %xmm0 ## xmm0 = mem[0],zero
addsd 8(%rdx,%rax,8), %xmm0
movsd %xmm0, 8(%rcx,%rax,8)
movsd 8(%r8,%rax,8), %xmm0 ## xmm0 = mem[0],zero
addsd 8(%r9,%rax,8), %xmm0
movsd %xmm0, 8(%r10,%rax,8)
addq $2, %rax
cmpq %rax, %rdi
jne LBB0_4
ここで使用されているSSE命令MOVSD、ADDSDは、１要素ずつ処理する命令です。これでは、ベクトル化（SIMD化）されたとは言えません。実際、xmmレジスタは、倍精度実数を２個保持でき、SSE命令MOVPDA、ADDPDであれば、２要素同時に計算できます。
SIMDに最適化できた場合
では、ストリップマイニングを行い、２つの処理を２つループに分割してみます。
include <stdlib.h>
#define STRIP_SIZE 8
#define MIN(x,y) ((x)<(y)?(x):(y))
void dxpy2( size_t m, const double* A, const double* B, double* C
, const double* D, const double* E, double* F )
{
size_t i,ii;
for( i=0; i<m; i+=STRIP_SIZE ){
for( ii=i; ii<MIN(m,i+STRIP_SIZE); ii++ ){
C[ii] = A[ii] + B[ii];
}
for( ii=i; ii<MIN(m,i+STRIP_SIZE); ii++ ){
F[ii] = D[ii] + E[ii];
}
}
}この例題の場合、ストリップマイニングを行わず、シンプルに２つループに分けてもSIMD化されます。ストリップマイニングは、もっと複雑な処理の場合に使うと良いでしょう。
ここで、STRIP_SIZEは８にしています。
では、上記のコードを、gcc -O3 -Sでコンパイルし、最深ループで計算に該当する箇所を抜き出します。（下記コード）
LBB0_24: ## Parent Loop BB0_2 Depth=1
## => This Inner Loop Header: Depth=2
movupd -48(%rbx,%rdx,8), %xmm0
movupd -32(%rbx,%rdx,8), %xmm1
movupd -48(%rcx,%rdx,8), %xmm2
addpd %xmm0, %xmm2
movupd -32(%rcx,%rdx,8), %xmm0
addpd %xmm1, %xmm0
movupd %xmm2, -48(%r10,%rdx,8)
movupd %xmm0, -32(%r10,%rdx,8)
movupd -16(%rbx,%rdx,8), %xmm0
movupd (%rbx,%rdx,8), %xmm1
movupd -16(%rcx,%rdx,8), %xmm2
addpd %xmm0, %xmm2
movupd (%rcx,%rdx,8), %xmm0
addpd %xmm1, %xmm0
movupd %xmm2, -16(%r10,%rdx,8)
movupd %xmm0, (%r10,%rdx,8)
addq $8, %rdx
addq $2, %r15
jne LBB0_24
使用されてるSSE命令が、MOVUPDとADDPDに変わっています。これは、倍精度実数を２要素同時に処理する命令で、ADDPDが4回呼ばれていることから、８要素が1回のループで計算されることがわかります。これは、STRIP_SIZEで指定した個数に等しいです。
上記のアセンブラコードを見ると、高性能化のポイントがいくつかあることが分かります。１つは、MOVUPD命令を使っていることです。これを全てMOVAPD命令にできれば、高速化されるはずです。そのためには、引数で入力されたポインタのアライメントによって場合分け行わないといけません。２つ目は、xmmレジスタが3つしか使われていません。xmmレジスタは、全部で16個あります。このままでは、xmm3~xmm15が遊んでいることになります。
まとめ
今回は、ストリップマイニングについてメモしました。しかし、ストリップマイニング技術は、ベクトル型コンピュータでもない限り、あまり使わないかもしれません。
次回は、キャッシュの利用効率を高めるブロッキングについて書こうと思います。