みなさん、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;
}