行列乗算高速化プロジェクト

2002/1/10 改訂

BLASやATLASなどのライブラリを使うと，行列乗算やLU分解などの計算が，自分で組んだプログラムより数段速いことがよく知られている。それは，それらのライブラリではメモリやキャッシュの管理が高度になされているからであり，また，ループのアンローリングによるオーバーヘッドの低減やパイプライニングによるスケジューリングの最適化がなされているからであろう。

しかし，自分で組んだものが遅い，というのは悔しいものがあるので，なんとかしてチューニングできないものだろうかと思い立ち，いろいろやってみようと思う。

行列乗算のテスト

ベンチマークによく使われる（と思う）のが，行列乗算である。これは，２つの行列Ａ,Ｂの掛け算Ｃ＝ＡＢを計算するもので，目的は単純だが，この計算を高速化しようといろいろなところで日々努力されている（はず）。行列乗算を普通に計算する手順をFortranで書くと，以下のようなプログラム（サブルーチン）になる。ここでは簡単のため，n×n の正方行列を考えることにする。 

      subroutine multmm(n, A, lda, B, ldb, C, ldc)
      implicit real*8 (a-h, o-z)
      dimension A(lda,*), B(ldb,*), C(ldc,*)
C
      do j=1,n
        do i=1,n
          C(i,j) = 0.0d0
        end do
      end do
C
      do i=1,n
        do j=1,n
          do k=1,n
            C(i,j) = C(i,j) + A(i,k)*B(k,j)
          end do
        end do
      end do
C
      return
      end
C言語では，以下のよう。 

void multmm(int n, double A[][lda], double B[][ldb], double C[][ldc])
{
	int i, j, k;

	for(i=0; i<n; i++){
		for(j=0; j<n; j++){
			C[i][j] = 0.0;
		}
	}

	for(i=0; i<n; i++){
		for(j=0; j<n; j++){
			for(k=0; k<n; k++){
				C[i][j] += A[i][k]*B[k][j];
			}
		}
	}
}
MATLABなどでは，以下のように簡潔である。 
>> C=A*B;

計算環境を以下に示す。

＜計算環境＞
PC : PC/AT互換機
CPU : Intel Pentium III 1066MHz
L1 Cache : 32KByte
L2 Cache : 512KByte
Main Memory : 640MByte
Compiler : DIGITAL Visual Fortran 6.0
まず，参考までに，MATLAB 6での実行結果を示しておく。

>> n=1000;
>> A=rand(n); B=rand(n);
>> tic; C=A*B; toc
n = 1000計算時間(秒)MFLOPS
C=A*B2.88693

流石に，ATLASで最適化されているだけあって高速である。これを目標としたい。

Fortranによるチューニング

さて，自分で手を加えるときに，最もやり易そうなのはFortranだと思ったので，Fortranでやってみた。

＜計算順序の見直し＞
プログラムをチューニングする上で，まず気をつけるのはメモリの参照方向であろう。上記のFortranのプログラムを見ると，ループの順番が外側から(i,j,k)の順になっている。この場合，順番は可変であるので，(j,i,k)や(k,j,i)などと順番を変更してみたところ，以下のような結果になった。尚，特別な記述がない限り，コンパイルオプションは「Full Optimizations（Releaseモードのデフォルト）」で実行してある。

n = 1000計算時間(秒)MFLOPS
(i,j,k)13.23151
(i,k,j)33.4559
(j,i,k)24.4675
(j,k,i)9.76204
(k,i,j)64.9730
(k,j,i)38.4452

この中では，(j,k,i)の順番が一番速い。

＜計算の内積化＞
まず，上記のプログラム（Fortran）では， C(i,j) = C(i,j) + A(i,k)*B(k,j) のように，C(i,j)を毎回ロード，ストアしているので，この部分を以下のように内積系計算にすることで，ロードおよびストアの回数を低減する。 

      do i=1,n
        do j=1,n
          s = 0.0d0
          do k=1,n
            s = s + A(i,k)*B(k,j)
          end do
          C(i,j) = s
        end do
      end do
(i,j,k)と(j,i,k)は，そのまま内積計算のように変更できる。実行結果は以下のようになった

n = 1000計算時間(秒)MFLOPS
(i,j,k)13.23151
(i,j,k)+内積化12.23163
(j,i,k)24.4675
(j,i,k)+内積化25.9277

若干，速くなった。

＜内積化におけるワークベクトルの利用＞
(i,j,k)型の計算では，内積化において以下のようにワークベクトルを１本使うことによって，メモリの非連続方向のデータ参照回数が減少し，さらに高速化が狙える。 

      do i=1,n
        do k=1,n
          W(k) = A(i,k)
        end do
        do j=1,n
          s = 0.0d0
          do k=1,n
            s = s + W(k)*B(k,j)
          end do
          C(i,j) = s
        end do
      end do
n = 1000計算時間(秒)MFLOPS
(i,j,k)13.23151
+内積化12.23163
+内積化+ワーク10.21195

とりあえず，この程度までは簡単に改善できることがわかる。

＜ループのアンローリング＞
チューニング技法のひとつにループのアンローリング（unrolling）がある。まず，アンローリングの効果は以下のようになる。 方法は，以下のよう。 

      do i=1,n
        do j=1,n
          s0 = 0.0d0
          s1 = 0.0d0
          s2 = 0.0d0
          s3 = 0.0d0
          do k=1,n,4
            s0 = s0 + A(i,k)*B(k,j)
            s1 = s1 + A(i,k+1)*B(k+1,j)
            s2 = s2 + A(i,k+2)*B(k+2,j)
            s3 = s3 + A(i,k+3)*B(k+3,j)
          end do
          C(i,j) = s0 + s1 + s2 + s3
        end do
      end do
この場合では，本当は，nが4の倍数でないときの処理を加えなければならないが，ここでは省略する。結果を以下に示す。

n = 1000計算時間(秒)MFLOPS
(i,j,k)13.23151
+内積化12.23163
+内積化+ワーク10.21195
+内積化+ワーク+unrolling(2)10.29194
+内積化+ワーク+unrolling(4)10.33193

現在のPC上ではほとんど結果は変わらないようなので，アンローリングは不要のようである。

＜コンパイラのコンパイルオプション＞
あまり効果が伸びなかったので，コンパイルオプションをFull Optimizationsから最高レベルの最適化であるMaximum Optimizationsに変えたところ，脅威のMFLOPS値が叩き出された。

n = 1000計算時間(秒)MFLOPS
(i,j,k), Full Opt.13.23151
(i,j,k), Max. Opt.5.16387

たったこれだけで，あっという間に２倍以上の高速化。今までの努力は何だったのだろうか・・。因みに，Max. Opt. では，上で示したチューニングを行ってもほとんど無効であった（逆に，遅くなる）。

（教訓） 自分で本格的に手を加える前に，まずコンパイルオプションを見直してみよう。

一体，Max. Opt.でコンパイラが何をしてくれたのかは定かではないが，キャッシュを利用したとしか考えられないので，そっちがその気なら・・，という意気込みで次へ。

＜行列のブロック化（blocking）＞
さて，このままでは悔しいので，通常では最高のコンパイルオプションである「Full Optimizations」において，なんとか「Maximum Optimizations」のMFLOPS値を上回ることができないものかと考えた。そのためには，キャッシュメモリ（特に，L2 Cache）を有効に利用することを考えなければならないだろう。そこで，まず，どれくらいの行列サイズまでキャッシュが効くのかをテストした。その結果（MFLOPS）を以下に示す。

type / n 20406080100120140160180200220240260280300
(i,j,k) 307363356368377381378304364338321213225195192
(i,j,k)+内積化 499444480491500491507453485459425258257224215
(i,j,k)+内積化+ワーク400444480491487491488483494481442294265246219
(i,k,j) 2853203453073573512711852241931891111079482
(j,i,k) 363347356359384388384338373338279143138127103
(j,i,k)+内積化 363470502526540526548499511481370162152140103
(j,k,i) 333400425433444451448446448439420289263265226
(k,i,j) 3073473353073503402711872271952001411087458
(k,j,i) 33342144244645445144845345545442528918911990
（単位：MFLOPS）

この結果から，行列サイズ 200×200（n = 200）くらいまでキャッシュが効いていることがわかる。そこで，行列Ａ,Ｂ共に10×10〜100×100のブロックに分割し，Ｃ＝ＡＢを計算する。ブロック行列同士の乗算には，キャッシュが有効なときに最も高速な結果となった「(j,i,k)+内積化」タイプを使用する。

（n = 1000）
ブロックサイズ計算時間(秒)MFLOPS
(10) 6.01332
(20) 5.00400
(30) 5.54361
(40) 4.95404
(50) 4.91407
(60) 5.47365
(70) 6.03331
(80) 6.94288
(90) 7.48267
(100)7.71259

これで，ようやく「何も手を加えない状態のMaximum Optimizations」に追いつくことができた。やれやれ，である。尚，Ｃ言語でも同様の結果が得られた（テストに使用したコンパイラはMicrosoft Visual C++ 6.0，もちろん「Maximum Optimizations」は存在しない）。

整理したプログラムのソースはここ（C言語用）とここ（Fortran用）で公開中。

＜Linuxでの実行結果＞
上記の作業をLinux上でも実行し（今までと同じPCを使用），さらにBLASのサブルーチンDGEMM（ATLAS最適化）での実行結果も示す。尚，コンパイル方法や実行方法などは以下の通り。
% g77 -O3 multmm.f xclock_ux.c -latlas
% ./a.out

n = 1000計算時間(秒)MFLOPS
(i,j,k) 12.95154
(i,j,k)+内積化 12.40161
(i,j,k)+内積化+ワーク9.89202
blocking(50) 5.30377
ATLAS最適化BLAS/DGEMM2.48806

BLASのDGEMMが圧倒的に高速である。このテストに使用したCPUはPentiumIII 1066MHzであるので，DGEMMは理論ピーク性能の実に80％を引き出していることになる。

Go back