テストプログラムの作成｜行列積高速化#3
この記事は、以下の記事を分割したものです。
［元の記事］行列積計算を分割してみる
一括で読みたい場合は、元の記事をご覧ください。
プログラムの高速化をしていると、間違った変更をしてしまいがちです。また、計算速度がどのくらいになっているのかを随時確認する必要もあります。そこで、以下の2つのテストプログラムを作成しました。
（１）計算結果を確認するプログラム：unit_test
（２）計算速度を確認するプログラム：speed_test, total_test
計算速度を確認するプログラムが2つありますが、speed_testは簡易的に確認するため、total_testはグラフ描画用に細かく測定するために作成しています。プログラム自体は、ループの刻み幅が粗いか細かいかを除いて、ほとんど同じです。
テストプログラムの実装を考える前に、いろいろ準備する必要があります。
3-1. 正解値の準備
計算結果を確認するには、正解値が必要です。今回は、正解としてNetlib.orgで公開されているBLASライブラリを使用します。もし、システムにBLASライブラリがパッケージとしてインストールされている場合は、それを使用しても構いません。
BLASは、Fortran言語用のライブラリです。しかし、今回はC言語インターフェースを使用したいので、上記ページからダウンロードできるcblas.tgzもコンパイルして使用します。
3-2. インターフェースの準備
行列積DGEMMは引数が多く、これを毎回記述するのは大変なので、次のような構造体を用意しておきます。
typedef void (*dgemm_func_t)
(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
const int K, const double alpha, const double *A,
const int lda, const double *B, const int ldb,
const double beta, double *C, const int ldc);
typedef struct _dgemm_test_t {
dgemm_func_t dgemm;
enum CBLAS_ORDER Order;
enum CBLAS_TRANSPOSE TransA;
enum CBLAS_TRANSPOSE TransB;
int M;
int N;
int K;
double alpha;
double* A;
int lda;
double* B;
int ldb;
double beta;
double* C;
int ldc;
} dgemm_test_t;dgemm_func_tは、CBLASのDGEMM関数の関数ポインタ型です。これは、cblas.hに記載されているCblas_dgemm関数を参考に作成しています。なぜ、関数ポインタ型を定義しているかというと、正解の関数と自作の関数を、テストプログラムで差し替えられるようにするためです。
dgemm_test_tは、テストプログラム内で利用するDGEMM関数と引数をセットにした構造体です。引数をそのままにして関数だけを差し替えるといった使い方ができるようにしています。
また、dgemm_test_t型を引数とした行列計算を簡単に実行できるように、次のような関数を用意しました。
void do_dgemm( dgemm_test_t* test ){
test->dgemm(test->Order, test->TransA, test->TransB, test->M, test->N, test->K,
test->alpha, test->A, test->lda, test->B, test->ldb,
test->beta , test->C, test->ldc);
}3-3. 行列初期化関数の準備
BLASルーチンでは、行列Aは、メモリ上の行列サイズがLDA×N、実際に使用する行列サイズをM×Nとして取り扱います。このLDA（及びLDB, LDC）を、Leading Dimensionと呼びます。これにより、連続一次元的に確保されたメモリ領域を、切り取って使うことができるようになっています。
Fortranの二次元配列A(M,N)は、１次元目Mがメモリ上は連続（Column Major）になっているので、行をLeading DimensionとしてA(LDA,N)という形にメモリ領域を確保します。しかし、C言語の二次元配列A[M][N]は、２次元目Nがメモリ上は連続（Row Major）になります。
Fortran用BLASでは、Leading Dimensionは自動的に１次元目を指すことになっています。しかし、C言語用BLASでは、１次元目か２次元目かを引数OrderでLeading Dimensionを選択できるようになっています。C言語用BLASでは、OrderがCblasColMajorの場合は１次元目（Fortranと同じ）を、CblasRowMajorの場合は２次元目を選んだことになります。
ただ、この仕組みがあるために、Leading Dimensionを考慮した初期化関数が必要になります。
計算速度の確認では、計算結果を問わないので、固定値で埋め尽くす初期化でも構いません。しかし、計算結果の確認は、精度比較で確認するので、倍精度の有効桁数の数字が入ったランダムな値が望ましいです。
そこで、下記の２種類の初期化関数を用意しました。
void init_matrix( const int m, const int n,
double* A, const int di, const int dj, double value )
{
double *a = A;
for( int j=0; j<n; j++ ){
for( int i=0; i<m; i++ ){
*a = value;
a = a + di;
}
a = a - m*di + dj;
}
}
void rand_matrix( const int m, const int n,
double* A, const int di, const int dj, uint64_t seed )
{
unsigned short xseed[3]={0};
if( seed != 0 ){
xseed[0] = (seed )&0xffff;
xseed[1] = (seed>>16)&0xffff;
xseed[2] = (seed>>32)&0xffff;
seed48(xseed);
}
double *a = A;
for( int j=0; j<n; j++ ){
for( int i=0; i<m; i++ ){
*a = drand48();
a = a + di;
}
a = a - m*di + dj;
}
}init_matrix関数は、1つの値で埋め尽くす初期化関数です。一方、rand_matrix関数は、倍精度の乱数で埋め尽くす初期化関数です。
行列Aの１次元目と２次元目のストライドサイズdi,djを指定できるようにすることで、CblasRowMajorとCblasColMajorの両方に対応させています。
ストライドってなあに？という方は、下記の記事をご覧ください。
3-4. エラーチェック関数の準備
2つの浮動小数の値の一致は、行列の各要素の正解値と計算値の相対誤差が、倍精度マシンイプシロンより小さいことで等しいと判定します。
相対誤差 ＝ |正解値ー計算値|÷|正解値| ≦ マシンイプシロン
実際には、|正解値|を右辺にうつして、下記の式で判定します。
|正解値ー計算値| ≦ マシンイプシロン×|正解値|
しかし、行列Cの各要素の計算では、行列要素abをK個足し上げた総和計算を行います。一般に、総和計算は、総和サイズKに比例した丸め誤差が発生します。特に、足し上げ順序を入れ替えると、丸め誤差が顕在化します。
実際には、丸め誤差の方向（多いのか少ないのか）が要素によって異なるため、√Kにおおよそ比例すると言われています。そこで、判定条件は以下にします。
|正解値ー計算値| ≦ マシンイプシロン×|正解値|×sqrt(K)
上式が成立しない場合が、エラーになります。｜正解値｜は、入力された値のどちらが正解値か分からないので、コード上は２つの値の絶対値の大きい方で代用することにします。
以上を踏まえて、行列の全要素を比較する関数は次のようになります。
static int check_matrix( const int m, const int n, const int k,
const double *C1, const int di1, const int dj1,
const double *C2, const int di2, const int dj2 )
{
for( int j=0; j<n; j++ ){
for( int i = 0; i<m; i++ ){
double c1 = fabs(*C1);
double c2 = fabs(*C2);
double cmax = ( c1>c2 ? c1 : c2 );
double cdiff = fabs(c1-c2);
double sqrtk = sqrt((double)k);
double eps = DBL_EPSILON * sqrtk * cmax;
if( cdiff > eps ){
fprintf(stderr,"[ERROR] An element C(%d,%d) is invalid : t1=%G, t2=%G, diff=%G, eps=%G\n",i,j,*C1,*C2,cdiff,eps);
return 1;
}
C1 = C1 + di1;
C2 = C2 + di2;
}
C1 = C1 - m*di1 + dj1;
C2 = C2 - m*di2 + dj2;
}
return 0;
}また、確保したメモリ領域を全て比較する必要はなく、比較範囲は設定によって変わるので、次のようなチェック関数を用意しておきます。
int check_error( const dgemm_test_t* t1, const dgemm_test_t* t2 ){
double* C1 = t1->C;
double* C2 = t2->C;
int ldc1 = t1->ldc;
int ldc2 = t2->ldc;
if( t1->M != t2->M ){
fprintf(stderr,"[ERROR] M is not equal : t1.M=%d, t2.M=%d \n",t1->M,t2->M);
}
if( t1->N != t2->N ){
fprintf(stderr,"[ERROR] N is not equal : t1.N=%d, t2.N=%d \n",t1->N,t2->N);
}
int m = t1->M;
int n = t1->N;
int error = 0;
if( t1->Order == CblasRowMajor && t2->Order == CblasRowMajor ){
error = check_matrix( m, n, C1, 1, ldc1, C2, 1, ldc2 );
}else if( t1->Order == CblasColMajor && t2->Order == CblasRowMajor ){
error = check_matrix( m, n, C1, ldc1, 1, C2, 1, ldc2 );
}else if( t1->Order == CblasRowMajor && t2->Order == CblasColMajor ){
error = check_matrix( m, n, C1, 1, ldc1, C2, ldc2, 1 );
}else if( t1->Order == CblasColMajor && t2->Order == CblasColMajor ){
error = check_matrix( m, n, C1, ldc1, 1, C2, ldc2, 1 );
}
return error;
}3-5. スピードチェック関数の準備
計算速度の確認は、DGEMMの呼び出し前後の経過時間を測定します。経過時間の測定については、下記の記事をご覧ください。
スピードチェック関数の実装は、下記のようになります。
double check_speed( dgemm_test_t* test ){
if( test->Order == CblasRowMajor ){
init_matrix( test->M, test->K, test->A, 1, test->lda, 1e0 );
init_matrix( test->K, test->N, test->B, 1, test->ldb, 1e0 );
init_matrix( test->M, test->N, test->C, 1, test->ldc, 0e0 );
}else{
init_matrix( test->M, test->K, test->A, test->lda, 1, 1e0 );
init_matrix( test->K, test->N, test->B, test->ldb, 1, 1e0 );
init_matrix( test->M, test->N, test->C, test->ldc, 1, 0e0 );
}
double t1,t2;
int error;
t1 = get_realtime();
do_dgemm( test );
t2 = get_realtime();
return (t2-t1);
}計算値はなんでも良いので、スピードチェック関数の中で行列を初期化してしまします。
3-6. 演算数の見積もり
BLASの行列積計算ルーチンDGEMMは、次のような公式を計算しています。
C = alpha * A * B + beta * C
これの最も簡単なプログラムは下記のようになります。
for( j=0; j<N; j++){
for( i=0; i<M; i++){
C[i][j] = beta*C[i][j];
for( k=0; k<K; k++){
C[i][j] = alpha*A[i][k]*B[k][j] + C[i][j];
}}}betaの乗算は、i,jのループに依存し、乗算が１つあるのでM*N回の乗算が行われます。alphaと行列A、Bの行列積は、i,j,kのループに依存し、乗算が２つあるので2*M*N*K回の乗算が行われます。また、加算は、１つのi,jの組みに対してK回行われるので、M*N*K回の加算が行われます。これを合計すると、次のようになります。
2*M*N*K + M*N + M*N*K = M*N*(3*K+1)
実際には、alphaの計算をkループの外側に出すことも可能です。
for( j=0; j<N; j++){
for( i=0; i<M; i++){
ab = 0;
for( k=0; k<K; k++){
ab = ab + A[i][k]*B[k][j];
}
C[i][j] = alpha*ab + beta*C[i][j];
}}この場合だと、最深ループの行列Aと行列Bの行列積は、加算１つと乗算１つで構成されています。これが、M*N*K回計算されるので、演算数は2*M*N*K個になります。また、最新ループの外側の行列Cに書き込む計算は、加算１つと乗算２つで構成されています。これが、M*N回計算されるので、演算数は3*M*N個となります。したがって、演算数の総数は、下記のようになります。
2*M*N*K + 3*M*N = M*N*(2*K+3)
一方、行列積の本質的な計算は、行列Aと行列Bの行列積部分なので、演算数として2*M*N*K個が使われることもあります。まとめると、次のようになります。
（１）2*M*N*K : プログラム間の性能比較でよく用いられる近似的な演算数
（２）M*N*(2*K+3) : 最小化された正確な演算数
（３）M*N*(3*K+1)：alpha乗算を冗長化した正確な演算数
今回は、他のライブラリの実装との性能比較ではなく、理論ピーク性能に対する比率を計算するので、正確な演算数が必要です。最適化を進めていくと、alpha乗算は最小ではなくなりますが、（２）の演算数M*N*(2*K+3)を用いることにします。
3-7. 計算結果を確認するプログラム
以上で準備ができたので、計算結果を確認するテストプログラムを作成します。行列サイズは、今後の場合分けに備えて1023×1023に設定します。今後、2^n毎に場合わけをする可能性があるので、全ての場合を通るように2^n-1のサイズにしています。
また、Order, TransA, TransBの全ての可能な設定の組み合わせについて確認するため、全部で８パターンをテストする必要があります。
テストプログラムを下記の通りになりました。
#define SIZE 1023
#define SEED 23091
int main(int argc, char** arv){
enum CBLAS_ORDER majors[2]={CblasRowMajor,CblasColMajor};
enum CBLAS_TRANSPOSE transposes[2]={CblasNoTrans,CblasTrans};
const char* cmajor[2]={"Row-major","Col-major"};
const char* ctrans[2]={"NoTrans","Trans"};
double *A, *B, *C, *D;
dgemm_test_t cblas = { cblas_dgemm, CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE, SIZE, SIZE, 1e0, NULL, SIZE, NULL, SIZE, 1e0, NULL, SIZE };
dgemm_test_t myblas= { cblas_dgemm, CblasRowMajor, CblasNoTrans, CblasNoTrans, SIZE, SIZE, SIZE, 1e0, NULL, SIZE, NULL, SIZE, 1e0, NULL, SIZE };
A = calloc( SIZE*SIZE, sizeof(double) );
B = calloc( SIZE*SIZE, sizeof(double) );
C = calloc( SIZE*SIZE, sizeof(double) );
D = calloc( SIZE*SIZE, sizeof(double) );
rand_matrix( SIZE, SIZE, A, 1, SIZE, SEED ); // RowMajor
rand_matrix( SIZE, SIZE, B, 1, SIZE, SEED ); // RowMajor
cblas.A = A;
cblas.B = B;
cblas.C = C;
myblas.A = A;
myblas.B = B;
myblas.C = D;
int error = 0;
for( int iorder=0; iorder<2; iorder++ ){
for( int itransa=0; itransa<2; itransa++ ){
for( int itransb=0; itransb<2; itransb++ ){
cblas.Order = majors[iorder];
cblas.TransA = transposes[itransa];
cblas.TransB = transposes[itransb];
myblas.Order = majors[iorder];
myblas.TransA = transposes[itransa];
myblas.TransB = transposes[itransb];
printf("Case : Order=%s, TransA=%s, TransB=%s ... ",cmajor[iorder],ctrans[itransa],ctrans[itransb]);
init_matrix( SIZE, SIZE, cblas.C , 1, SIZE, 0e0 );
init_matrix( SIZE, SIZE, myblas.C, 1, SIZE, 0e0 );
do_dgemm( &cblas );
do_dgemm( &myblas );
error = check_error( &cblas, &myblas );
if( error ){
printf("NG\n");
break ;
}else{
printf("OK\n");
}
}
}
}
free(A);
free(B);
free(C);
free(D);
return error;
}dgemm_test_t構造体のmyblas変数の第一変数dgemmは、後ほど自作プログラムの関数ポインタに置き換えます。
上記のテストプログラムを実行し、問題がない場合は、下記のように表示されます。
Case : Order=Row-major, TransA=NoTrans, TransB=NoTrans ... OK
Case : Order=Row-major, TransA=NoTrans, TransB=Trans ... OK
Case : Order=Row-major, TransA=Trans, TransB=NoTrans ... OK
Case : Order=Row-major, TransA=Trans, TransB=Trans ... OK
Case : Order=Col-major, TransA=NoTrans, TransB=NoTrans ... OK
Case : Order=Col-major, TransA=NoTrans, TransB=Trans ... OK
Case : Order=Col-major, TransA=Trans, TransB=NoTrans ... OK
Case : Order=Col-major, TransA=Trans, TransB=Trans ... OK3-8. 計算速度を確認するプログラム
計算性能は行列サイズによって変わってくるので、粗くサイズを変えて計算性能を確認するようにしました。テストプログラムは、下記のようになりました。
#define MAX_SIZE 2048
int main(int argc, char** arv){
int error = 0;
flops_info_t cpu;
double *A, *B, *C;
dgemm_test_t myblas = { cblas_dgemm, CblasRowMajor, CblasNoTrans, CblasNoTrans, MAX_SIZE, MAX_SIZE, MAX_SIZE, 1e0, NULL, MAX_SIZE, NULL, MAX_SIZE, 1e0, NULL, MAX_SIZE };
A = calloc( MAX_SIZE*MAX_SIZE, sizeof(double) );
B = calloc( MAX_SIZE*MAX_SIZE, sizeof(double) );
C = calloc( MAX_SIZE*MAX_SIZE, sizeof(double) );
myblas.A = A;
myblas.B = B;
myblas.C = C;
peak_flops( &cpu );
double mpeak = cpu.mflops_double_max / cpu.num_cores;
double bpeak = cpu.mflops_double_base / cpu.num_cores;
printf("Max Peak MFlops per Core: %G MFlops \n",mpeak);
printf("Base Peak MFlops per Core: %G MFlops \n",bpeak);
printf("size , elapsed time[s], MFlops, base ratio[%%], max ratio[%%] \n");
for( size_t n=16; n<=MAX_SIZE; n=n*2 ){
myblas.M = n;
myblas.N = n;
myblas.K = n;
myblas.lda = n;
myblas.ldb = n;
myblas.ldc = n;
double nflop = (2*((double)myblas.K)+3) * ((double)myblas.M) * ((double)myblas.N);
double dt = check_speed( &myblas );
double mflops = nflop / dt / 1000 / 1000;
printf("%6u, %15G, %15G, %15G, %15G \n",n,dt,mflops,mflops/bpeak*100,mflops/mpeak*100);
}
free(A);
free(B);
free(C);
return error;
}現在は、Netlib.orgのリファレンスBLASの測定をしています。
これを実行すると、下記のように表示されます。
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.86102E-06, 3131.75, 6.74945, 5.93134
32, 1.40667E-05, 4877.34, 10.5115, 9.23738
64, 0.000114918, 4669.22, 10.063, 8.84322
128, 0.000982046, 4321.04, 9.31258, 8.18378
256, 0.00650811, 5186, 11.1767, 9.82196
512, 0.052299, 5147.74, 11.0943, 9.74951
1024, 0.489076, 4397.33, 9.47701, 8.32828
2048, 4.04466, 4250.65, 9.16089, 8.05048 ここから、Netlib.orgからダウンロードできるBLASは、基本周波数の場合で理論ピーク性能比9~11%が出ていることがわかります。
次の記事
元の記事はこちらです。
ソースコードはGitHubで公開しています。
この記事が気に入ったらサポートをしてみませんか？