ループ交換1回目｜行列積高速化#7

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

L1キャッシュブロッキングループの順序を、下記のようにKループを外側に変更しました。こうすると、行列B2の最初のメモリロード回数を減らすことができるためです。

＜変更前＞
     // L1 cache
     for( size_t j1=j2; j1<j2+N2; j1+=MIN(N-j1,MYBLAS_TILE_N ) ){
     for( size_t i1=i2; i1<i2+M2; i1+=MIN(M-i1,MYBLAS_TILE_M ) ){
     for( size_t k1=k2; k1<k2+K2; k1+=MIN(K-k1,MYBLAS_TILE_K ) ){

         C = C + ldc*j1 + i1;
         A2 = A2 + lda2*(k1-k2) + (i1-i2);
         B2 = B2 + ldb2*(j1-j2) + (k1-k2);

＜変更後＞
     // L1 cache
     for( size_t k1=k2; k1<k2+K2; k1+=MIN(K-k1,MYBLAS_TILE_K ) ){
     for( size_t j1=j2; j1<j2+N2; j1+=MIN(N-j1,MYBLAS_TILE_N ) ){
         B2 = B2 + ldb2*(j1-j2) + (k1-k2);
     for( size_t i1=i2; i1<i2+M2; i1+=MIN(M-i1,MYBLAS_TILE_M ) ){

         C = C + ldc*j1 + i1;
         A2 = A2 + lda2*(k1-k2) + (i1-i2);

このように、ブロッキングループは、行列Aと行列Bが共通して依存するKループを外側に持っていくと、メモリアクセスが減少します。

これに合わせて、前節で作成したコピー処理でも、同様にループ順序を入れ替え、Kループを外側へ移動する必要があります。常に、行列積計算ループとコピーループの順序は整合させておかないと、後々困ったことになります。

ループ交換をしたコピー処理は次のようになります。

// On L2-Cache Copy for A
A = A + lda*k2 + i2;
for( size_t k1=k2; k1<k2+K2; k1+=MIN(K-k1,MYBLAS_TILE_K ) ){
 size_t K1 = MIN(MYBLAS_TILE_K ,K-k1);
 for( size_t i1=i2; i1<i2+M2; i1+=MIN(M-i1,MYBLAS_TILE_M ) ){
   size_t M1 = MIN(MYBLAS_TILE_M ,M-i1);
   for( size_t i =i1; i < i1+M1; i++ ){
     for( size_t k =k1; k < k1+K1; k++ ){
       (*A2) = (*A);
       A += lda ;
       A2+= lda2;
     }
     A  = A  - lda *K1 + 1;  
     A2 = A2 - lda2*K1 + 1;  
   }
 }
 A = A  - M2 + lda *K1;  
 A2= A2 - M2 + lda2*K1;  
}
A = A  - lda *K2;
A2= A2 - lda2*K2;
A = A - lda*k2 - i2; // return to head 

// On L2-Cache Copy for B
B = B + ldb*j2 + k2;
for( size_t k1=k2; k1<k2+K2; k1+=MIN(K-k1,MYBLAS_TILE_K ) ){
 size_t K1 = MIN(MYBLAS_TILE_K ,K-k1);
 for( size_t j1=j2; j1<j2+N2; j1+=MIN(N-j1,MYBLAS_TILE_N ) ){
   size_t N1 = MIN(MYBLAS_TILE_N ,N-j1);
   for( size_t j =j1; j < j1+N1; j++ ){
     for( size_t k =k1; k < k1+K1; k++ ){
       *B2=*B;
       B++;
       B2++;
     }
     B  = B  - K1 + ldb ;    
     B2 = B2 - K1 + ldb2;    
   }
 }
 B  = B  - ldb *N2 + K1; 
 B2 = B2 - ldb2*N2 + K1; 
}
B  = B  - K2;
B2 = B2 - K2;
B  = B  - ldb*j2 - k2; // return to head 

では、この計算速度を測定してみましょう。今回はメモリロードがわずかに減少するだけなので、あまり速度は変化しないことが予想されます。

<ループ交換前＞

Max  Peak MFlops per Core: 52800 MFlops 
Base Peak MFlops per Core: 46400 MFlops 
size  , elapsed time[s],          MFlops,   base ratio[%],    max ratio[%] 
   16,     5.10216E-05,         175.612,        0.378474,        0.332598 
   32,     0.000152111,         451.039,        0.972067,         0.85424 
   64,     0.000288963,          1856.9,         4.00194,         3.51686 
  128,      0.00190306,         2229.81,         4.80563,         4.22313 
  256,        0.015748,         2143.19,         4.61895,         4.05908 
  512,        0.137929,         1951.88,         4.20665,         3.69675 
 1024,         1.01311,          2122.8,         4.57501,         4.02046 
 2048,         7.90767,         2174.15,         4.68566,          4.1177 

＜ループ交換後＞

Max  Peak MFlops per Core: 52800 MFlops 
Base Peak MFlops per Core: 46400 MFlops 
size  , elapsed time[s],          MFlops,   base ratio[%],    max ratio[%] 
   16,     3.40939E-05,         262.804,        0.566388,        0.497735 
   32,     0.000128031,         535.871,         1.15489,         1.01491 
   64,     0.000282049,         1902.42,         4.10004,         3.60307 
  128,      0.00202417,         2096.39,         4.51808,         3.97043 
  256,       0.0162041,         2082.87,         4.48894,         3.94483 
  512,        0.131743,         2043.54,         4.40418,         3.87034 
 1024,        0.977956,         2199.11,         4.73945,         4.16497 
 2048,         7.88368,         2180.76,         4.69992,         4.13024 

測定の結果、予想通り、速度はあまり変化はありませんでした。しかし、キャッシュチューニングでは、削減できるメモリアクセスは徹底的に減らしておく必要があります。

次の記事

元の記事はこちらです。

ソースコードはGitHubで公開しています。

jyam45/dgemm