行列積とループ・ブロッキング。

インテル日本語資料の「IA-32 インテルアーキテクチャー最適化」にも記されている、ループ・ブロッキングの有効性を検証してみた。
ループ・ブロッキングとは、キャッシュサイズを考慮してループをブロックすることである。
反復のたびにおきてしまうキャッシュミスを軽減することができ、キャッシュとメモリーのパフォーマンスを向上させることができる。

実験環境は Xeon 3.0GHz であるが、今回はスレッド数は１に固定した。ちなみにキャッシュサイズは 4096 KB となっている。
今回は 512 x 512 と 1024 x 1024 の２種類のみであるが、次の４パターンについてプロファイルしてみた。
(0)　GotoBLAS （あくまでも参考）
(1)　通常の積演算
(2)　通常の積演算にループ・ブロッキングを適用
(3)　事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算
(4)　事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算に対してループ・ブロッキングを適用


ループの大きさによってかなりのバラツキがあることにおどろいた。
通常のループであれば、ブロックサイズは２～８程度と小さめで、
事前にアクセス方向を修正しておけば、２５６と大きめがよいようだ。


512 x 512
　　　　1　　　2　　4　　 8　　　16　　32　　64　 　128　 256 　　512
(0)　　35
(1)　1700
(2)　1500　765　825　1000　1327　1400　1420　1450　1450　1500
(3)　 190
(4) 　920　365　257　　228　　211　200　　196　　202　198　　196


1024 x 1024
　　　　1　　　　2　　　　4　　　8　　　16　　　32　　　64　 　128　 　256 　　512　　1024
(0) 235
(1) 60000
(2) 65000　23300　10700　9600　11200　11500　11500　11800　11700　17500　54500
(3) 2550
(4) 7600　　3050　　2120　1800　　1710　1625　　1580　　1626　1610　　2050 　2420


以下はソースコードの一部である。


　/* (1)通常の積演算 */
　for (m = 0; m ＜ M; m++) {
　　for (n = 0; n ＜ N; n++) {
　　　for (k = 0; k ＜ K; k++) {
　　　　C[m * N + n] += A[m * K + k] * B[k * N + n];
　　　}
　　}
　}


　/* (2)通常の積演算にループ・ブロッキングを適用 */
　for (m = 0; m ＜ M; m+=block) {
　　for (n = 0; n ＜ N; n+=block) {
　　　for (k = 0; k ＜ K; k+=block) {
　　　　
　　　　for (mm = m; mm ＜ m+block; mm++) {
　　　　　for (nn = n; nn ＜ n+block; nn++) {
　　　　　　for (kk = k; kk ＜ k+block; kk++) {
　　　　　　　C[mm * N + nn] += A[mm * K + kk] * B[kk * N + nn];
　　　　　　}
　　　　　}
　　　　}
　　　　
　　　}
　　}
　}

　/* (3)事前に行列Ｂに対して転置処理を行い、アクセス方向を修正した積演算 */
　for (m = 0; m ＜ M; m++) {
　　for (n = 0; n ＜ N; n++) {
　　　for (k = 0; k ＜ K; k++) {
　　　　C[m * N + n] += A[m * K + k] * B[n * K + k];
　　　}
　　}
　}

　/* (4)事前に行列Ｂに対して転置処理を行い、
　　　アクセス方向を修正した積演算に対してループ・ブロッキングを適用 */
　for (m = 0; m ＜ M; m+=block) {
　　for (n = 0; n ＜ N; n+=block) {
　　　for (k = 0; k ＜ K; k+=block) {
　　　　
　　　　for (mm = m; mm ＜ m+block; mm++) {
　　　　　for (nn = n; nn ＜ n+block; nn++) {
　　　　　　for (kk = k; kk ＜ k+block; kk++) {
　　　　　　　C[mm * N + nn] += A[mm * K + kk] * B[nn * K + kk];
　　　　　　}
　　　　　}
　　　　}
　　　　
　　　}
　　}
　}