C++ OpenMP矩阵乘法无性能提升？如何缩小与NumPy的差距

云上老周

2026-7-11

问题：C++矩阵乘法性能远低于NumPy，寻求优化方案

我是C++新手，写了一个能正常运行的矩阵乘法程序，但基准测试显示速度比NumPy慢很多。一开始用OpenMP加速时没加/openmp编译选项，性能没变化；加上之后速度提升到原版本的1/4，但还是比NumPy慢一个数量级。以下是我的代码、编译参数和测试结果，求进一步缩小性能差距的方法。

代码

#include <algorithm>
#include <chrono>
#include <iostream>
#include <omp.h>
#include <string>
#include <vector>

using std::vector;
using std::chrono::high_resolution_clock;
using std::chrono::duration;
using std::chrono::duration_cast;
using std::chrono::microseconds;
using std::cout;
using line = vector<double>;
using matrix = vector<line>;

void fill(line &l) {
    std::generate(l.begin(), l.end(), []() { return ((double)rand() / (RAND_MAX)); });
}

matrix random_matrx(int64_t height, int64_t width) {
    matrix mat(height, line(width));
    std::for_each(mat.begin(), mat.end(), fill);
    return mat;
}

matrix dot_product(const matrix &mat0, const matrix &mat1) {
    size_t h0, w0, h1, w1;
    h0 = mat0.size();
    w0 = mat0[0].size();
    h1 = mat1.size();
    w1 = mat1[0].size();
    if (w0 != h1) {
        throw std::invalid_argument("matrices cannot be cross multiplied");
    }

    matrix out(h0, line(w1));
    for (int y = 0; y < h0; y++) {
        for (int x = 0; x < w1; x++) {
            double s = 0;
            for (int z = 0; z < w0; z++) {
                s += mat0[y][z] * mat1[z][x];
            }
            out[y][x] = s;
        }
    }

    return out;
}

matrix dot_product_omp(const matrix& mat0, const matrix& mat1) {
    size_t h0, w0, h1, w1;
    h0 = mat0.size();
    w0 = mat0[0].size();
    h1 = mat1.size();
    w1 = mat1[0].size();
    if (w0 != h1) {
        throw std::invalid_argument("matrices cannot be cross multiplied");
    }

    matrix out(h0, line(w1));
    omp_set_num_threads(4);
    #pragma omp parallel for schedule(dynamic)
    for (int y = 0; y < h0; y++) {
        for (int x = 0; x < w1; x++) {
            double s = 0;
            for (int z = 0; z < w0; z++) {
                s += mat0[y][z] * mat1[z][x];
            }
            out[y][x] = s;
        }
    }

    return out;
}

int main()
{
    matrix a, b;
    a = random_matrx(16, 9);
    b = random_matrx(9, 24);
    auto start = high_resolution_clock::now();
    for (int64_t i = 0; i < 65536; i++) {
        dot_product(a, b);
    }
    auto end = high_resolution_clock::now();
    duration<double, std::nano> time = end - start;
    double once = time.count() / 65536000;
    cout << "mat(16, 9) * mat(9, 24): " + std::to_string(once) + " microseconds\n";
    a = random_matrx(128, 256);
    b = random_matrx(256, 512);
    start = high_resolution_clock::now();
    for (int64_t i = 0; i < 512; i++) {
        dot_product(a, b);
    }
    end = high_resolution_clock::now();
    time = end - start;
    once = time.count() / 512000;
    cout << "mat(128, 256) * mat(256, 512): " + std::to_string(once) + " microseconds\n";
    start = high_resolution_clock::now();
    for (int64_t i = 0; i < 512; i++) {
        dot_product_omp(a, b);
    }
    end = high_resolution_clock::now();
    time = end - start;
    once = time.count() / 512000;
    cout << "mat(128, 256) * mat(256, 512) omp: " + std::to_string(once) + " microseconds\n";
}

初始测试结果（未加`/openmp`）

PS D:\MyScript> C:\Users\Xeni\source\repos\matmul\x64\Release\matmul.exe
mat(16, 9) * mat(9, 24): 5.200116 microseconds
mat(128, 256) * mat(256, 512): 30128.739453 microseconds
mat(128, 256) * mat(256, 512) omp: 30116.103125 microseconds

编译参数

使用Visual Studio 2022、C++20标准，编译主参数：

/permissive- /ifcOutput "x64\Release\" /GS /GL /W3 /Gy /Zc:wchar_t /Zi /Gm- /O2 /Ob2 /sdl /Fd"x64\Release\vc143.pdb" /Zc:inline /fp:precise /D "NDEBUG" /D "_CONSOLE" /D "_UNICODE" /D "UNICODE" /errorReport:prompt /WX- /Zc:forScope /std:c17 /Gd /Oi /MD /std:c++20 /FC /Fa"x64\Release\" /EHsc /nologo /Fo"x64\Release\" /Ot /Fp"x64\Release\matmul.pch" /diagnostics:column

附加参数：/arch:AVX2 /fp:fast

添加`/openmp`后的测试结果

PS D:\MyScript> C:\Users\Xeni\source\repos\matmul\x64\Release\matmul.exe
mat(16, 9) * mat(9, 24): 5.126476 microseconds
mat(128, 256) * mat(256, 512): 30999.137109 microseconds
mat(128, 256) * mat(256, 512) omp: 8574.475195 microseconds

NumPy测试结果

In [374]: a = np.random.random((128, 256))
In [375]: b = np.random.random((256, 512))
In [376]: %timeit a @ b
382 µs ± 19.6 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

优化方案

1. 优化内存布局与循环顺序，解决缓存命中问题

当前用vector<vector<double>>的嵌套结构存储矩阵，访问mat1[z][x]时是按列遍历，会导致大量CPU缓存失效——CPU缓存按连续内存块加载，列遍历会跳过大量缓存行，这是性能瓶颈的核心原因。NumPy用连续一维内存存储矩阵，且通过循环重排保证缓存友好的访问模式。

改进方式：

调整循环顺序：将原有的y->x->z改为y->z->x，让对mat1的访问变为连续行访问
（可选）改用连续内存存储矩阵，比如用vector<double>模拟二维矩阵，通过y * width + x计算索引

修改后的核心乘法循环示例：

matrix dot_product_optimized(const matrix& mat0, const matrix& mat1) {
    size_t h0 = mat0.size();
    size_t w0 = mat0[0].size();
    size_t w1 = mat1[0].size();
    matrix out(h0, line(w1, 0.0));

    #pragma omp parallel for
    for (int y = 0; y < h0; y++) {
        const auto& row0 = mat0[y];
        auto& row_out = out[y];
        for (int z = 0; z < w0; z++) {
            double val = row0[z];
            const auto& row1 = mat1[z];
            for (int x = 0; x < w1; x++) {
                row_out[x] += val * row1[x];
            }
        }
    }
    return out;
}

2. 显式引导编译器生成SIMD指令

虽然开启了/arch:AVX2，但编译器自动向量化可能因循环结构限制效果不佳，可通过显式指令优化：

在内存访问连续的内层循环添加#pragma omp simd，强制编译器生成AVX2向量指令
确保循环内无复杂分支，变量依赖关系清晰

示例：

#pragma omp simd
for (int x = 0; x < w1; x++) {
    row_out[x] += val * row1[x];
}

3. 调整OpenMP并行策略

移除手动设置omp_set_num_threads(4)，让OpenMP自动使用系统全部核心，提升并行规模
将schedule(dynamic)改为schedule(static)，动态调度会带来额外开销，静态调度更适合计算量均匀的矩阵乘法

修改后的并行区域：

#pragma omp parallel for schedule(static)

4. 减少内存分配开销

每次调用乘法函数都创建新矩阵，频繁的内存分配释放会产生额外开销。改为预先分配输出矩阵，重复使用：

void dot_product_inplace(const matrix& mat0, const matrix& mat1, matrix& out) {
    size_t h0 = mat0.size();
    size_t w0 = mat0[0].size();
    size_t w1 = mat1[0].size();

    #pragma omp parallel for schedule(static)
    for (int y = 0; y < h0; y++) {
        std::fill(out[y].begin(), out[y].end(), 0.0);
        const auto& row0 = mat0[y];
        auto& row_out = out[y];
        for (int z = 0; z < w0; z++) {
            double val = row0[z];
            const auto& row1 = mat1[z];
            #pragma omp simd
            for (int x = 0; x < w1; x++) {
                row_out[x] += val * row1[x];
            }
        }
    }
}

测试时预先创建好输出矩阵，避免重复分配。

5. 直接调用BLAS库

NumPy底层依赖优化后的BLAS库（如OpenBLAS、Intel MKL），这些库利用了CPU全特性（AVX-512、多线程、缓存优化等），性能远超手动实现。可以直接调用BLAS的dgemm函数：

示例（Intel MKL）：

#include <mkl.h>

void dot_product_mkl(const matrix& mat0, const matrix& mat1, matrix& out) {
    size_t h0 = mat0.size();
    size_t w0 = mat0[0].size();
    size_t w1 = mat1[0].size();
    const double alpha = 1.0;
    const double beta = 0.0;
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                h0, w1, w0, alpha,
                mat0[0].data(), w0,
                mat1[0].data(), w1,
                beta, out[0].data(), w1);
}

需配置MKL的头文件和链接库路径。

6. 编译参数微调

确保/fp:fast开启，允许编译器进行浮点优化，精度和NumPy保持一致
添加/Qvec-report:2编译参数，查看编译器向量化报告，针对未成功向量化的循环调整结构

内容的提问来源于stack exchange，提问作者Ξένη Γήινος

C++ OpenMP矩阵乘法无性能提升？如何缩小与NumPy的差距

代码 #

初始测试结果（未加/openmp） #

编译参数 #

添加/openmp后的测试结果 #

NumPy测试结果 #

1. 优化内存布局与循环顺序，解决缓存命中问题 #

改进方式： #

2. 显式引导编译器生成SIMD指令 #

3. 调整OpenMP并行策略 #

4. 减少内存分配开销 #

5. 直接调用BLAS库 #

6. 编译参数微调 #

代码