前回（第1弾・第2弾）では、C++ 4種 + Rust 4種の行列積ベンチマーク比較を行い、OpenBLAS が圧倒的に速いという結論に至った。
では、そのOpenBLASをさらに速くすることはできるのか？
AMD Ryzen 5 7530U（Zen 3 / WSL2環境）で3つのアプローチを試した。
環境
- CPU: AMD Ryzen 5 7530U（Zen 3、6コア12スレッド）
- 実行環境: WSL2（2コア4スレッドが割り当て）
- OS: Ubuntu on WSL2
- コンパイラ: GCC、-O3 -march=native
- OpenBLAS: 0.3.20（ZENカーネル入り）
- 行列サイズ: N = 512 / 1024 / 2048（double精度、各3回中央値）
ベースライン（前回結果）：
| 実装 | N=2048 |
|---|---|
| OpenBLAS DGEMM（C++） | 206 ms |
| OpenBLAS cblas FFI（Rust） | 178 ms |
アプローチ 1：OpenBLASのスレッド数チューニング
動機
OpenBLASはデフォルトで利用可能なコア数を使う。しかしスレッド数が多ければ速いとは限らない。WSL2のオーバーヘッド・ハイパーバイザー越しのスケジューリングにより、スレッド競合が発生する可能性がある。
openblas_set_num_threads() でスレッド数を明示的に制御して計測した。
extern "C" void openblas_set_num_threads(int num_threads);
// 計測前に設定
openblas_set_num_threads(threads);
結果
| スレッド数 | N=512 | N=1024 | N=2048 |
|---|---|---|---|
| 1 | 11.8 ms | 49.7 ms | 386.1 ms |
| 2 | 3.4 ms | 28.3 ms | 229.2 ms |
| 4 | 1.9 ms | 28.0 ms | 170.9 ms |
考察
4スレッドで 170.9ms（従来C++比 1.21倍速）。
デフォルト設定（C++の206ms）との差が生まれた理由は、前回のC++実装がスレッド数を指定していなかったためで、OpenBLASの自動検出が最適でなかった可能性がある。
ポイント:
openblas_set_num_threads()を明示すると再現性が上がり、WSL2では4スレッドが現状最速。
アプローチ 2：Eigen の BLAS バックエンド（EIGEN_USE_BLAS）
動機
前回、EigenはOpenBLASより約2.3倍遅かった（476ms vs 206ms）。
EIGEN_USE_BLAS を定義すると、EigenのGEMM計算を外部BLASに委譲できる。
「Eigenの使いやすいAPIを保ちつつ、OpenBLASの速度を得られるか？」を検証した。
// cblas.h の前に定義する
#define EIGEN_USE_BLAS
#include <x86_64-linux-gnu/cblas.h>
#include <Eigen/Dense>
// 以降は通常の Eigen コード
Eigen::MatrixXd C = A * B; // 内部で cblas_dgemm が呼ばれる
CMakeLists.txt での設定：
add_executable(bench_eigen_blas src/matmul_eigen_blas.cpp)
target_include_directories(bench_eigen_blas PRIVATE ${CBLAS_INCLUDE_DIR})
target_link_libraries(bench_eigen_blas PRIVATE Eigen3::Eigen ${OPENBLAS_LIB})
結果
| 実装 | N=512 | N=1024 | N=2048 |
|---|---|---|---|
| Eigen（純粋） | 9.4 ms | 59.1 ms | 475.6 ms |
| Eigen + BLAS | 4.4 ms | 23.3 ms | 229.8 ms |
| OpenBLAS直接（t4） | 1.9 ms | 28.0 ms | 170.9 ms |
考察
EIGEN_USE_BLAS でEigenは約2倍速くなったが、OpenBLAS直接呼び出しには届かなかった（229ms vs 171ms）。
理由として考えられるのは：
- Eigen が内部でスレッド数を制御せず、デフォルト設定のまま動く
- 行列配置（列優先 vs 行優先）の変換コスト
ポイント:
EIGEN_USE_BLASは純粋Eigenの倍速だが、直接cblas呼び出しの代替にはならない。「Eigenで書いてBLASで動かす」用途に限定すべき。
アプローチ 3：float32 SGEMM（精度を落として速度を得る）
動機
AVX2 は 256bit幅のSIMDレジスタを持つ。
-
double（64bit）: 1レジスタに 4要素 -
float（32bit）: 1レジスタに 8要素 ← 2倍！
理論上、cblas_sgemm（単精度）は cblas_dgemm（倍精度）の 2倍のスループット を発揮できる。
精度が落ちるというトレードオフがあるが、機械学習推論や物理シミュレーションではfloat32で十分な場合が多い。
cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
N, N, N,
1.0f, A.data(), N, // float*
B.data(), N,
0.0f, C.data(), N);
結果
| 実装 | N=512 | N=1024 | N=2048 |
|---|---|---|---|
| OpenBLAS DGEMM（t4） | 1.9 ms | 28.0 ms | 170.9 ms |
| OpenBLAS SGEMM | 3.3 ms | 12.4 ms | 104.8 ms |
考察
N=2048 で 104.8ms。DGEMM比で 約1.97倍速。ほぼ理論値の2倍を達成。
N=512ではSGEMMが若干遅い（3.3ms vs 1.9ms）。小行列ではSIMD効率よりメモリ帯域・起動コストが支配的になるためと考えられる。
ポイント: 精度を float32 に落とせる用途なら SGEMM で約2倍の高速化が可能。機械学習・グラフィクス・科学技術計算の多くはfloat32で十分。
総合比較（N=2048）
| 実装 | N=2048 (ms) | 対DGEMM基準 |
|---|---|---|
| Eigen 純粋 | 475.6 | 0.43x |
| OpenBLAS DGEMM（1スレッド） | 386.1 | 0.53x |
| OpenBLAS DGEMM（2スレッド） | 229.2 | 0.90x |
| Eigen + BLAS backend | 229.8 | 0.90x |
| OpenBLAS DGEMM（従来C++） | 206.5 | 基準 |
| OpenBLAS DGEMM（4スレッド明示） | 170.9 | 1.21x |
| 🏆 OpenBLAS SGEMM（float32） | 104.8 | 1.97x |
まとめ
「BLASより速くできるか？」への回答：
| アプローチ | 効果 | 推奨度 |
|---|---|---|
| スレッド数チューニング | 1.21x 速化 | ⭐⭐⭐ 常に試す価値あり |
| Eigen + BLAS backend | Eigen比 2x、BLAS直接比 0.9x | ⭐⭐ Eigenの利便性を保つ用途向け |
| SGEMM（float32） | 1.97x 速化 | ⭐⭐⭐ 精度許容できるなら最有力 |
最も効果的な答えは「float32 SGEMM」。
AMDのAVX2実装がSIMD幅を余すことなく使い切り、理論値に近い2倍を達成した。
次の選択肢として以下が考えられる：
-
AOCL（AMD Optimizing CPU Libraries）：AMD公式のBLIS実装でZen向けに最適化されたCPUライブラリ。WSL2でも動作するが、AMD公式サイトでのアカウント登録・EULA同意が必要なためコンテナ環境では自動インストールが困難。ホスト環境であれば
.debを手動インストールして計測できる。 - ROCm / rocBLAS（GPU）：WSL2からAMD iGPU（Radeon Vega）へのROCmアクセスは非対応のため今回は断念。ホスト環境かつ専用GPU環境であれば有効な選択肢。
コード
（src/matmul_openblas_mt.cpp, src/matmul_eigen_blas.cpp, src/matmul_sgemm.cpp を参照）