Help us understand the problem. What is going on with this article?

[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;
}
t--k
arara
「アイディアとテクノロジーでみんながハッピーになる社会を創る。」 アララ株式会社のエンジニアチームです。
https://www.arara.com
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away