もう少し高速になりました。｜行列積高速化＃２９
前回の記事で、OpenBLASよりも7~8%性能が劣っていました。その後、いろいろ試していたのですが、2~4%ほど高速になったので、性能が向上したポイントを記録しておこうと思います。
ピーク性能比率でみると95~97%、FLOPS値でみるとおよそ45GFLOPSくらいなりました。それでも、OpenBLASよりも劣っているので、まだ工夫の余地はあるのでしょう・・・。
プリフェッチの見直し
カーネル関数のプリフェッチに関しては第18回で書いていますが、第25回でカーネル関数を書き換えたときにプリフェッチ命令を挿入していました。
しかし、この挿入では不完全だったらしく、いろいろ試す中でより高速になるポイントが見えてきたので、記録しておこうと思います。ここでのプリフェッチは、全てL１キャッシュメモリへのプリフェッチです。つまり、PREFETCHT0命令を使用しています。
プリフェッチのポイント
（１）カーネル関数の開始時に、行列Aと行列Bの先頭データをプリフェッチ
（２）最深ループ開始前に、行列Cの書き込み先のデータをプリフェッチ
（３）最深ループ開始時に、３ループ先で使用する行列A・Bをプリフェッチ
（４）最深ループ計算中に、行列A・Bを64Bずつずらして最深ループ内で使用するデータ量分をプリフェッチ
（１）先頭データのプリフェッチ
（１）は、カーネル計算の最深ループの最初のロード命令を少しでも早く処理できるようにするために挿入しました。データがレジスタに格納されていなければ何も計算できないため、処理の初めの方にはどうしてもロード命令が多くなります。
しかし、最初に読み込むデータをCPUは予測できないため、キャッシュメモリに存在しない可能性が高いです。キャッシュメモリに存在しないと、メインメモリからデータをロードすることになり、膨大な時間がかかってしまいます。
結果として、メインメモリからの読み込みを待ってから計算が始まるため、効率が悪化してしまいます。そこで、カーネル関数全体の開始時に、先頭データをプリフェッチしておきます。ループ処理などを行っている間に、L1キャッシュメモリにデータを持ってくることが狙いです。
void myblas_dgemm_kernel_detail(
size_t M, size_t N, size_t K,
double alpha, const double *A, const double *B,
double *C, size_t ldc )
{
double *c0 = C;
double *c1 = C + ldc;
size_t ldc2 = ldc * 2 * sizeof(double);
double alpha4[4] = {alpha,alpha,alpha,alpha};
size_t NQ = N/6;
size_t NR = N%6;
__asm__ __volatile__ (
"\n\t"
"prefetcht0 0*8(%[a])\n\t"
"prefetcht0 0*8(%[b])\n\t"
"prefetcht0 4*8(%[a])\n\t"
"prefetcht0 4*8(%[b])\n\t"
"prefetcht0 8*8(%[b])\n\t"
"\n\t"
:[a]"+r"(A),[b]"+r"(B)
:);
if( NQ ){
size_t n4 = NQ; // unrolling Nこのプリフェッチは、第25回のときにも挿入してあったのですが、これが無いとピーク性能比率が1~2%程度低下する印象です。
（２）行列Cのプリフェッチ
行列Cのプリフェッチは、第18回でも挿入しました。このプリフェッチを挿入すると、ほぼ確実にピーク性能比率1~2%程度の速度向上になるので、必須だと考えた方が良いようです。
このプリフェッチの挿入箇所は、カーネルがN-M-Kループの場合、Mループの開始直後かつKループの開始前です。Kループの終了後、Mループの最終段階で行列Cへの足し込みが行われるので、Kループが実行されている間に、行列Cの保存先をキャッシュに読み込んでおくことが狙いです。
while( n4-- ){
if( M >> 2 ){
size_t m4 = ( M >> 2 ); // unrolling M
while( m4-- ){
__asm__ __volatile__ (
"prefetcht0 0*8(%[c0] )\n\t"
"prefetcht0 0*8(%[c1] )\n\t"
"prefetcht0 0*8(%[c0],%[ldc2] )\n\t"
"prefetcht0 0*8(%[c1],%[ldc2] )\n\t"
"prefetcht0 0*8(%[c0],%[ldc2],2)\n\t"
"prefetcht0 0*8(%[c1],%[ldc2],2)\n\t"
:[c0]"+r"(c0),[c1]"+r"(c1)
:[ldc2]"r"(ldc2)
);
if( K >> 3 ){
size_t k8 = ( K >> 3 ); // Unrolling K
while( k8-- ){行列C は、単に書き込むだけでなく、行列積ABの結果を加算するために既存の値を読み込む必要があります。この読み込み速度を高速にするために、L1キャッシュにデータを配置しておく必要がありました。
（３）最深ループ開始時のプリフェッチ
今回の実装では、最深ループのプリフェッチは、ループ３つ分先を読み込むとちょうど良いようです。例えば、ループ変数をkとすると、k=1のときはk=4でロードするデータを、k=2のときはk=5でロードするデータを事前にプリフェッチしておくというイメージです。ただ、遅くなる場合もあったので、必ずしも上手くいくとは限りませんでした。
現状、最速のMUxNUxKUxLU=4x6x4x2のカーネル関数（参考：カーネル関数の性能推定と改良）では、行列Aを4x4x2=32要素、行列Bを6x4x2=48要素使用します。そのため、３ループ先のデータとは、行列Aは32x3=96要素先、行列Bは48x3=144要素先のデータを指すことになります。
コードは、次のような感じになります。
if( K >> 3 ){
size_t k8 = ( K >> 3 ); // Unrolling K
while( k8-- ){
__asm__ __volatile__ (
"\n\t"
"prefetcht0 96*8(%[a])\n\t"
"prefetcht0 144*8(%[b])\n\t"（４）64Bずつプリフェッチ
最深ループ開始時のプリフェッチだけでは、ループ反復で最初にロードするデータだけがプリフェッチされている状態なので、ループ反復内で使用する全てのデータをプリフェッチしておく必要がありました。例えば、4x6x4x2カーネルの場合、行列Aを32要素、行列Bを48要素だけプリフェッチしておきます。
プリフェッチは２キャッシュライン＝32x2B=64BのデータをL1キャッシュに配置させるので、64Bを単位としてプリフェッチ命令を挿入します。例えば、倍精度浮動小数点の場合、１要素が８Bなので、64B/8B=8要素ずつずらしてプリフェッチを挿入します。4x6x4x2カーネルの場合、１反復内で、行列Aは32要素使用するので、32要素/8要素＝4回のプリフェッチが必要でした。同様に、行列Bは6回のプリフェッチが必要でした。
具体的には、下記のようなプリフェッチを挿入しました。
"prefetcht0 96*8(%[a])\n\t"
"prefetcht0 104*8(%[a])\n\t"
"prefetcht0 112*8(%[a])\n\t"
"prefetcht0 120*8(%[a])\n\t"
"\n\t"
"prefetcht0 144*8(%[b])\n\t"
"prefetcht0 152*8(%[b])\n\t"
"prefetcht0 160*8(%[b])\n\t"
"prefetcht0 168*8(%[b])\n\t"
"prefetcht0 176*8(%[b])\n\t"
"prefetcht0 184*8(%[b])\n\t"ブロッキングサイズの見直し
4x6x4x2カーネルの場合、Nループは６個ずつ計算に使用していきます。ところが、L1キャッシュブロッキングのためのN方向のブロッキングサイズTILE_Nは16に設定していました。これだと、少し遅いN=4のループを毎回使用してしまうため、6の倍数になるように変更しました。
当初のブロッキングサイズは、記事「ブロッキングサイズの改良」で設定したものになっていました。
TILE_N=12とTILE_N=6を比較したところ、TILE_N=6の方が最高速度が上だったため、結果的にTILE_N=6を選択しました。ただし、TILE_N=6の場合は、L1キャッシュ容量が余るため、TILE_Kを2倍にしています。具体的には、次のように変更しました。
変更前
#define MYBLAS_TILE_M 16
#define MYBLAS_TILE_N 16
#define MYBLAS_TILE_K 128変更後
#define MYBLAS_TILE_M 16
#define MYBLAS_TILE_N 6
#define MYBLAS_TILE_K 256おそらく、TILE_Kを長くしたことで、C行列への足し込み計算が減少して、高速になったと考えています。
アセンブラ命令順序の見直し
最後に、こちらの記事で行ったように、カーネル関数のアセンブラ命令の順序をもう一度見直しました。
検討の様子は、途中までですが、下図のようになります。
上図は、上が検討前、下が検討後になっています。基本的な考え方は、レジスタの依存関係を見ながら空白を消していくことです。空白が少なくなったのが分かるかと思います。
しかしながら、やってみたものの高速化効果はほとんどありませんでした。
性能測定結果
上記の変更を行なって性能を測定したところ、ベストケースではピーク性能比率97%代が出るようになりました。
前回の測定結果と比較すると、ややOpenBLASの結果に近付いたのが分かるかと思います。
まとめ
プリフェッチの改良とブロッキングサイズの変更、およびアセンブラ命令順序の変更でピーク性能比率が2~4%ほど向上しました。しかしながら、性能が向上したものの、まだOpenBLASの性能には追いつけません。
OpenBLASは速いですね。工夫の余地が、もう一手か二手ある印象です。
今回の改良は、下記のGitHubのソースコードに反映してあります。現在のコードがどんなものかは下記のsrc/myblas_dgemm_kernel_detail.cをご覧ください。
この記事が気に入ったらサポートをしてみませんか？