行列スケール関数の最適化｜行列積高速化#23

この記事は、以下の記事を分割したものです。
［元の記事］行列積計算を高速化してみる
一括で読みたい場合は、元の記事をご覧ください。

スケール関数の最適化も、コピー関数の最適化と同じの手続きを進めます。ただ、スケール関数はまだアンローリングを行っていないため、アンローリングから始めます。また、計算時間は非常に短いため、今回はアライメントの場合分けを見送ります。

最適化の手続き
（１）アンローリング
（２）ベクトライズ、アセンブラ化
（３）プリフェッチの挿入
（４）ループシフト／命令順序の入れ換え

23-1. アンローリング

行列Cにbetaをかけるスケール計算は、計算量がO(N^2)なので行列積計算ではほとんど負荷になりませんが、可能な限り高速化するため、まずはアンローリングします。

アンローリングは、まず段数を決める必要があります。アンローリングで処理量を増やしレジスタを使いまわせなくなると効率が低下するため、アンロール段数はレジスタ数によって制約されます。

x86_64アーキテクチャでは、YMMレジスタが16個装備されています。betaの値を保持するのに1個使いたいので、行列Cのデータには残りの15個を使うことができます。

また、YMMレジスタには、倍精度の値を４つ格納できるので、行列Cの連続メモリ方向のアンロール段数は、4の倍数にすると好都合です。行列Cは、M×N行列でM方向に連続ですので、M方向のアンロール段数は4の倍数に限定します。

さらに、アンロール段数が2の倍数になっていると、条件分岐が軽量なビット演算でできるようになります。そこで、M方向もN方向もアンロール段数は2の倍数に設定したいところです。

以上の条件から、アンロール段数の設定候補として、下記の２つを検討しました。

（１）8x4段アンロール　… レジスタ数の限界値。レジスタはかなり余る。
（２）16x4段アンロール ... レジスタが１つ足りない。使い回しでカバーできる。

例えば、8x4段アンロールの場合のプログラムコードは次のような形になります。

＜アンローリング前＞

size_t n = N;
while( n-- ){
    size_t m = M;
    while( m-- ){
        *C=beta*(*C);
         C++;
    }
    C = C - M + ldc;
}
C = C - ldc*N; // retern to head of pointer.


＜アンローリング後＞

double ymm0 ,ymm1 ,ymm2 ,ymm3 ;
double ymm4 ,ymm5 ,ymm6 ,ymm7 ;
double ymm8 ,ymm9 ,ymm10,ymm11;
double ymm12,ymm13,ymm14,ymm15;

ymm15 = beta;

size_t n = N;
if( n >> 2 ){
   size_t n4 = ( n >> 2 ); // unrolling with 4 elements
   while( n4-- )
   {
       size_t m = M;
       if( m >> 3 ){
           size_t m8 = ( m >> 3 );
           while( m8-- ){
               size_t m8a = 4; // YMM Vector-length
               size_t m8b = 4; // YMM Vector-length
               while( m8a-- ){ // vectorizing
                   ymm0 = *(C+0*ldc);
                   ymm1 = *(C+1*ldc);
                   ymm2 = *(C+2*ldc);
                   ymm3 = *(C+3*ldc);
                   ymm0 *= ymm15;
                   ymm1 *= ymm15;
                   ymm2 *= ymm15;
                   ymm3 *= ymm15;
                   *(C+0*ldc) = ymm0;
                   *(C+1*ldc) = ymm1;
                   *(C+2*ldc) = ymm2;
                   *(C+3*ldc) = ymm3;
                   C++;
               }
               while( m8b-- ){ // vectorizing
                   ymm8 = *(C+0*ldc);
                   ymm9 = *(C+1*ldc);
                   ymm10= *(C+2*ldc);
                   ymm11= *(C+3*ldc);
                   ymm8  *= ymm15;
                   ymm9  *= ymm15;
                   ymm10 *= ymm15;
                   ymm11 *= ymm15;
                   *(C+0*ldc) = ymm8 ;
                   *(C+1*ldc) = ymm9 ;
                   *(C+2*ldc) = ymm10;
                   *(C+3*ldc) = ymm11;
                   C++;
               }

           }
       }
       if( m & 0x4 ){
           size_t m4 = 4;
           while( m4-- )
           {
               ymm0 = *(C+0*ldc);
               ymm1 = *(C+1*ldc);
               ymm2 = *(C+2*ldc);
               ymm3 = *(C+3*ldc);
               ymm0 *= ymm15;
               ymm1 *= ymm15;
               ymm2 *= ymm15;
               ymm3 *= ymm15;
               *(C+0*ldc) = ymm0;
               *(C+1*ldc) = ymm1;
               *(C+2*ldc) = ymm2;
               *(C+3*ldc) = ymm3;
               C++;
           }
       }
       if( m & 0x3 ){
           size_t mr = ( m & 0x3 );
           while( mr-- ){
               ymm0 = *(C+0*ldc);
               ymm1 = *(C+1*ldc);
               ymm2 = *(C+2*ldc);
               ymm3 = *(C+3*ldc);
               ymm0 *= ymm15;
               ymm1 *= ymm15;
               ymm2 *= ymm15;
               ymm3 *= ymm15;
               *(C+0*ldc) = ymm0;
               *(C+1*ldc) = ymm1;
               *(C+2*ldc) = ymm2;
               *(C+3*ldc) = ymm3;
               C++;
           }
       }
       C = C - M + 4*ldc;
   }
}
...以下省略...

M≧８のループ内では、YMMレジスタのベクトル長に合わせて、４要素ずつのループ２つ（m8aとm8bのループ）に分割しています。このコードでは、レジスタがymm4~ymm7、ymm12~ymm14と7個も使用されていません。これを使い切るために16x4段アンロールにすると、今度はレジスタが１つ不足します。

以降は、有料記事にさせていただきます。

次の記事