最近BLAS(Basic Linear Algebra Subprograms)について調べたので簡単な使い方等をまとめます。

はじめに

BLAS (Basic Linear Algebra Subprograms) は行列やベクトルに関する基本的な演算の仕様を定めたもので、この仕様に従った実装がいくつか存在します。私の使っているIntelのCPUを搭載したラップトップPCで利用可能な代表的な実装としては、

Reference BLAS: リファレンス実装。遅い。
OpenBLAS: 高速なオープンソース実装。
ATLAS: 自動チューニングするオープンソース実装。
Intel MKL: Intelによる高速な実装。有料。

があるようです。今回は、OpenBLASとATLASを実際にインストールして使ってみたので、その結果を記録します。

今回の実験環境は以下のとおりです。

OS: Debian GNU Linux 7.8 64bit
CPU: Intel Core i7-3667U 2.00GHz
メモリ: 8GB

OpenBLASのインストール

シェルで以下のように入力するとインストールできます。

$ git clone https://github.com/xianyi/OpenBLAS.git
$ make
$ sudo make install

ATLASのインストール

ATLASはコンパイル時にチューニングを行うので、コンパイルに非常に時間がかかります。

まず、Source Forge のページからソースコードをダウンロードしてきます。

ダウンロードしてきたソースコードを展開して、以下のようにconfigure、make、make installします。

$ tar jxvf atlas3.10.2.tar.bz2
$ mkdir build
$ cd build
$ ../ATLAS/configure
$ make
$ sudo make install

行列積を計算するルーチンの呼び出し

BLASのルーチンのうち、行列積を計算するC言語インターフェースは、 cblas_dgemm関数で、以下のプロトタイプになっている。

void cblas_dgemm(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);

これらのうち、CBLAS_ORDERとCBLAS_TRANSPOSEは、

enum CBLAS_ORDER     {CblasRowMajor=101, CblasColMajor=102};
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};

となっている。

これらは、cblas.hというヘッダーファイルに書かれている。

cblas_dgemmは、

C = αAB + βC

を計算する。引数の意味だが、以下のようになっている。

Order: 行列のメモリ上での配置。C言語であればCblasRowMajorを指定しておけば良い。
TransA, TransB: 入力行列が転置されているかどうか。
A, B, C: 行列
alpha, beta: 上のαとβの値
M, N, K: 行列のサイズ。A(M×N), B(N×K), C(M×K)
lda, ldb, ldc: 行列の下の行までの距離。大きな行列の一部分だけを取り出すときに使う。

時間の計測

今回はルーチンの呼び出しにどのくらい時間がかかるのかを測定する。 Linuxで時間を測定する方法はたくさんあるらしいが、ここでは、clock_gettime関数を使う。

clock_gettime関数を簡単に使うために、以下の関数を用意する。

static struct timespec ts_start, ts_stop;

void start_timer(void)
{
  clock_gettime(CLOCK_MONOTONIC_RAW, &ts_start);
}

void stop_timer(void)
{
  clock_gettime(CLOCK_MONOTONIC_RAW, &ts_stop);
}

void print_time(void)
{
  double time =
    ((double)ts_stop.tv_sec + (double)ts_stop.tv_nsec * 1e-9) - 
    ((double)ts_start.tv_sec + (double)ts_start.tv_nsec * 1e-9);
  printf("time = %f\n", time);
}

行列の割り当てと初期化

行列の割り当てと初期化を行う以下のルーチンを用意した。

double *alloc_matrix(int height, int width)
{
    return malloc(sizeof(double) * height * width);
}

void free_matrix(double *matrix)
{
    free(matrix);
}

void set_random_matrix(double *matrix, int height, int width)
{
    int i;
    for (i = 0; i < height * width; i++)
        matrix[i] = (double)rand() / RAND_MAX;
}

void clear_matrix(double *matrix, int height, int width)
{
    memset(matrix, 0, sizeof(double) * height * width);
}

行列積を計算ルーチンを呼び出すプログラム

上で用意した関数を使って、以下のプログラムを記述した。

#define N 2000

int main(void)
{
    double *a, *b, *c;
    
    a = alloc_matrix(N, N);
    b = alloc_matrix(N, N);
    c = alloc_matrix(N, N);
    
    set_random_matrix(a, N, N);
    set_random_matrix(b, N, N);
    clear_matrix(c, N, N);
    
    start_timer();
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
        N, N, N, 1.0, a, N, b, N, 0.0, c, N);
    stop_timer();
    print_time();

    return 0;
}

OpenBLASでのコンパイルと実行

上の手順でOpenBLASをインストールすると、OpenBLASは/opt/OpenBLAS以下に保存される。

実際にコンパイルするには以下のコマンドを入力する。ここではCのコンパイラとして、 Clang (Version 3.6) を利用し、ライブラリはスタティックリンクした。

$ clang -static -o dgemm_test -I/opt/OpenBLAS/include -L/opt/OpenBLAS/lib \
  dgemm_test.c -lopenblas -lpthread -lrt

OpenBLASでの実行結果

以下は上のプログラムを5回実行した結果である。

回数	秒
1回目	1.603888
2回目	1.588852
3回目	1.574097
4回目	1.627662
5回目	1.610352 (中央値)

ATLASでのコンパイルと実行

上の手順でATLASをインストールすると、ATLASは/usr/local/atlas以下に保存される。

$ clang -static -o dgemm_test -I/usr/local/atlas/include -L/usr/local/atlas/lib \
  dgemm_test.c -lcblas -latlas -lrt

ATLASでの実行結果