行列行列積の高速化手法

本記事では、行列行列積に焦点を当てて、高速化を行う上で重要な事項に関して、既存の資料よりもできるだけわかりやすく説明することを心がけて記載する。
アーキテクチャの詳細など細かい技術的なことが知りたい人はHPCで検索すれば溢れかえるほどの資料が出てくるはずなので、ここでは言及しない。

*内容が間違っている場合は指摘いただけると嬉しいです。
追記：ループアンローリングでは、最外側のループからアンロールを行なった方が高速なので、プログラムを変更しました。

#序論
コンピュータの高速化を行う際にどのような計算に対して高速な処理を行うのかというのは、いろいろな議論が起こっている。

TOP500は世界のコンピュータの中で高速なコンピュータの上位500位までをランク付けするもので、年に2回行われている。

かつて、日本の富士通が開発したK computerが世界1位を取ったということでニュースになったのもこのTOP500によるランク付けでの1位ということである。

このTOP500による性能指標にはLINPACKというベンチマークソフトが利用されている。
LINPACKでは密行列の連立1次方程式を計算している。

#行列行列積
行列行列積は密行列と密行列の積を計算することをさしており、
通常、連立1次方程式よりもキャッシュの再利用効率が良いので、ピーク性能に近い性能が出しやすいことから性能評価に利用される。

今回は、LINPACKの密行列の連立1次方程式よりもさらに性能が出やすい行列行列積を利用して性能を評価してみる。

#行列行列積の最適化
行列行列積では単純に計算部分を示すと以下のようになる。

行列行列積計算部分.c

for (i = 0; i < N; i++){
    for (j = 0; j < N; j++){
        for (k = 0; k < N; k++){
            c[i*N+j]+=a[i*N+k]*b[k*N+j];
        }
    }
}

上記のコードでも行列行列積がしたいだけなら十分なのだが、これだけでは実行効率という意味では効率化は何も行われていないことになる。

##i, j, kループ交換法
for文のループ変数の順番を交換する方法
これを利用することで、メモリのアクセス方向を連続にすることができる。キャッシュの空間的局所性を利用している。
ポイントとしては、C言語などでは行方向にメモリが連続なので行方向の連続性を保持できる方法を考える必要がある。
転置なしの行列行列積である場合は、i,k,jループが最もアクセスが連続になるループになる。

行列行列積ループ交換法適用.c

for (i = 0; i < N; i++){
    for (k = 0; k < N; k++){
        for (j = 0; j < N; j++){
            c[i*N+j]+=a[i*N+k]*b[k*N+j];
        }
    }
}

##ループアンローリング
ループアンローリングでは、アンロールという名前の通り、巻かれたものを広げるという意味である。
for文などのループを行う際には分岐命令という命令が発行されているが、分岐命令は普通の命令に比べて処理が多く、時間がかかる。
高速なプログラムではできるだけループを減らす工夫がされている。
今回のプログラムで分岐命令を減少させるには積和演算部分を数段に展開することである。
ループアンローリングをiループに対して4段・kループに対して4段、
実施したプログラムを以下に示す。

行列行列積ループアンローリング適用.c

for (i = 0; i < N; i+=4){
    for (k = 0; k < N; k+=4){
        for (j = 0; j < N; j++){
            c[i*N+j]+=a[i*N+k]*b[k*N+j];
            c[i*N+j]+=a[i*N+(k+1)]*b[(k+1)*N+j];
            c[i*N+j]+=a[i*N+(k+2)]*b[(k+2)*N+j];
            c[i*N+j]+=a[i*N+(k+3)]*b[(k+3)*N+j];

            c[(i+1)*N+j]+=a[(i+1)*N+k]*b[k*N+j];
            c[(i+1)*N+j]+=a[(i+1)*N+(k+1)]*b[(k+1)*N+j];
            c[(i+1)*N+j]+=a[(i+1)*N+(k+2)]*b[(k+2)*N+j];
            c[(i+1)*N+j]+=a[(i+1)*N+(k+3)]*b[(k+3)*N+j];

            c[(i+2)*N+j]+=a[(i+2)*N+k]*b[k*N+j];
            c[(i+2)*N+j]+=a[(i+2)*N+(k+1)]*b[(k+1)*N+j];
            c[(i+2)*N+j]+=a[(i+2)*N+(k+2)]*b[(k+2)*N+j];
            c[(i+2)*N+j]+=a[(i+2)*N+(k+3)]*b[(k+3)*N+j];

            c[(i+3)*N+j]+=a[(i+3)*N+k]*b[k*N+j];
            c[(i+3)*N+j]+=a[(i+3)*N+(k+1)]*b[(k+1)*N+j];
            c[(i+3)*N+j]+=a[(i+3)*N+(k+2)]*b[(k+2)*N+j];
            c[(i+3)*N+j]+=a[(i+3)*N+(k+3)]*b[(k+3)*N+j];
        }
    }
}

##ループアンローリングを適用した場合の応用
ループアンローリングでは、上記のプログラムのように演算部分を展開するが、展開した段数の加算に無駄な演算が行われているので、その部分をループ内で1回のみ行うようにすると計算量が減少する。
ループアンローリングを適用した場合の応用を行なったプログラムを以下に示す。

行列行列積ループアンローリング適用の応用.c

for (i = 0; i < N; i+=4){
    i1 = i + 1;
    i2 = i + 2;
    i3 = i + 3;
    for (k = 0; k < N; k+=4){
        k1 = k + 1;
        k2 = k + 2;
        k3 = k + 3;

        for (j = 0; j < N; j++){
            c[i*N+j]+=a[i*N+k]*b[k*N+j];
            c[i*N+j]+=a[i*N+k1]*b[k1*N+j];
            c[i*N+j]+=a[i*N+k2]*b[k2*N+j];
            c[i*N+j]+=a[i*N+k3]*b[k3*N+j];

            c[i1*N+j]+=a[i1*N+k]*b[k*N+j];
            c[i1*N+j]+=a[i1*N+k1]*b[k1*N+j];
            c[i1*N+j]+=a[i1*N+k2]*b[k2*N+j];
            c[i1*N+j]+=a[i1*N+k3]*b[k3*N+j];

            c[i2*N+j]+=a[i2*N+k]*b[k*N+j];
            c[i2*N+j]+=a[i2*N+k1]*b[k1*N+j];
            c[i2*N+j]+=a[i2*N+k2]*b[k2*N+j];
            c[i2*N+j]+=a[i2*N+k3]*b[k3*N+j];

            c[i3*N+j]+=a[i3*N+k]*b[k*N+j];
            c[i3*N+j]+=a[i3*N+k1]*b[k1*N+j];
            c[i3*N+j]+=a[i3*N+k2]*b[k2*N+j];
            c[i3*N+j]+=a[i3*N+k3]*b[k3*N+j];
        }
    }
}

##キャッシュブロッキング
キャッシュブロッキングはキャッシュの時間的局所性をカバーできる非常に効率的な高速化手法である。キャッシュブロッキングによって実行効率を向上させるためには、キャッシュサイズを意識する必要がある。
キャッシュに収まるデータサイズは各CPU毎に異なるため、CPUに応じたキャッシュサイズに収まるようにブロッキングする。

以下にBLOCKサイズが10の場合 (実際には10[ブロックの大きさ]*8[double1つ8Byte]16[4段4段のアンロール段数]=1ブロック1280Byte) でのプログラムを示す。

行列行列積キャッシュブロッキング適用.c

BLOCK = 10;

for (i = 0; i <N; i+=BLOCK){
    for(k = 0; k < N; k+=BLOCK){
        for(j = 0; j < N; j+=BLOCK){
            for(ii = i; ii < (i + BLOCK); ii+=4){
                ii1 = ii + 1;
                ii2 = ii + 2;
                ii3 = ii + 3;
                for (kk = k; kk < (k + BLOCK); kk+=4){
                    kk1 = kk + 1;
                    kk2 = kk + 2;
                    kk3 = kk + 3;

                    for (jj = j; jj < (j + BLOCK); jj++){
                        c[ii*N+jj]+=a[ii*N+kk]*b[kk*N+jj];
                        c[ii*N+jj]+=a[ii*N+kk1]*b[kk1*N+jj];
                        c[ii*N+jj]+=a[ii*N+kk2]*b[kk2*N+jj];
                        c[ii*N+jj]+=a[ii*N+kk3]*b[kk3*N+jj];

                        c[ii1*N+jj]+=a[ii1*N+kk]*b[kk*N+jj];
                        c[ii1*N+jj]+=a[ii1*N+kk1]*b[kk1*N+jj];
                        c[ii1*N+jj]+=a[ii1*N+kk2]*b[kk2*N+jj];
                        c[ii1*N+jj]+=a[ii1*N+kk3]*b[kk3*N+jj];

                        c[ii2*N+jj]+=a[ii2*N+kk]*b[kk*N+jj];
                        c[ii2*N+jj]+=a[ii2*N+kk1]*b[kk1*N+jj];
                        c[ii2*N+jj]+=a[ii2*N+kk2]*b[kk2*N+jj];
                        c[ii2*N+jj]+=a[ii2*N+kk3]*b[kk3*N+jj];

                        c[ii3*N+jj]+=a[ii3*N+kk]*b[kk*N+jj];
                        c[ii3*N+jj]+=a[ii3*N+kk1]*b[kk1*N+jj];
                        c[ii3*N+jj]+=a[ii3*N+kk2]*b[kk2*N+jj];
                        c[ii3*N+jj]+=a[ii3*N+kk3]*b[kk3*N+jj];
                    }
                }
            }
        }
    }
}

##転置処理
i, j, kループ交換法の説明を行なったが、i, j, kループ交換を行わなくても、片方の行列の行と列を入れ替えることによって、アクセス方向を連続にすることができる。
この場合は、演算部分を少し工夫する必要があり、その後のブロッキングなどにも影響してくるが、そこまでの高速化に影響がない(アクセス方向を変えるだけでも転置と同様の効果が出る)ので、今回は取り上げない・・・時間があったら追記するかも！

#まとめ
・i, j, kループ交換法
・ループアンローリング
・キャッシュブロッキング
・転置処理

これら上記で説明した最適化を元に基本的には高速化を行なっていくことが通常は多くなっている。(むしろ、他の効率的な高速化があれば逆に教えていただきたい・・・)
アンロール段数やキャッシュブロッキングではチューニングによって実行効率が大きく向上するため、いろいろなパラメータを指定してみると良い。

今回取り上げた内容は検索や以下の参考文献からより詳しく解説された資料を見ることができる。

#参考
・TOP500-wikipedia
・LINPACK-wikipedia
・TOP500のランキングルールが変わる? - LINPACKに変わるベンチマークが提案
・top500
・第1講 プログラム高速化の基礎