数値計算レポート 第２回

g99p1457 矢島伸吾

出題日： 2001 年 5 月 18 日提出日： 2001 年 6 月 2 日提出期限： 2001 年 6 月 2 日

1 問題

1000x1000 行列 (double) の積を高速に計算するプログラムを作れ。 10 回別の行列のペアに適用して，平均の計算時間を出せ。

2 プログラムの概要

この問題を解くにあたり、以下のような関数を作成した。

mul() 行列の積を計算する関数（メイン）

multest() 行列の積を基本通りに計算する関数（デバッグ用）

init() 行列の初期化をする関数

check() mul() から multest() を引き、計算ミスを調べる関数

また、 N(= 1000) を表す為に、 M, L を定義し、 L ∗ M = N となるようにした。デバッグ時にはマクロにより N を小さくするとともに、 multest(), check() を用いて計算にミスがないかどうか確かめられるようにしてある。

init() は、行列用の配列 a,b,c,d （N ∗ N 行列）の初期化を行う。 a,b は全要素にランダムなデータを入力し、 c,d は全要素を 0.0 で初期化している。 c, d は a ∗ b の格納に使われる。

multest() は、なにも高速化を行わず、基本通りに行列計算を行う関数である。 check() により、高速化した関数の結果からこの関数の結果を要素毎に引き、その絶対値の最も大きなものを返す。この値が十分に小さければ、計算が正しいとする事ができる。今回最終的に作ったプログラムでは、 3 ∗ 10^{−10} 程度の誤差が出る。

main() 関数では、まず行列の初期化 init() を実行している。その後に clock() 関数を用いて現在のクロックカウントを得て、 mul() により行列の積を計算し、 clock() で計算後のクロックカウントを得ている。その２数の差を得、 CLOCKS PER SEC で割ってやる事により、計算にかかった時間を求めている。 for 文によりそれを１０回実行し、計算にかかった時間の和を１０で割ることにより平均の時間を求めている。また、マクロにより、 DEBUG を定義した時には計算に誤りがないかを確認できるようにしてある。

3 行列の積の計算方法について

3.1 基本

行列の式の基本的なアルゴリズムは、次のような物である。

    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            for(k=0; k<N; k++){
                d[i][j] += a[i][k]*b[k][j];
            }
        }
    }
ここで、 a[N][N], b[N][N], d[N][N] であり、 d が結果となる。

multest() 内ではこのアルゴリズムを用いて正解を求めている。

3.2 キャッシュの利用

ＣＰＵにはキャッシュメモリが組み込まれており、ここに一度メモリから読み込んでおいた値を保存しておき、次にその値が必要になった時メモリからではなくキャッシュから読み込む事で高速化する事が出来る。キャッシュ内のデータは古い物から消されて行ってしまうので、行列 a,b の積を求めるにあたってキャッシュを有効利用する為には、行列 a,b がキャッシュ内に収まりきるような十分小さくなければならない。

1000 ∗ 1000 行列そのままでは大きすぎるので、

A∗B = ( \begin{matrix}A_{11} & … & A_{1L} \\ ⋮ & ... & ⋮ \\ A_{L1} & … & A_{LL} \end{matrix} ) ( \begin{matrix}B_{11} & … & B_{1L} \\ ⋮ & ... & ⋮ \\ B_{L1} & … & B_{LL} \end{matrix} ) = ( \begin{matrix}A_{11}B_{11} + … + A_{1L}B_{L1} & … & A_{1L}B_{1L} + … + A_{1L}B_{LL} \\ ⋮ & ... & ⋮ \\ A_{L1}B_{11} + … + A_{LL}B_{L1} & … & A_{L1}B_{1L} + … + A_{LL}B_{LL} \end{matrix} )

を利用して、小さな行列の積に直す。このプログラムでは、 A_{IJ} を M ∗ M 行列とし、 L ∗ M = N(= 1000) となるように設定してある。

ここで、 I′ = (I − 1)M, J′ = (J − 1)M とおくと、

A_{IJ} = ( \begin{matrix}a_{I′+1,J′+1} & … & A_{I′+1,J′+M} \\ ⋮ & ... & ⋮ \\ A_{I′+M,J′+1} & … & A_{I′+M,J′+M} \end{matrix} )

と表される。これをプログラムで表すと、次の様になる。

    for(I=0; I<N; I+=M){
        for(K=0; K<N; K+=M){
            for(J=0; J<N; J+=M){
                for(i=0; i<M; i++){

                    for(j=0; j<M; j++){
                        for(k=0; k<M; k++){
                            c[I+i][J+j] += a[I+i][K+k]*b[K+k][J+j];
                        }
                    }
                }
            }
        }
    }

外側の３ループで

C_{IJ}+ = A_{IK} ∗ B_{KJ}

を指定し、内側の３ループでそれぞれの行列の積を求めている。

3.3 内積化

上のプログラムでは、行列 c の要素についても毎回ロード、ストアを行っている（キャッシュに読みこんでいる）ので、 a,b へのキャッシュが効きにくく、速さも損なわれる。

よって、一時変数 s に最内部の総和を保存し、計算が終了してから c に加える事で、ロード、ストアを少なくし、キャッシュの使用量も少なくしている。

プログラムは、上記のアルゴリズムの最内部を以下の様に変更する。

                        s = 0.0;
                        for(k=0; k<M; k++){
                            s += a[I+i][K+k]*b[K+k][J+j];
                        }
                        c[I+i][J+j] += s;

3.4 計算順序

I,J,K のループの順序を変え、どれが最も高速か調べた。学校のパソコン、デバッグモード（N = 400）で、結果は 0.55〜0.58(sec) とばらつき、一番速かったのは (I,J,K),(J,I,K) の順だった。よってここでは (I,J,K) の順序を採用する事にする。

3.5 最適なキャッシュの値

今回は L*M=1000 にならなければならないので、候補として上げられる M の値は M = 20, 25, 40, 50 であろう。どの値が最も適しているかを確かめてみた。使用したマシンは自宅の Celeron600MHz である。

M 25 40 50
計算時間 (sec) 12.13 10.38 28.61

これより、 M = 40 が最も適している事がわかった。

3.6 ループのアンローリング

最内部の for ループを書き下す事によって、分岐処理の回数を減らす事が出来、高速化できる。今回は M = 40 としてあるので、 k ループの計算 40 個を展開し、全て書き下した。

このプログラムを最終版とし、レポート末に記載する。

3.7 その他の高速化

1 回目の実行では、若干計算速度が遅くなる傾向があった。これは、最初の実行時は命令がキャッシュに入っておらず、メモリに読みこまれる為であると思われるので、まず 1 回実行した後で、 10 回の計測を行っている。

また、 mul() を関数にせず、 main() の中にいれてしまえば関数呼出にかかる時間を削る事が出来るが、ほとんど計測不能なほどしか変わらないだろう。

4 実行結果

実行環境は、自宅の自作パソコン（CPU:Celeron 600MHz, memory:256MB) である。

4.1 基本

63.99 sec
66.02 sec
66.24 sec
64.82 sec
64.70 sec
64.43 sec
65.25 sec
65.85 sec
64.37 sec
65.09 sec
average: 65.08

4.2 キャッシュの利用

17.20 sec
17.69 sec
17.24 sec
17.35 sec
18.12 sec
18.12 sec
17.96 sec
17.52 sec

17.41 sec
17.41 sec
average: 17.60

4.3 ＋内積化

10.44 sec
10.33 sec
10.44 sec
10.60 sec
10.49 sec
10.38 sec
10.32 sec
10.27 sec
10.27 sec
10.27 sec
average: 10.38

4.4 ＋アンローリング

7.14 sec
7.14 sec
7.19 sec
7.20 sec
7.19 sec
7.20 sec
7.19 sec
7.14 sec
7.20 sec
7.41 sec
average: 7.20

4.5 結果

7.20(sec) で計算する事が出来た。（時によりもう少し速いこともある）

4.6 補足

学校のパソコン (CPU:Intel PentiumIII 600MHz, memory:256MB) での、＋内積化での実行結果は 9.10sec だった。周波数では同じでも実行結果には開きがあることがわかる。

また、最適化していないアンローリングプログラムを学校で実行した所、 6.44(sec) という値が出た。

4.7 ＋アンローリング（デバッグモード）での実行結果

DEBUG mode N = 400
0.50 sec (max_dif:3.055902e-010)
0.44 sec (max_dif:2.619345e-010)
0.44 sec (max_dif:2.473826e-010)
0.49 sec (max_dif:2.764864e-010)
0.49 sec (max_dif:2.910383e-010)
0.44 sec (max_dif:3.201421e-010)
0.50 sec (max_dif:2.910383e-010)
0.44 sec (max_dif:2.619345e-010)
0.50 sec (max_dif:2.619345e-010)
0.44 sec (max_dif:2.910383e-010)
average: 0.47

これにより、ほぼ正確に行列計算ができている事がわかる。

5 プログラム

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
//#define DEBUG
#ifdef DEBUG
#define L 10
#define M 40
#else
#define L 25
#define M 40
#endif
#define N L*M
double a[N][N], b[N][N], c[N][N], d[N][N];
void init(double a[N][N], double b[N][N]);
double check(double c[N][N], double d[N][N]);
void multest(double a[N][N], double b[N][N], double d[N][N]);
void mul(double a[N][N], double b[N][N], double c[N][N]);

void mul(double a[N][N], double b[N][N], double c[N][N])
{
    int I, J, K, i, j;
    double s;
    for(I=0; I<N; I+=M){
        for(J=0; J<N; J+=M){
            for(K=0; K<N; K+=M){
                for(i=0; i<M; i++){
                    for(j=0; j<M; j++){
                        s = 0.0;
                        s += a[I+i][K]*b[K][J+j];
                        s += a[I+i][K+1]*b[K+1][J+j];
                        s += a[I+i][K+2]*b[K+2][J+j];
                        s += a[I+i][K+3]*b[K+3][J+j];
                        s += a[I+i][K+4]*b[K+4][J+j];

s += a[I+i][K+5]*b[K+5][J+j];
s += a[I+i][K+6]*b[K+6][J+j];
s += a[I+i][K+7]*b[K+7][J+j];
s += a[I+i][K+8]*b[K+8][J+j];
s += a[I+i][K+9]*b[K+9][J+j];
s += a[I+i][K+10]*b[K+10][J+j];
s += a[I+i][K+11]*b[K+11][J+j];
s += a[I+i][K+12]*b[K+12][J+j];
s += a[I+i][K+13]*b[K+13][J+j];
s += a[I+i][K+14]*b[K+14][J+j];
s += a[I+i][K+15]*b[K+15][J+j];
s += a[I+i][K+16]*b[K+16][J+j];
s += a[I+i][K+17]*b[K+17][J+j];
s += a[I+i][K+18]*b[K+18][J+j];
s += a[I+i][K+19]*b[K+19][J+j];
s += a[I+i][K+20]*b[K+20][J+j];
s += a[I+i][K+21]*b[K+21][J+j];
s += a[I+i][K+22]*b[K+22][J+j];
s += a[I+i][K+23]*b[K+23][J+j];
s += a[I+i][K+24]*b[K+24][J+j];
s += a[I+i][K+25]*b[K+25][J+j];
s += a[I+i][K+26]*b[K+26][J+j];
s += a[I+i][K+27]*b[K+27][J+j];
s += a[I+i][K+28]*b[K+28][J+j];
s += a[I+i][K+29]*b[K+29][J+j];
s += a[I+i][K+30]*b[K+30][J+j];
s += a[I+i][K+31]*b[K+31][J+j];
s += a[I+i][K+32]*b[K+32][J+j];
s += a[I+i][K+33]*b[K+33][J+j];
s += a[I+i][K+34]*b[K+34][J+j];
s += a[I+i][K+35]*b[K+35][J+j];
s += a[I+i][K+36]*b[K+36][J+j];
s += a[I+i][K+37]*b[K+37][J+j];
s += a[I+i][K+38]*b[K+38][J+j];
s += a[I+i][K+39]*b[K+39][J+j];
c[I+i][J+j] += s;
}
}
}
}
}
}
int main()
{
int i;
double t, sum = 0.0;
clock_t start, end;
#ifdef DEBUG
printf("DEBUG mode N = %d\n",N);
#endif
init(a, b);

mul(a, b, c); // １回処理を実行して命令等をキャッシュに

for(i = 0; i < 10; i++){
init(a, b);
start = clock();
mul(a, b, c);
end = clock();
t = (double)(end - start)/CLOCKS_PER_SEC;
#ifdef DEBUG
multest(a, b, d);   // 基本のやり方でやる
printf("%.2f sec (max_dif:%e)\n", t, check(c, d));
#else
printf("%.2f sec\n", t);
#endif
sum += t;
}
printf("average: %.2f\n", sum / 10);
return 0;
}
void multest(double a[N][N], double b[N][N], double d[N][N])
{
int i, j, k;
for(i=0; i<N; i++){

    for(j=0; j<N; j++){
        for(k=0; k<N; k++){
            d[i][j] += a[i][k]*b[k][j];
        }
    }
}
}
void init(double a[N][N], double b[N][N])
{
int i, j;
srand( (unsigned)time( NULL ) );
for(i = 0; i < N; i++){
    for(j = 0; j < N; j++){
        a[i][j] = (double)rand() / 1000.0;
        b[i][j] = (double)rand() / 1000.0;
        c[i][j] = 0.0;
        d[i][j] = 0.0;
    }
}
}
double check(double c[N][N], double d[N][N])
{
int i, j;
double temp, max_dif;
max_dif = 0.0;
for(i=0; i<N; i++){
    for(j=0; j<N; j++){
        temp = fabs(c[i][j] - d[i][j]); // 絶対値を取る
        if(temp > max_dif) max_dif = temp;
    }
}
return max_dif;
}