はじめに

「行列積ってどの実装が速いの？」という素朴な疑問から、C++ 4種と Rust で同じ N×N 行列積を実装してベンチマーク比較しました。

比較対象は以下の5種です：

実装	方式
`naive`	C++ i-k-j ループ（キャッシュフレンドリー）
`openmp`	C++ + OpenMP 並列化
`eigen`	C++ + Eigen::MatrixXd
`openblas`	C++ + cblas_dgemm（DGEMM）
`rust`	pure Rust i-k-j ループ

Dev Container（Ubuntu 22.04）内で完結しており、bash benchmark/run_all.sh 1コマンドで全5実装を一括計測できます。

環境

項目	バージョン
OS	Ubuntu 22.04 (Dev Container)
g++	11.4.0
cmake	3.22.1
Eigen	3.4.0
OpenBLAS	0.3.20
rustc / cargo	1.94.1

実装概要

ベンチマーク条件

行列サイズ: 512×512 / 1024×1024 / 2048×2048
各サイズを 3回計測して中央値を採用
乱数シードは seed=42 で固定（再現性確保）
単位: ミリ秒（ms）

共通の出力フォーマット

[実装名]  N=512   median=XXX.XX ms
[実装名]  N=1024  median=XXX.XX ms
[実装名]  N=2048  median=XXX.XX ms

各実装の解説

1. C++ naive（i-k-j ループ）

最もシンプルな実装です。ループ順を i-k-j にすることで B のアクセスが行方向に連続し、キャッシュフレンドリーになります（i-j-k より大幅に速い）。

static void matmul(const std::vector<double>& A,
                   const std::vector<double>& B,
                   std::vector<double>&       C,
                   int N)
{
    for (int i = 0; i < N; ++i) {
        for (int k = 0; k < N; ++k) {
            double a_ik = A[i * N + k];   // ← キャッシュに乗せる
            for (int j = 0; j < N; ++j) {
                C[i * N + j] += a_ik * B[k * N + j];
            }
        }
    }
}

2. C++ + OpenMP

naive の外側ループ（i）に #pragma omp parallel for を付けるだけで並列化できます。

#pragma omp parallel for schedule(static)
for (int i = 0; i < N; ++i) {
    for (int k = 0; k < N; ++k) {
        double a_ik = A[i * N + k];
        for (int j = 0; j < N; ++j) {
            C[i * N + j] += a_ik * B[k * N + j];
        }
    }
}

実行時のスレッド数は omp_get_max_threads() で確認できます（今回は4スレッド）。

3. C++ + Eigen

Eigen::MatrixXd の演算子 * で呼び出すだけ。Eigen は内部で SIMD 最適化を自動適用します。

Eigen::MatrixXd A(N, N), B(N, N);
// ... 初期化 ...

auto start = std::chrono::high_resolution_clock::now();
Eigen::MatrixXd C = A * B;   // ← SIMD 最適化が入る
auto end   = std::chrono::high_resolution_clock::now();

コードが非常にシンプルで、しかも高速という優れた選択肢です。

4. C++ + OpenBLAS（cblas_dgemm）

BLAS の DGEMM 関数を直接呼び出します。C = α·A·B + β·C という一般形で、α=1.0, β=0.0 に設定。

// C = 1.0 * A * B + 0.0 * C
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            N, N, N,
            1.0, A.data(), N,
                 B.data(), N,
            0.0, C.data(), N);

引数が多く見た目は複雑ですが、パフォーマンスは圧倒的です。

5. pure Rust（i-k-j ループ）

C++ naive と同アルゴリズムを Rust で実装。std::hint::black_box でコンパイラの最適化除去を防いでいます。

fn matmul(a: &[f64], b: &[f64], c: &mut [f64], n: usize) {
    for i in 0..n {
        for k in 0..n {
            let a_ik = a[i * n + k];
            for j in 0..n {
                c[i * n + j] += a_ik * b[k * n + j];
            }
        }
    }
}

// 計測時は black_box でラップ
matmul(black_box(&a), black_box(&b), black_box(&mut c), n);

ベンチマーク結果

実装	N=512 (ms)	N=1024 (ms)	N=2048 (ms)
naive	17.94	161.95	2475.91
openmp (4threads)	7.70	76.08	1508.45
eigen	9.87	59.19	463.94
openblas	1.84	31.25	200.53
rust	25.47	256.32	3750.35

速度傾向

OpenBLAS >> Eigen >> OpenMP >> naive ≈ Rust(pure)

N=2048 での対 naive 比:

実装	naive 比
OpenBLAS	12.3× 高速
Eigen	5.3× 高速
OpenMP (4threads)	1.6× 高速
naive	1.0×（基準）
Rust	0.66×（低速）

考察

OpenBLAS が圧倒的。高度に最適化された BLAS 実装の力を見せつける結果に。
Eigen は使いやすさと速さのバランスが良く、コードもシンプルで実用的。
OpenMP は 4スレッドでも naive の 1.6× 止まり。N=2048 でも naive と同じアルゴリズムなので、スレッドオーバーヘッドが効いている。
Rust は C++ naive より遅いという結果に。LCG による擬似乱数生成が MT19937_64 より遅い可能性があるほか、LTO を有効にしても C++ の -march=native による SIMD 最適化には届かなかった。

ハマりどころ

libopenblas-dev が apt でインストールできない

Ubuntu 22.04 の Dev Container 環境では libopenblas-dev が apt-get install できない問題に遭遇しました。

$ sudo apt-get install -y libopenblas-dev
E: Package 'libopenblas-dev' has no installation candidate

解決方法: libopenblas0 は既にインストール済みで共有ライブラリは存在しています。.so シンボリックリンクを手動で作成することで解決できました。

sudo ln -sf /usr/lib/x86_64-linux-gnu/libopenblas.so.0 \
            /usr/lib/x86_64-linux-gnu/libopenblas.so

cblas.h は libblas-dev に付属する /usr/include/x86_64-linux-gnu/cblas.h を使用します。

Rust の最適化除去

Rust はリリースビルドで積極的に最適化するため、何も返さない計算はコンパイラに丸ごと除去される可能性があります。std::hint::black_box を使ってこれを防ぎました。

use std::hint::black_box;

matmul(black_box(&a), black_box(&b), black_box(&mut c), n);
let _ = black_box(c[0]); // 結果を参照して除去防止

まとめ

最速は OpenBLAS。N=2048 で naive の 12.3 倍高速。BLAS を使える場面では積極的に採用するべき。
Eigen は「速さ × 使いやすさ」のベストバランス。A * B と書くだけで SIMD 最適化が入る。
pure Rust は C++ naive と同等かやや遅い。Rust でも BLAS バインディング（blas crate）を使えば OpenBLAS と同等の速度が出るはず。
OpenMP は手軽だが劇的な高速化は難しい。アルゴリズム改善（BLAS / Eigen）の方が効果的。

参考