ループ交換2回目｜行列積高速化#10

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

さて、カーネル関数を切り出したため、myblas_dgemm_main関数はだいぶスッキリしました。

// scaling beta*C
block2d_info_t infoC = {M,N,1,1};
myblas_dgemm_scale2d(beta,C,ldc,&infoC);

double*   A2 = calloc( MYBLAS_BLOCK_M*MYBLAS_BLOCK_K, sizeof(double) );
size_t  lda2 = MYBLAS_BLOCK_M; 
double*   B2 = calloc( MYBLAS_BLOCK_K*MYBLAS_BLOCK_N, sizeof(double) );
size_t  ldb2 = MYBLAS_BLOCK_K; 

// L3 cache
for( size_t j3=0 ; j3<N; j3+=MIN(N-j3,MYBLAS_PANEL_N) ){
for( size_t i3=0 ; i3<M; i3+=MIN(M-i3,MYBLAS_PANEL_M) ){
for( size_t k3=0 ; k3<K; k3+=MIN(K-k3,MYBLAS_PANEL_K) ){
   size_t M3 = MIN(MYBLAS_PANEL_M ,M-i3);
   size_t N3 = MIN(MYBLAS_PANEL_N ,N-j3);
   size_t K3 = MIN(MYBLAS_PANEL_K ,K-k3);
   // L2 cache
   for( size_t j2=j3; j2<j3+N3; j2+=MIN(N-j2,MYBLAS_BLOCK_N) ){
   for( size_t i2=i3; i2<i3+M3; i2+=MIN(M-i2,MYBLAS_BLOCK_M) ){
   for( size_t k2=k3; k2<k3+K3; k2+=MIN(K-k2,MYBLAS_BLOCK_K) ){
       size_t M2 = MIN(MYBLAS_BLOCK_M ,M-i2);
       size_t N2 = MIN(MYBLAS_BLOCK_N ,N-j2);
       size_t K2 = MIN(MYBLAS_BLOCK_K ,K-k2);

       // On L2-Cache Copy for A
       block2d_info_t infoA = {M2,K2,MYBLAS_TILE_M,MYBLAS_TILE_K};
       myblas_dgemm_copy_t(A+lda*k2+i2,lda,A2,&infoA);

       // On L2-Cache Copy for B
       block2d_info_t infoB = {K2,N2,MYBLAS_TILE_K,MYBLAS_TILE_N};
       myblas_dgemm_copy_n(B+ldb*j2+k2,ldb,B2,&infoB);

       // L1 cache     
       block3d_info_t info3d = {M2,N2,K2,MYBLAS_TILE_M,MYBLAS_TILE_N,MYBLAS_TILE_K};
       myblas_dgemm_kernel(alpha,A2,B2,C+ldc*j2+i2,ldc,&info3d);

   }}}

}}}

free(A2);
free(B2);

ここでは、L2キャッシュをブロッキングしている３重ループの順番について考えます。

現在、ループ変数k2のループ（k2ループ）が最内で反復していますが、k2は行列Aと行列Bの両方に関係しているため、毎回コピーを行わなければならなくなっています。

もし、k2ループを最も外側に持っていき、k2->i2->j2の順序にするとどうでしょうか？

そうすると、行列Aはj2と無関係なため、行列Aのコピー関数呼び出しは、最内ループの外側に移すことができます。したがって、コピーの回数が減り、メモリアクセス回数も減るため、高速化されることが期待されます。

行列Aと行列Bのコピー関数のどちらを外側に移すかという問題は、行列Aのコピー関数の呼び出し回数を減らすほうが正解です。これは、行列Aのコピー関数は、ストライドアクセスのため、やや遅い処理になっているためです。

順序を入れ替えたプログラムコードは次の通りです。

// scaling beta*C
block2d_info_t infoC = {M,N,1,1};
myblas_dgemm_scale2d(beta,C,ldc,&infoC);

double*   A2 = calloc( MYBLAS_BLOCK_M*MYBLAS_BLOCK_K, sizeof(double) );
double*   B2 = calloc( MYBLAS_BLOCK_K*MYBLAS_BLOCK_N, sizeof(double) );

// L3 cache
for( size_t j3=0 ; j3<N; j3+=MIN(N-j3,MYBLAS_PANEL_N) ){
for( size_t i3=0 ; i3<M; i3+=MIN(M-i3,MYBLAS_PANEL_M) ){
for( size_t k3=0 ; k3<K; k3+=MIN(K-k3,MYBLAS_PANEL_K) ){
   size_t M3 = MIN(MYBLAS_PANEL_M ,M-i3);
   size_t N3 = MIN(MYBLAS_PANEL_N ,N-j3);
   size_t K3 = MIN(MYBLAS_PANEL_K ,K-k3);

   // L2 cache
   for( size_t k2=k3; k2<k3+K3; k2+=MIN(K-k2,MYBLAS_BLOCK_K) ){
     size_t K2 = MIN(MYBLAS_BLOCK_K ,K-k2);

     for( size_t i2=i3; i2<i3+M3; i2+=MIN(M-i2,MYBLAS_BLOCK_M) ){
       size_t M2 = MIN(MYBLAS_BLOCK_M ,M-i2);

       // On L2-Cache Copy for A
       block2d_info_t infoA = {M2,K2,MYBLAS_TILE_M,MYBLAS_TILE_K};
       myblas_dgemm_copy_t(A+lda*k2+i2,lda,A2,&infoA);

       for( size_t j2=j3; j2<j3+N3; j2+=MIN(N-j2,MYBLAS_BLOCK_N) ){
         size_t N2 = MIN(MYBLAS_BLOCK_N ,N-j2);

         // On L2-Cache Copy for B
         block2d_info_t infoB = {K2,N2,MYBLAS_TILE_K,MYBLAS_TILE_N};
         myblas_dgemm_copy_n(B+ldb*j2+k2,ldb,B2,&infoB);

         // L1 cache     
         block3d_info_t info3d = {M2,N2,K2,MYBLAS_TILE_M,MYBLAS_TILE_N,MYBLAS_TILE_K};
         myblas_dgemm_kernel(alpha,A2,B2,C+ldc*j2+i2,ldc,&info3d);

       }
     }
   }

}}}

free(A2);
free(B2);

この変更によって、計算速度は次のようになりました。

＜L2ループ交換前＞

Max  Peak MFlops per Core: 52800 MFlops 
Base Peak MFlops per Core: 46400 MFlops 
size  , elapsed time[s],          MFlops,   base ratio[%],    max ratio[%] 
   16,     2.00272E-05,         447.392,        0.964208,        0.847334 
   32,     0.000152111,         451.039,        0.972067,         0.85424 
   64,     0.000251055,         2137.29,         4.60622,         4.04789 
  128,      0.00179315,         2366.49,         5.10019,         4.48198 
  256,       0.0152731,         2209.84,         4.76258,          4.1853 
  512,        0.130602,         2061.39,         4.44265,         3.90415 
 1024,        0.916094,         2347.61,          5.0595,         4.44623 
 2048,         7.26869,         2365.28,         5.09758,         4.47969 

＜L2ループ交換後＞

Max  Peak MFlops per Core: 52800 MFlops 
Base Peak MFlops per Core: 46400 MFlops 
size  , elapsed time[s],          MFlops,   base ratio[%],    max ratio[%] 
   16,     2.00272E-05,         447.392,        0.964208,        0.847334 
   32,     0.000138044,             497,         1.07112,        0.941287 
   64,     0.000266075,         2016.63,         4.34619,         3.81938 
  128,      0.00193405,         2194.08,         4.72861,         4.15545 
  256,       0.0143108,         2358.43,         5.08281,         4.46671 
  512,         0.11483,         2344.53,         5.05286,         4.44039 
 1024,        0.893248,         2407.65,          5.1889,         4.55995 
 2048,         7.20221,         2387.11,         5.14463,         4.52104 

結果としては、行列サイズ2048×2048の場合、ループ交換前後で比較すると、2365 MFLOPS → 2387 MFLOPS、基本周波数の理論ピーク性能比で5.10% → 5.14%と、わずかに速度が上がりました。

次の記事

元の記事はこちらです。

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

jyam45/dgemm