2022/05/13(金)Ubuntu 22.04のBLAS (dgemm) をベンチマーク
BLAS, LAPACKとは何か、個々のBLASの特徴などは、20.04のときの記事を見て下さい。前回と同様、Reference BLAS, ATLAS, OpenBLAS, BLIS, MKLの5つのパッケージをaptでインストールし、dgemmを呼び出して計算時間を計測するプログラムを使って、update-alternativesの機能でlibblas.soを切り替えながらベンチマークを取り、計算時間を比較します。前回と大雑把な傾向はほとんど変わっていないのですが、MKLのパッケージが怪しげな挙動を示した(ような気がした)ので備忘録として書いています。
各BLASのインストール方法
- Reference BLAS
apt install libblas-dev
- ATLAS
apt install libatlas-base-dev
- OpenBLAS
apt install libopenblas-base apt install libopenblas-dev
- BLIS
apt install libblis-dev
- MKL
apt install intel-mkl最後のMKLですが、aptの実行中に
Use libmkl_rt.so as the default alternative to BLAS/LAPACK? [yes/no]と聞かれます。これには「no」と答えました。「yes」と答えると、「-lblas」を付けてコンパイルしたものはlibblas.soがリンクされずにMKL(libmkl_rt.so)がリンクされてしまってupdate-alternativesでlibblas.soの実体を切り替える作戦が不可能になってしまいました。何か自分が勘違いしてるのかも知れません。
これらのライブラリは、
update-alternatives --config libblas.so-x86_64-linux-gnuとすると
There are 5 choices for the alternative libblas.so-x86_64-linux-gnu (providing /usr/lib/x86_64-linux-gnu/libblas.so). Selection Path Priority Status ------------------------------------------------------------ 0 /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so 100 auto mode 1 /usr/lib/x86_64-linux-gnu/atlas/libblas.so 35 manual mode 2 /usr/lib/x86_64-linux-gnu/blas/libblas.so 10 manual mode * 3 /usr/lib/x86_64-linux-gnu/blis-openmp/libblas.so 80 manual mode 4 /usr/lib/x86_64-linux-gnu/libmkl_rt.so 1 manual mode 5 /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so 100 manual mode Press <enter> to keep the current choice[*], or type selection number:と表示され、libblas.soの実体をsymbolic linkで切り替えることができます。libblas.so.3の方も同様に
update-alternatives --config libblas.so.3-x86_64-linux-gnu
There are 5 choices for the alternative libblas.so.3-x86_64-linux-gnu (providing /usr/lib/x86_64-linux-gnu/libblas.so.3). Selection Path Priority Status ------------------------------------------------------------ 0 /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 100 auto mode 1 /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3 35 manual mode 2 /usr/lib/x86_64-linux-gnu/blas/libblas.so.3 10 manual mode * 3 /usr/lib/x86_64-linux-gnu/blis-openmp/libblas.so.3 80 manual mode 4 /usr/lib/x86_64-linux-gnu/libmkl_rt.so 1 manual mode 5 /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 100 manual mode Press <enter> to keep the current choice[*], or type selection number:と切り替えることができ、前回と同様、両方を同じものに切り替えて実験を行いました。
ベンチマーク
dgemmの実行時間計測に使ったプログラムは、以下のようなものです。前回の20.04のときと同じです。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
// prototype declaration
void dgemm_(char *transA, char *transB, int *m, int *n, int *k, double *alpha, double *A, int *ldA, double *B, int *ldB, double *beta, double *C, int *ldC);
int main(int argc, char **argv)
{
	int i, j;
	int m, n, k;
	int size;
	double *a, *b, *c;
	double alpha, beta;
	int lda, ldb, ldc;
	struct timespec ts1, ts2;
	size = atoi(argv[1]);
	m = size;
	n = size;
	k = size;
	a = (double *)malloc(sizeof(double) * m * k); // m x k matrix
	b = (double *)malloc(sizeof(double) * k * n); // k x n matrix
	c = (double *)malloc(sizeof(double) * m * n); // m x n matrix
	for (i=0; i<m; i++) {
		for (j=0; j<k; j++) {
			a[i + m * j] = rand() / (1.0 + RAND_MAX);
		}
	}
	for (i=0; i<k; i++) {
		for (j=0; j<n; j++) {
			b[i + k * j] = rand() / (1.0 + RAND_MAX);
		}
	}
	for (i=0; i<m; i++) {
		for (j=0; j<n; j++) {
			c[i + m * j] = 0;
		}
	}
	alpha = 1.;
	beta = 0.;
	lda = m; 
	ldb = k; 
	ldc = m; 
	// dgemm_(TransA, TransB, M, N, K, alpha, A, LDA, B, LDB, beta, C, LDC)
	// C = alpha * A * B + beta * C
	// A=M*K, B=K*N, N=M*N
	// Trans: "N"/"T"/"C"
	// LDA = number of row of A
	clock_gettime(CLOCK_REALTIME, &ts1);
	dgemm_("N", "N", &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
	clock_gettime(CLOCK_REALTIME, &ts2);
	printf("%g\n", (ts2.tv_sec - ts1.tv_sec) + (ts2.tv_nsec - ts1.tv_nsec) / 1e9);
	free(a);
	free(b);
	free(c);
	return 0;
}
単なる比較のための、単純な三重forループによる行列積のプログラムは、以下の通り。これも前回と同じ。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
double **alloc_matrix(int n, int m)
{
	double **a;
	int i;
	a = (double **)malloc(sizeof(double *) * n);
	a[0] = (double *)malloc(sizeof(double) * n * m);
	for (i=1; i<n; i++) {
		a[i] = a[0] + i * m;
	}
	return a;
}
void free_matrix(double **a)
{
	free(a[0]);
	free(a);
}
// a: n x m, b: m x s, c: n x s
void m_m_mul(double **a, double **b, double **c, int n, int m, int s)
{
	int i, j, k;
	int i1, j1, k1;
	for (i=0; i<n; i++) {
		for (j=0; j<s; j++) {
			c[i][j] = 0.0;
		}
	}
	for (i=0; i<n; i++) {
		for (j=0; j<s; j++) {
			for (k=0; k<m; k++) {
				c[i][j] += a[i][k] * b[k][j];
			}
		}
	}
}
int main(int argc, char **argv)
{
	int i, j, n;
	double **a;
	double **b;
	double **c;
	struct timespec ts1, ts2;
	n = atoi(argv[1]);
	a = alloc_matrix(n,n);
	b = alloc_matrix(n,n);
	c = alloc_matrix(n,n);
	for (i=0; i<n; i++) {
		for (j=0; j<n; j++) {
			a[i][j] = rand()/(1.0 + RAND_MAX);
			b[i][j] = rand()/(1.0 + RAND_MAX);
		}
	}
	clock_gettime(CLOCK_REALTIME, &ts1);
	m_m_mul(a, b, c, n, n, n);
	clock_gettime(CLOCK_REALTIME, &ts2);
	printf("%g\n", (ts2.tv_sec - ts1.tv_sec) + (ts2.tv_nsec - ts1.tv_nsec) / 1e9);
	free_matrix(a);
	free_matrix(b);
	free_matrix(c);
	return 0;
}
これらは、cc -O3 dgemm.c -lblas cc -O3 forloop.cのようにコンパイルします。ただし、dgemm.cの方のコンパイルをするときに、update-alternativesの選択をMKL以外にする必要がありました。MKLの状態でコンパイルすると、libblas.soがリンクされずにlibmkl_rt.soがリンクされてしまって、update-alternativesによる事後切り替えが不可能なbinaryができてしまいました。これもまた自分の勘違いかも。
さて、ベンチマークの結果は以下のようになりました。使ったCPUはcore i7 1195G7 (ノートPC) で、前回と異なるので絶対値には意味がないかも。

前回と似たような結果ですね。最速はOpenBLAS, BLIS, MKL。小サイズでMKLが妙に遅いのも同じ。