VS2012用 OpenBLAS + LAPACK をFORTRAN無しでビルドする方法

概要

SIMD を用いて非常に効率よく計算できる線形代数用のライブラリ OpenBLAS を Visual Studio で使うため、ビルドする方法を調べました。 ビルドにあたり、下記の方針をとりました。

• Visual Studio 2012 はもう古いですが、辛うじてビルドできたので、今回も 2012 をターゲットにします。 2017 等にバージョンを上げる分には読み替えで対応できると思います。
• OpenBLAS のビルドには、Visual Studio 2012 の他に Clang、 cmake、 Ninja を使います。
• OpenBLAS の中に LAPACK も組み込まれているのですが、ビルドするためには FORTRAN が必要となります。Windows版 Flang が見つからず、OpenBLAS の CMakeに組み込まれている手順では LAPACKがビルドできないため、別途 CLapack をダウンロードしてビルドします。LAPACKE は OpenBLAS に含まれているものを使用します。
• Windows版 OpenBLAS についてはサポート外のようですが、スタティックライブラリの作成をトライします。
• 各ライブラリのバージョンは、 OpenBLAS : 0.3.19、 CLAPACK : 3.2.1 を使います。

ソース一式

• OpenBLAS + CLAPACK ソース

Visual Studio 2012 でビルドできるように、ファイルを追加したり、一部ソースを書き換えています。変更内容は"元ソースからの変更内容"項に記載します。

上記ソースからのビルド手順

• Clang、 cmake、 Ninja の Windows用バイナリをそれぞれ入手する。
• ダウンロードした上記 "OpenBLAS + CLAPACK ソース" zip ファイルの内容をビルド作業用のフォルダに展開する。
• Visual Studio 2012 x64 Cross Tools コマンドプロンプトを開き、 Clang、cmake、 ninja へのパスを通す。
• コマンドプロンプトのカレントディレクトリを <展開したフォルダ>\OpenBLAS-0.3.19\build_vs2012 に移動する。
• cmake_args.bat をメモ帳等で開き、適宜設定を書き換える。スタティックライブラリ or dll の切り替え、DYNAMIC_ARCH (OpenBLASの複数コアタイプ自動切換え機能) ON/OFF の切り替え、対象コアタイプの設定、Cランタイムもスタティックにするかどうか、等を切り替えできます。(なお、DYNAMIC_ARCH=ONのときは組み込みたいコアタイプを DYNAMIC_LIST にセミコロン区切りで指定することになります。コアタイプ名は getarch.c 等を見ればわかるようです)
• "cmake_args.bat" を実行する。→ ビルドの前処理が実行されます。
• "cmake --build ." を実行する。→ ビルドの本処理が実行されます。ビルド成功すると、build_vs2012/lib/Release/openblas.lib が生成されます。
• "create_include.bat" を実行する。 → ライブラリを使う際に必要となるインクルードファイルを build_vs2012/include フォルダ下にまとめます。
• clapack-3.2.1-CMAKE\build_vs2012\lapack.sln を Visual Studio 2012 で開き、Cランタイムもスタティックにするかどうかの設定を cmake_args.bat の設定に合わせ、 (-DMSVC_STATIC_CRT=ON なら /MT、 -DMSVC_STATIC_CRT=OFF なら /MD)、 ビルドする。→ ビルド成功すると、lib\Release_x64\lapack.lib が生成されます。この中に LAPACK、LAPACKE が入っています。

元ソースからの変更内容

• OpenBLAS-0.3.19\build_vs2012 フォルダの作成、cmake_args.bat 追加、create_include.bat 追加、f77blas.h 追加 (f77blas.h は作り方がわからなかったため、OpenBLAS バイナリ配布パッケージからコピーしてきました。)
• OpenBLAS-0.3.19\utest\utest_main2.cの <complex.h> をインクルードする行のコメントアウト (VS2012が C99に対応していない問題に対する処置。VS2013以降であれば不要？)
• OpenBLAS-0.3.19\utest\utest.hの <inttypes.h> をインクルードする行のコメントアウト、及びintmax_t、uintmax_t のtypedef追加 (VS2012が C99に対応していない問題に対する処置。VS2013以降であれば不要？)
• clapack-3.2.1-CMAKE\build_vs2012 フォルダの作成、及びソリューション一式の作成
  • 空プロジェクトを作成し、OpenBLAS-0.3.19\lapack-netlib\LAPACKE\src、OpenBLAS-0.3.19\lapack-netlib\LAPACKE\utils\のcソース一式、clapack-3.2.1-CMAKE\SRC\のCソース一式、clapack-3.2.1-CMAKE\INSTALL\ のメイン関数があるもの以外、かつ、second.c、dsecnd.c 以外のCソース一式、clapack-3.2.1-CMAKE\F2CLIBS\libf2c\のCソース一式をフィルタで分けて追加。
  • etime_.c、dtime_.c の extern "C" 関連の記述位置の修正 (恐らくバグ)
  • lapacke 下のソース一式、及び libf2c下の c_cos.c、c_exp.c、c_log.c、c_sin.c、c_sqrt.c を 「C++ としてビルド」の設定に変更 (VS2012が C99に対応していない問題に対する処置。VS2013以降であれば不要？)
• clapack-3.2.1-CMAKE\F2CLIBS\libf2c\arith.h (空ファイル)の追加

試しに使ってみた結果

処理速度や生成される実行ファイルのサイズが気になりますので、下記検証コードで試しに使ってみました。 画像処理ライブラリ MIST に行列計算クラスや、BLAS/Lapackラップ関数が入っていますので、それも使って検証しています。

検証用ソースコード

// 参考にしたサイト:
// * https://qiita.com/t--k/items/69c43a667a1283578012
// * https://auewe.hatenablog.com/entry/2013/11/25/024149

#include <windows.h>
#include <cstdlib>
#include <cstdio>
#include <cmath>

#define USE_CBLAS
#define USE_F77BLAS
#define USE_MIST
#define EXAMINE_LAPACK

#ifdef USE_CBLAS
#include "cblas.h"
#endif

#ifdef USE_F77BLAS
#include "f77blas.h"
#endif

#ifdef USE_MIST
#include "mist/numeric.h"
#endif

void simple_multiply(double *a, double *b, int m, int n, int k, double *c){
	for (int i3 = 0; i3 < m; ++i3) {
		for (int i2 = 0; i2 < n; ++i2) {
			for (int i1 = 0; i1 < k; ++i1) {
				c[i3 + i2 * m] += a[i3 + i1 * m] * b[i1 + i2 * k];
			}
		}
	}
}

#ifdef USE_CBLAS
void cblas_multiply(double *a, double *b, int m, int n, int k, double *c){
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1.0, a, m, b, k, 0.0, c, m);
}
#endif

#ifdef USE_F77BLAS
void f77blas_multiply(double *a, double *b, int m, int n, int k, double *c){
	double alpha = 1.0, beta = 0.0;
	char tra[2] = "N", trb[2] = "N";
	int lda = m, ldb = k, ldc = m;
	dgemm_(tra, trb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);
}
#endif

#ifdef USE_MIST
void mist_multiply(double *a, double *b, int m, int n, int k, double *c, LARGE_INTEGER &c1, LARGE_INTEGER &c2){
	mist::matrix<double> ma(m, k), mb(k, n);
	for (size_t i = 0; i < ma.size(); ++i) ma[i] = a[i];
	for (size_t i = 0; i < mb.size(); ++i) mb[i] = b[i];
	::QueryPerformanceCounter(&c1);
	mist::matrix<double> mc = ma * mb;
	::QueryPerformanceCounter(&c2);
	for (size_t i = 0; i < mc.size(); ++i) c[i] = mc[i];
}

void mist_numeric_multiply(double *a, double *b, int m, int n, int k, double *c, LARGE_INTEGER &c1, LARGE_INTEGER &c2){
	mist::matrix<double> ma(m, k), mb(k, n);
	for (size_t i = 0; i < ma.size(); ++i) ma[i] = a[i];
	for (size_t i = 0; i < mb.size(); ++i) mb[i] = b[i];
	mist::matrix<double> mc(m, n);
	::QueryPerformanceCounter(&c1);
	mist::multiply(ma, mb, mc, false, false, 1.0, 0.0);
	::QueryPerformanceCounter(&c2);
	for (size_t i = 0; i < mc.size(); ++i) c[i] = mc[i];
}

#endif

double diff(double *c1, double *c2, int count){
	double sum = 0;
	for (int i = 0; i < count; ++i){
		sum += std::abs(c1[i] - c2[i]);
	}
	return sum;
}

void dump(double *c, int m, int n){
	for (int j = 0; j < m; ++j){
		for (int i = 0; i < n; ++i){
			std::printf("%f ", c[j*n + i]);
		}
		std::printf("\n");
	}
}

int main(){
	// OpenBLAS の関数が内部で使用するスレッド数
#ifdef USE_CBLAS
	openblas_set_num_threads(1);
//	openblas_set_num_threads(4);
#endif

	// A x B = C
	// A : m 行 k 列
	// B : k 行 n 列
	// C : m 行 n 列
	const int m = 1500, k = 1000, n = 4000;

	double *a = static_cast<double *>(std::malloc(m * k * sizeof(double)));
	double *b = static_cast<double *>(std::malloc(k * n * sizeof(double)));
	double *c_simple = static_cast<double *>(std::malloc(m * n * sizeof(double)));
	double *c_cblas = static_cast<double *>(std::malloc(m * n * sizeof(double)));
	double *c_f77blas = static_cast<double *>(std::malloc(m * n * sizeof(double)));
	double *c_mist = static_cast<double *>(std::malloc(m * n * sizeof(double)));
	double *c_mist_numeric = static_cast<double *>(std::malloc(m * n * sizeof(double)));

	for (int i = 0; i < m * k; ++i) a[i] = static_cast<double>(std::rand()) * 40.0 / static_cast<double>(RAND_MAX - 1);
	for (int i = 0; i < k * n; ++i) b[i] = static_cast<double>(std::rand()) * 40.0 / static_cast<double>(RAND_MAX - 1);

	LARGE_INTEGER c[8], freq;

	::QueryPerformanceFrequency(&freq);

	::QueryPerformanceCounter(&c[0]);

	simple_multiply(a, b, m, n, k, c_simple);

	::QueryPerformanceCounter(&c[1]);

#ifdef USE_CBLAS
	cblas_multiply(a, b, m, n, k, c_cblas);
#endif

	::QueryPerformanceCounter(&c[2]);

#ifdef USE_F77BLAS
	f77blas_multiply(a, b, m, n, k, c_f77blas);
#endif

	::QueryPerformanceCounter(&c[3]);

#ifdef USE_MIST
	mist_multiply(a, b, m, n, k, c_mist, c[4], c[5]); // MIST独自実装。普通にC++で書いてある。
#endif

#ifdef USE_MIST
	mist_numeric_multiply(a, b, m, n, k, c_mist_numeric, c[6], c[7]); // MIST のOpenBLAS/LAPACKラップ機能。OpenBLAS の dgemm をf77のAPI経由で呼び出している。
#endif

	::QueryPerformanceCounter(&c[5]);


	std::printf("simple: %f msec\n", static_cast<double>(c[1].QuadPart - c[0].QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart));

#ifdef USE_CBLAS
	std::printf("cblas: %f msec\n", static_cast<double>(c[2].QuadPart - c[1].QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart));
#endif

#ifdef USE_F77BLAS
	std::printf("f77blas: %f msec\n", static_cast<double>(c[3].QuadPart - c[2].QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart));
#endif

#ifdef USE_MIST
	std::printf("mist: %f msec\n", static_cast<double>(c[5].QuadPart - c[4].QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart));
	std::printf("mist numeric: %f msec\n", static_cast<double>(c[7].QuadPart - c[6].QuadPart) * 1000.0 / static_cast<double>(freq.QuadPart));
#endif


#ifdef USE_CBLAS
	std::printf("diff (cblas - simple): %f\n", diff(c_cblas, c_simple, m * n));
#endif

#ifdef USE_F77BLAS
	std::printf("diff (f77blas - simple): %f\n", diff(c_f77blas, c_simple, m * n));
#endif

#ifdef USE_MIST
	std::printf("diff (mist - simple): %f\n", diff(c_mist, c_simple, m * n));
	std::printf("diff (mist_numeric - simple): %f\n", diff(c_mist_numeric, c_simple, m * n));
#endif

#ifdef EXAMINE_LAPACK
	mist::matrix<double> mat(2, 2);
	mat(0,0) = 2.0, mat(0,1) = 1.0;
	mat(1,0) = 3.0, mat(1,1) = 0.5;
	mist::matrix<double> imat = mist::inverse(mat);
	mist::matrix<double> iimat = mist::inverse(imat);

	printf("%+8.5lf %+8.5lf\n", imat(0,0), imat(0,1));
	printf("%+8.5lf %+8.5lf\n", imat(1,0), imat(1,1));

	printf("%+8.5lf %+8.5lf\n", iimat(0,0), iimat(0,1));
	printf("%+8.5lf %+8.5lf\n", iimat(1,0), iimat(1,1));

#endif
	return 0;
}

double配列の確保方法としては、元の参考サイトの書き方のように std::vector<double> を使うのが普通ですが、std::vectorはテンプレートライブラリなのでちょっと実行ファイルサイズが大きくなる傾向があります。 今回は実行ファイルサイズの検証をしたいため、サイズがあまり大きくならない std::malloc で書き換えてます。 行列の定義方法は、MISTも f77blas も「列優先」になっているので、そちらに合わせています。

MIST の numeric.h は、 OpenBLAS にリンクできるようにするため一部書き換えています。 (115行目付近～)。 中央の USE_OPENBLAS_LIBRARY の条件分けを追加し、コンパイルオプションに /DUSE_OPENBLAS_LIBRARY を追加しています。

// インテルのMKLとの互換性を保つための，関数名の変換マクロ
#if defined(_USE_INTEL_MATH_KERNEL_LIBRARY_) && _USE_INTEL_MATH_KERNEL_LIBRARY_ != 0
	#define LPFNAME( name ) name	// LAPACK用
	#define BLFNAME( name ) name	// BLAS用
#elif defined(USE_OPENBLAS_LIBRARY)
	#define LPFNAME( name ) name ## _	// LAPACK用
	#define BLFNAME( name ) name ## _	// BLAS用
#else
	#define LPFNAME( name ) name ## _	// LAPACK用
	#define BLFNAME( name ) f2c_ ## name	// BLAS用
#endif

処理速度比較結果

CPU: Core i7-6700、DYNAMIC_ARCH:OFF、TARGET_CORE:HASWELL での試行結果です。最適化オプションは /O2 にしています。

• openblas_set_num_threads(1) のとき

  simple: 6931.020700 msec
  cblas: 222.667100 msec
  f77blas: 221.349100 msec
  mist: 2834.766600 msec
  mist numeric: 217.342400 msec
  diff (cblas - simple): 0.001709
  diff (f77blas - simple): 0.001709
  diff (mist - simple): 0.000000
  diff (mist_numeric - simple): 0.001709
  

• openblas_set_num_threads(4) のとき

  simple: 6867.846000 msec
  cblas: 69.926100 msec
  f77blas: 69.732300 msec
  mist: 2588.654300 msec
  mist numeric: 64.678300 msec
  diff (cblas - simple): 0.001709
  diff (f77blas - simple): 0.001709
  diff (mist - simple): 0.000000
  diff (mist_numeric - simple): 0.001709
  

単純に forループで書いたもの(simple)に比べて、OpenBLAS は 4スレッド設定で約100倍、1スレッド設定で約32倍の処理速度となりました。 4スレッド設定か、1スレッド設定かは、BLAS関数を使う処理をマルチスレッド化できるかどうかで使い分ける感じが良いのでしょうか？。 MIST(numericではない方)がSIMDもマルチスレッドも使っていないことを思うと結構速いですね(OpenBLAS / MIST の速度比：4スレッド設定で約37倍、1スレッド設定で約13倍)。 OpenBLASを使わない場合はこういったライブラリを使うことになると思うので、速度比はこちらで見た方が実情に合うかもしれません。 HASWELL なので AVX2 を使って高速化しているのだと思います。AVX512 に対応しているマシンを持ってないので検証できませんが、AVX512を有効化すると一体どれくらいの速さになるのでしょうか・・・さらに2倍くらい？

実行ファイルサイズ比較結果

DYNAMIC_ARCH=OFF、OpenBLAS:スタティックリンク、Cランタイム:スタティックリンク で検証した結果です。

 * CBLASのみ有効化 : 146,432 bytes
 * F77BLASのみ有効化 : 145,920 bytes
 * CBLAS、F77BLAS を有効化 : 147,968 bytes
 * MISTのみ有効化 : 171,520 bytes
 * EXAMINE_LAPACK以外の全てを有効化 : 251,392 bytes
 * 全てを有効化 : 756,736 bytes

LAPACKを入れなければ 100kb ～ 200kb の範囲に収まりました。MISTのような C++ テンプレートライブラリを組み込むことが実行ファイルサイズの面で許容できる条件であれば、問題ない範囲かと思います。 ちなみに、DYNAMIC_ARCH:ON にすると2～4MBくらいに跳ね上がります。 OpenBLASは、DYNAMIC_ARCH:ONの場合、全てのBLAS関数をコアタイプごとにテーブルに入れて、実行時にコアタイプを自動判定してテーブルを切り替える作りになっているらしく、BLAS関数1つをリンクすると、各テーブルに入っている全てのBLAS関数を丸ごとリンクすることになる、ということでファイルサイズが2～4MBになってしまうようです。 (それでも、DLL よりは小さく済む傾向があるようですので、活かせないこともない？DLLは、コア種類を 4種くらいに絞って 9MB程度です。) LAPACK はちょっと大きいですね・・・。