Edited at

[VC++] OpenBLASを使ってみた

More than 1 year has passed since last update.

みなさん、BLASって知ってますか? はい、そうですね。大きな行列演算をするときに使うとめっちゃ速い、あのBLASですね。

けど、実際どれくらい速くなるの?

よくわからないので、試しにOpenBLASで行列の掛け算をやってみました。


Visual StudioでOpenBLASを使う

Visual StudioでOpenBLASを使うには、だいたい3つくらい方法があります。


1. 自分でビルドする

ここにやり方が載ってます。

https://github.com/xianyi/OpenBLAS/wiki/How-to-use-OpenBLAS-in-Microsoft-Visual-Studio

正直、やりたくないです。


2. SourceForgeにあるコンパイル済みバイナリを落とす

ReadMeに書いてありますが、Windows用バイナリが公式に配布されています。

https://sourceforge.net/projects/openblas/files/

この記事を書いた時点では、最新のv0.2.20にはWindows用バイナリが見当たらず、一つ古いv0.2.19だとバイナリがダウンロードできます。


3. NuGetから取ってくる

公式なものじゃないですが、実はありました。

https://www.nuget.org/packages/OpenBLAS/

一番簡単ですが、バージョンが古く、更新されていないようです。


バイナリをダウンロードしてきた

今回は2の方法をとることにしました。SourceForgeにあるv0.2.19を使用します。

インクルードディレクトリ、ライブラリディレクトリの追加を忘れないでください。

リンカーのプロパティにある、追加の依存ファイルの設定も必要です。 libopenblas.dll.a を追加してください。似た名前の libopenblas.a もありますが、こちらはGCC用なので、Visual Studioでコンパイルするときは使いません。

また、libopenblas.dllはexeの作られるフォルダにコピーする(か、パスを通しておく)必要があります。


行列積を求める

とりあえず、乱数で行列を作ってみます。


行列を作ろう

    const int m = 1000, k = 2000, n = 3000;

std::vector<double> a(m * k);
std::vector<double> b(k * n);
std::vector<double> c_blas(m * n);
std::vector<double> c_loop(m * n);

std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<> dist(0.0, 10.0);

for (double& d : a) {
d = dist(mt);
}
for (double& d : b) {
d = dist(mt);
}


1000x2000の行列と、2000x3000の行列を掛け算しよう、という魂胆です。

ざっくり、2000回の掛け算を、1000x3000回しないといけないので、60億回の掛け算ですね。

今使っているパソコンは、4コアのようなので、4コアに分散させても、15億回の掛け算ですね。

ちなみに、CPUの周波数は2.6GHz、つまり、秒間26億クロックです。

普通に、ループで行列積を求めると、こうなります。

ちなみに、設定でOpenMPを有効にしてマルチスレッド化しています。

今回は深く触れませんが、for (int i = 0; i < m; i++)のループを並列演算しています。


ループで行列積を求める

#pragma omp parallel for

for (int i = 0; i < m; i++) {
for (int t = 0; t < k; t++) {
for (int j = 0; j < n; j++) {
c_loop[i * n + j] += a[i * k + t] * b[t * n + j];
}
}
}

続いて、BLASでは、このように書きます。


gemmを使う

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 0.0, a.data(), k, b.data(), n, 1.0, c_blas.data(), n);


gemmというのが行列積を求める関数で、dgemmはdouble型のgemmです。

もう少し厳密には、単なる行列積よりは少し高機能で、

C = \alpha AB + \beta C

が計算できます。

引数が多くて嫌になっちゃいますが、各引数の説明は、以下の通りです。


  • CblasRowMajor


    • データの並びが横優先(行を並べる形)になっていることを示します。C/C++ではこう扱うことが多いです。



  • CblasNoTrans (2回)


    • それぞれA, Bの行列を、転置せずにそのまま使うことを表します。



  • m, n, k


    • Aの行列はm行k列、Bの行列はk行n列、Cの行列はm行n列であることを表します。



  • 1.0 (α)


    • 上式のαの値です。単に行列積を求めるときは、1.0を指定します。



  • a.data()



    • double*型で配列Aを渡します。



  • k


    • leading dimensionと呼ばれる数で、通常は、横優先なら列数を、縦優先なら行数を指定します。



  • b.data(), n


    • 上記a.data(), kと同様です。



  • 0.0 (β)


    • 上式のβの値です。今回、Cはあらかじめ0に初期化されているので何でもいいのですが、0.0を指定しました。



  • c.data(), n


    • 上記a.data(), kと同様です。



また、計算結果の確認にも、BLASを使います。


計算結果確認

    // そのために、まず、c_blasからc_loopを引き算する。

cblas_daxpy(c_blas.size(), -1.0, c_loop.data(), 1, c_blas.data(), 1);
// 続いて、引き算した各要素の絶対値の和をとる。
auto diff = cblas_dasum(c_blas.size(), c_blas.data(), 1);
std::cout << "Diff: " << diff << std::endl;

axpyは、

Y = \alpha X + Y

を計算する関数です。αを-1にすることで、BLASで計算した結果から、ループで計算した結果を引き算しています。

また、asumは、各要素の絶対値の和を返す関数です。

BLASとループで同じ計算結果ならば、asumした結果はほぼゼロになるので、ちゃんと計算ができていることを確認できます。


動かしてみた

Releaseビルドで、最適化かけてやってます。

CPUはCore i7 6500Uで、メモリは16GBです。

> .\test_openblas.exe

Loop:
elapsed time: 2.327 sec.
BLAS:
elapsed time: 0.168 sec.
Diff: 0.000149777

え、うそやん・・・・ なにこれ。めっちゃ速い。


結論

OpenBLASめっちゃ速い。


コード全体

一応貼っておきます。


コード全体

#include <vector>

#include <iostream>
#include <random>
#include <chrono>
#include <cblas.h>

int main(void)
{
const int m = 1000, k = 2000, n = 3000;

std::vector<double> a(m * k);
std::vector<double> b(k * n);
std::vector<double> c_blas(m * n);
std::vector<double> c_loop(m * n);

std::random_device rd;
std::mt19937 mt(rd());
std::uniform_real_distribution<> dist(0.0, 10.0);

for (double& d : a) {
d = dist(mt);
}
for (double& d : b) {
d = dist(mt);
}

// 普通のループで行列の積を求める
std::cout << "Loop:" << std::endl;
auto begin = std::chrono::high_resolution_clock::now();

#pragma omp parallel for
for (int i = 0; i < m; i++) {
for (int t = 0; t < k; t++) {
for (int j = 0; j < n; j++) {
c_loop[i * n + j] += a[i * k + t] * b[t * n + j];
}
}
}

auto end = std::chrono::high_resolution_clock::now();
auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << " elapsed time: " << elapsed / 1000.0 << " sec." << std::endl;

// OpenBLAS で行列の積を求める
std::cout << "BLAS:" << std::endl;
begin = std::chrono::high_resolution_clock::now();

cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0, a.data(), k, b.data(), n, 0.0, c_blas.data(), n);

end = std::chrono::high_resolution_clock::now();
elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - begin).count();
std::cout << " elapsed time: " << elapsed / 1000.0 << " sec." << std::endl;

// 一応、求めた答えがだいたい同じか確かめる。
// そのために、まず、c_blasからc_loopを引き算する。
cblas_daxpy(c_blas.size(), -1.0, c_loop.data(), 1, c_blas.data(), 1);
// 続いて、引き算した各要素の絶対値の和をとる。
auto diff = cblas_dasum(c_blas.size(), c_blas.data(), 1);
std::cout << "Diff: " << diff << std::endl;
}