高速memcpyの実装例
超高性能プログラミング技術のメモ（６）
高性能プログラミングの技術を忘れないようにメモしています。今回は、memcpyを高性能化してみます。
仕様の確認
C言語プログラマならほとんど知っていると思いますが、まずmemcpy関数のインターフェース仕様を確認します。ここでは、manコマンドで確認してみます。
$> man 3 memcpy
表示されたマニュアルから、SYNOPSISを見ると次のように表示されます。
SYNOPSIS
#include <string.h>
void *
memcpy(void *restrict dst, const void *restrict src, size_t n);memcpy関数は、string.hで定義され、引数にコピー先ポインタdst、コピー元ポインタsrc、コピーサイズnを渡し、コピー後のポインタが返却されてきます。
最もシンプルな実装は、次ようなコードになります。
void* memcpy( void* dst, const void* src, size_t n ){
const unsigned char* x = (const unsigned char*) src;
unsigned char* y = (unsigned char*) dst;
while( n-- ){ *(y++) = *(x++); }
return src;
}このコードは、シンプルですが、１バイトずつコピーするため高速ではありません。
高性能化のポイント
アセンブラ命令で考えると、メモリコピーとはメモリロード命令とメモリストア命令だけで構成されます。演算が全くないので、高性能化するポイントは多くはありません。実際に考えられるポイントは、次の３つくらいです。
①区画の境界にアライメントし、高速な命令を使用する
②１バイトずつではなく、数バイトまとめてコピーする
③ループをずらして、メモリロード時間を稼ぐ
以降では、x86_64アーキテクチャを前提とし、拡張命令AVXとYMMレジスタを使用した実装を考えます。
①アライメント
引数に渡される２つポインタdst,srcは、区画の境界にアライメントされているかどうか分かりません。そのため、アライメントを意識した実装では、４つ場合に分けて実装する必要があります。
ケース０ dstとsrcの両方がアライメントされている。
ケース１ dstはアライメントされていないが、srcはアライメントされている。
ケース２ dstはアライメントされているが、srcはアライメントされていない。
ケース３ dstとsrcの両方がアライメントされていない。
これをコードにすると、次のように書くことができます。
#define ALIGN_SIZE 32 // bytes for YMM
#define ALIGN_CHECK 0x1f // 11111
void* fast_memcpy(void*restrict dst, const void* restrict src, size_t n){
const unsigned char* x = (const unsigned char*) src;
unsigned char* y = (unsigned char*) dst;
uint64_t xalign = ((uint64_t)x)&ALIGN_CHECK;
uint64_t yalign = ((uint64_t)y)&ALIGN_CHECK;
unsigned int align = ( ( ( (yalign>0)?1:0) << 1 ) | ((xalign>0)?1:0) );YMMレジスタは256ビット（＝32バイト）のため、32バイトの区画境界にアライメントしている必要があります。そのため、上記のコードでは、渡されたポインタが32B境界上に在るかを確認し、align=0~3のケース番号を作成しています。AVX最適化の詳細は、Intel最適化ガイドラインの第12章をご覧ください。
②まとめてコピー
キャッシュメモリへ一度にデータ転送するラインサイズは64B（古いアーキテクチャだと32B）でした。このキャッシュラインがロードされると、次の64Bが自動的にプリフェッチされます。そのため、1Bずつロードしているとプリフェッチしたデータが無駄になります。結果、YMMレジスタ（256 bits=32 bytes）全体を使って、32Bを一度にコピーするのが効率的です。
また、x86_64アーキテクチャはYMMレジスタを16個持っており、16個のコピーを同時に処理することができます。メモリロードとメモリストアには、AVX命令であるVMOVDQA命令を使用します。この命令に渡されるメモリアドレスは、32B区画境界にアライメントされていなければいけません。アライメントが保証できない場合は、VMOVDQU命令を使います。
YMMレジスタを全部使った場合、32B*16個＝512Bをまとめてコピーすることができます。下記のコードは、512Bをまとめてコピーする部分の抜粋です。プログラミング言語は、C言語とインラインアセンブラ（AT&T形式）で記述しています。
xx = x + 512;
yy = y + 512;
/*************************************************************************
CASE 0 : DIST and SRC are aligned
**************************************************************************/
/* 512-bytes block */
while( i > 512 ){
__asm__ __volatile__(
"\n\t"
"vmovdqa -16*32(%[x]), %%ymm0 \n\t"
"vmovdqa -15*32(%[x]), %%ymm1 \n\t"
"vmovdqa -14*32(%[x]), %%ymm2 \n\t"
"vmovdqa -13*32(%[x]), %%ymm3 \n\t"
"vmovdqa -12*32(%[x]), %%ymm4 \n\t"
"vmovdqa -11*32(%[x]), %%ymm5 \n\t"
"vmovdqa -10*32(%[x]), %%ymm6 \n\t"
"vmovdqa -9*32(%[x]), %%ymm7 \n\t"
"vmovdqa -8*32(%[x]), %%ymm8 \n\t"
"vmovdqa -7*32(%[x]), %%ymm9 \n\t"
"vmovdqa -6*32(%[x]), %%ymm10\n\t"
"vmovdqa -5*32(%[x]), %%ymm11\n\t"
"vmovdqa -4*32(%[x]), %%ymm12\n\t"
"vmovdqa -3*32(%[x]), %%ymm13\n\t"
"vmovdqa -2*32(%[x]), %%ymm14\n\t"
"vmovdqa -1*32(%[x]), %%ymm15\n\t"
"\n\t"
"vmovdqa %%ymm0 , -16*32(%[y])\n\t"
"vmovdqa %%ymm1 , -15*32(%[y])\n\t"
"vmovdqa %%ymm2 , -14*32(%[y])\n\t"
"vmovdqa %%ymm3 , -13*32(%[y])\n\t"
"vmovdqa %%ymm4 , -12*32(%[y])\n\t"
"vmovdqa %%ymm5 , -11*32(%[y])\n\t"
"vmovdqa %%ymm6 , -10*32(%[y])\n\t"
"vmovdqa %%ymm7 , -9*32(%[y])\n\t"
"vmovdqa %%ymm8 , -8*32(%[y])\n\t"
"vmovdqa %%ymm9 , -7*32(%[y])\n\t"
"vmovdqa %%ymm10 , -6*32(%[y])\n\t"
"vmovdqa %%ymm11 , -5*32(%[y])\n\t"
"vmovdqa %%ymm12 , -4*32(%[y])\n\t"
"vmovdqa %%ymm13 , -3*32(%[y])\n\t"
"vmovdqa %%ymm14 , -2*32(%[y])\n\t"
"vmovdqa %%ymm15 , -1*32(%[y])\n\t"
"\n\t"
"subq $-16*32, %[x]\n\t"
"subq $-16*32, %[y]\n\t"
"\n\t"
:[x]"=r"(xx),[y]"=r"(yy)
:"0"(xx),"1"(yy)
);
i-=512;
}
MOV命令系は、メモリやレジスタ間のデータ転送命令群で、左のオペランドから右のオペランドへデータを転送します。左にメモリアドレス（-16*32(%[x])等）で右にレジスタ（%%ymm0等）が書かれていればロード命令になります。逆に、左にレジスタ（%%ymm0等）で右にメモリアドレス（-16*32(%[y])等）が書かれていればストア命令になります。
上記のコードは、前半で512Bデータをロードし、後半で512Bデータをストアしています。アクセスするメモリアドレスは、必ず32の倍数だけシフトしたアドレス（-n*32(%[x])は、x-n*32の意味）なので、先頭アドレスが32B境界にアライメントされていれば、シフトしたアドレスも必ずアライメントできています。最後に、subq命令（減算）で使用していたポインタを512B分だけシフトし、次の反復に備えています。
③ループずらし、時間稼ぎ
最初のYMM0へのロード命令（VMOVDQA）から、YMM0のストア命令（VMOVDQA）までは、無関係な15個のVMOVDQA命令があります。Intel最適化ガイドラインの付録Cによると、スループット時間は0.25cycleです。これは、４命令が同時に実行できることを示しています。しかし、マイクロアーキテクチャのブロック図からは、2個のロード演算器と2個のストア演算器となっています。すなわち、ロード命令は2個が1cycleで実行でき、16個の連続したロード命令は8cycleで実行できることになります。同様に、16個の連続したストア命令も8cycleで実行されるので、上記のコードは16cycleかかると予測できます。
ロード命令２個とストア命令２個が同時実行可能ということは、ロード命令とストア命令は2個ずつ交互に並べた方が良いことを示しています。パイプラインハザードが起こらないと仮定すると、ロード命令16個とストア命令16個を2個ずつ交互に並べた場合、ロード命令２個とストア命令２個で1cycleなので、全スループット時間は8cycleとなり連続的に並べた場合の半分になります。しかし、ロード命令→ストア命令の順はRAWハザードを発生させてしまうため、命令を近づけるとストア演算器が停止してしまいます。そこで、命令間の依存関係を逆順にするために、ループずらしを行います。
/*************************************************************************
CASE 0 : DIST and SRC are aligned
**************************************************************************/
/* 512-bytes block */
if( i > 512 ){
__asm__ __volatile__(
"\n\t"
"vmovdqa -16*32(%[x]), %%ymm0 \n\t"
"vmovdqa -15*32(%[x]), %%ymm1 \n\t"
"vmovdqa -14*32(%[x]), %%ymm2 \n\t"
"vmovdqa -13*32(%[x]), %%ymm3 \n\t"
"vmovdqa -12*32(%[x]), %%ymm4 \n\t"
"vmovdqa -11*32(%[x]), %%ymm5 \n\t"
"vmovdqa -10*32(%[x]), %%ymm6 \n\t"
"vmovdqa -9*32(%[x]), %%ymm7 \n\t"
"vmovdqa -8*32(%[x]), %%ymm8 \n\t"
"vmovdqa -7*32(%[x]), %%ymm9 \n\t"
"vmovdqa -6*32(%[x]), %%ymm10\n\t"
"vmovdqa -5*32(%[x]), %%ymm11\n\t"
"vmovdqa -4*32(%[x]), %%ymm12\n\t"
"vmovdqa -3*32(%[x]), %%ymm13\n\t"
"vmovdqa -2*32(%[x]), %%ymm14\n\t"
"vmovdqa -1*32(%[x]), %%ymm15\n\t"
"\n\t"
:[x]"=r"(xx),[y]"=r"(yy)
:"0"(xx),"1"(yy)
);
while( i > 512*2 ){
//printf("i=%zu\n",i);
__asm__ __volatile__(
"\n\t"
"vmovdqa %%ymm0 , -16*32(%[y])\n\t"
"vmovdqa 0*32(%[x]), %%ymm0 \n\t"
"vmovdqa %%ymm1 , -15*32(%[y])\n\t"
"vmovdqa 1*32(%[x]), %%ymm1 \n\t"
"vmovdqa %%ymm2 , -14*32(%[y])\n\t"
"vmovdqa 2*32(%[x]), %%ymm2 \n\t"
"vmovdqa %%ymm3 , -13*32(%[y])\n\t"
"vmovdqa 3*32(%[x]), %%ymm3 \n\t"
"vmovdqa %%ymm4 , -12*32(%[y])\n\t"
"vmovdqa 4*32(%[x]), %%ymm4 \n\t"
"vmovdqa %%ymm5 , -11*32(%[y])\n\t"
"vmovdqa 5*32(%[x]), %%ymm5 \n\t"
"vmovdqa %%ymm6 , -10*32(%[y])\n\t"
"vmovdqa 6*32(%[x]), %%ymm6 \n\t"
"vmovdqa %%ymm7 , -9*32(%[y])\n\t"
"vmovdqa 7*32(%[x]), %%ymm7 \n\t"
"vmovdqa %%ymm8 , -8*32(%[y])\n\t"
"vmovdqa 8*32(%[x]), %%ymm8 \n\t"
"vmovdqa %%ymm9 , -7*32(%[y])\n\t"
"vmovdqa 9*32(%[x]), %%ymm9 \n\t"
"vmovdqa %%ymm10 , -6*32(%[y])\n\t"
"vmovdqa 8*32(%[x]), %%ymm8 \n\t"
"vmovdqa %%ymm9 , -7*32(%[y])\n\t"
"vmovdqa 9*32(%[x]), %%ymm9 \n\t"
"vmovdqa %%ymm10 , -6*32(%[y])\n\t"
"vmovdqa 10*32(%[x]), %%ymm10\n\t"
"vmovdqa %%ymm11 , -5*32(%[y])\n\t"
"vmovdqa 11*32(%[x]), %%ymm11\n\t"
"vmovdqa %%ymm12 , -4*32(%[y])\n\t"
"vmovdqa 12*32(%[x]), %%ymm12\n\t"
"vmovdqa %%ymm13 , -3*32(%[y])\n\t"
"vmovdqa 13*32(%[x]), %%ymm13\n\t"
"vmovdqa %%ymm14 , -2*32(%[y])\n\t"
"vmovdqa 14*32(%[x]), %%ymm14\n\t"
"vmovdqa %%ymm15 , -1*32(%[y])\n\t"
"vmovdqa 15*32(%[x]), %%ymm15\n\t"
"\n\t"
"subq $-16*32, %[x]\n\t"
"subq $-16*32, %[y]\n\t"
"\n\t"
:[x]"=r"(xx),[y]"=r"(yy)
:"0"(xx),"1"(yy)
);
i-=512;
}
__asm__ __volatile__(
"\n\t"
"\n\t"
"vmovdqa %%ymm0 , -16*32(%[y])\n\t"
"vmovdqa %%ymm1 , -15*32(%[y])\n\t"
"vmovdqa %%ymm2 , -14*32(%[y])\n\t"
"vmovdqa %%ymm3 , -13*32(%[y])\n\t"
"vmovdqa %%ymm4 , -12*32(%[y])\n\t"
"vmovdqa %%ymm5 , -11*32(%[y])\n\t"
"vmovdqa %%ymm6 , -10*32(%[y])\n\t"
"vmovdqa %%ymm7 , -9*32(%[y])\n\t"
"vmovdqa %%ymm8 , -8*32(%[y])\n\t"
"vmovdqa %%ymm9 , -7*32(%[y])\n\t"
"vmovdqa %%ymm10 , -6*32(%[y])\n\t"
"vmovdqa %%ymm11 , -5*32(%[y])\n\t"
"vmovdqa %%ymm12 , -4*32(%[y])\n\t"
"vmovdqa %%ymm13 , -3*32(%[y])\n\t"
"vmovdqa %%ymm14 , -2*32(%[y])\n\t"
"vmovdqa %%ymm15 , -1*32(%[y])\n\t"
"\n\t"
"subq $-16*32, %[x]\n\t"
"subq $-16*32, %[y]\n\t"
"\n\t"
:[x]"=r"(xx),[y]"=r"(yy)
:"0"(xx),"1"(yy)
);
i-=512;
}上記のコードは、前半のロード命令16個をループの前に出し、後半のストア命令16個をループの後ろに出しました。ループ内では、ロード命令16個をストア命令の背後に移動し、512B先のロード命令に書き換えます（例：-16*32(%[x])を0*32(%[x])に書き換える）。元のループでは、これは次の反復のロード命令になります。
ストア命令→ロード命令はWAR型のため、パイプラインハザードが発生しません。そのため、同じYMM0レジスタを使っていても、ストア命令の直後までロード命令を移動することができます。移動することで、メモリロードの時間稼ぎも最大化でき、一石二鳥になります。（上記コードの場合は、subq命令分しか追加されません。）
効果測定
では、作成した関数をfast_memcpyと名付けて、組み込みのmemcpyと速度比較をしてみます。コピー時間は、データサイズに依存するので、64Bから2倍していき２GBまで測定してみます。キャッシュを利用しないように、memcpyとfast_memcpyは同じサイズの別のデータを使用しました。測定の結果が、下記の表になります。
memcpy time列とfast_memcpy time列が測定された時間です。memcpy performance列とfast_memcpy performance列は、Datasizeを測定時間で割った値で、データ転送速度（スループット）を表します。speed-up ratioは、memcpyの測定時間をfast_memcpyの測定時間で割った値で、fast_memcpyが何倍高速化されたかを表します。
speed-up ratioを見ると、16KB〜1MBは10倍以上、4MB〜64MBまでは2〜5倍、128MB〜512MBは1〜2倍と少々落ち着き、1GB〜2GBでは再び2〜3倍に高速化されていることが分かります。測定時間とデータ転送速度の違いは、数値では分かりにくいのでグラフにしてみます。
上の左図はデータサイズに対する測定時間の比較（縦軸対数）、右図はデータサイズに対するデータ転送速度の比較です。左図からは、fast_memcpyがどのデータサイズに対しても高速であることが分かります。一方、右図からは、fast_memcpyのデータ転送速度が段階的に落ちてきていることが分かります。これは、データがキャッシュに残っていた可能性が考えられます。
データサイズがキャッシュ容量を超えると、下層のキャッシュを使い始めるため、このような段階構造が見えることがあります。最下層に到達した128MB以上は、アルゴリズムの純粋な比較と考えられます。これによると、だいたい1.2~3.2倍高速化できたと考えて良さそうです。
まとめ
今回は、標準関数であるmemcpyを高速化してみました。
サンプルコードは、GitHubアップロードしてあります。（2020/5/23 ダウロード先を変更）