ループ交換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で公開しています。