みなさん、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では、このように書きます。
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;
}