为何C语言中BLAS的cblas_sgemm比NumPy的np.dot慢9倍？

我对Python NumPy与C语言OpenBLAS进行了500×500矩阵乘法的基准测试，发现np.dot的运行速度几乎是cblas_sgemm的9倍，想确认是否存在操作失误。

测试结果 #

• Python：0.001830秒
• C：0.016374秒（使用gcc -O3 -march=native编译后耗时0.016420秒）

Python代码 #

import numpy as np
import time

np.random.seed(42)
matrix_a = np.random.rand(500, 500)
matrix_b = np.random.rand(500, 500)

start_time = time.time()
result = np.dot(matrix_a, matrix_b)
end_time = time.time()

elapsed_time = end_time - start_time
print(f"{elapsed_time:.6f} seconds")

C语言代码 #

#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include <time.h>

#define N 500

int main(void) {
    srand(42);

    float *matrix_a = (float *)malloc(N * N * sizeof(float));
    float *matrix_b = (float *)malloc(N * N * sizeof(float));
    float *result = (float *)malloc(N * N * sizeof(float));

    for (int i = 0; i < N * N; ++i) {
        matrix_a[i] = (float)rand() / RAND_MAX;
        matrix_b[i] = (float)rand() / RAND_MAX;
    }

    clock_t start_time = clock();

    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                N, N, N,
                1.0, matrix_a, N, matrix_b, N,
                0.0, result, N);

    clock_t end_time = clock();

    float elapsed_time = (float)(end_time - start_time) / CLOCKS_PER_SEC;
    printf("%f seconds\n", elapsed_time);

    free(matrix_a);
    free(matrix_b);
    free(result);

    return 0;
}

编译命令 #

cc -lopenblas main.c

环境说明 #

MacOS 10.15.7，通过brew install openblas安装OpenBLAS，np.show_config()输出如下：

openblas64__info:
    libraries = ['openblas64_', 'openblas64_']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
    runtime_library_dirs = ['/usr/local/lib']
blas_ilp64_opt_info:
    libraries = ['openblas64_', 'openblas64_']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None)]
    runtime_library_dirs = ['/usr/local/lib']
openblas64__lapack_info:
    libraries = ['openblas64_', 'openblas64_']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None), ('HAVE_LAPACKE', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_ilp64_opt_info:
    libraries = ['openblas64_', 'openblas64_']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None), ('BLAS_SYMBOL_SUFFIX', '64_'), ('HAVE_BLAS_ILP64', None), ('HAVE_LAPACKE', None)]
    runtime_library_dirs = ['/usr/local/lib']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

核心问题排查方向 #

1. 数据类型不匹配 #

NumPy的np.random.rand默认生成双精度浮点数（float64），而C代码使用的是单精度浮点数（float32）。np.dot会用float64计算，cblas_sgemm是单精度实现，理论上单精度运算更快，但这里结果相反，大概率是NumPy链接的OpenBLAS对float64做了更充分的SIMD指令优化，而你的C代码单精度版本未利用到最优指令集。

2. OpenBLAS版本差异（ILP64 vs LP64） #

从np.show_config()可知，NumPy使用的是ILP64版本的OpenBLAS（openblas64_），该版本用64位整数做索引，大矩阵下性能更优。而你编译C代码时默认链接的是LP64版本，索引方式的差异会直接导致性能差距。

3. 编译参数与库链接正确性 #

即使加了-O3 -march=native，如果未显式链接ILP64版本的OpenBLAS，性能无法和NumPy对齐。编译时需指定链接openblas64库：

gcc -O3 -march=native main.c -lopenblas64

4. 计时方式不一致 #

Python用time.time()统计墙钟时间，C代码用clock()统计CPU时间。如果测试C代码时系统有其他进程占用CPU，clock()会累加所有CPU耗时，导致结果偏高。建议C代码改用墙钟时间计时：

#include <time.h>

struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
// 执行矩阵乘法
clock_gettime(CLOCK_MONOTONIC, &end);
double elapsed = (end.tv_sec - start.tv_sec) + (end.tv_nsec - start.tv_nsec) / 1e9;
printf("%f seconds\n", elapsed);

5. 内存对齐问题 #

NumPy数组会自动做内存对齐（通常64字节），OpenBLAS对对齐内存的访问有明显优化。而C代码中malloc的内存可能未对齐，可改用posix_memalign分配对齐内存：

float *matrix_a;
posix_memalign((void**)&matrix_a, 64, N*N*sizeof(float));

验证建议 #

• 修改C代码，将float替换为double，调用cblas_dgemm（双精度版本）重新测试；
• 显式链接ILP64版本的OpenBLAS编译C代码；
• 对齐C代码的计时方式与Python一致；
• 使用内存对齐的分配方式。

内容的提问来源于stack exchange，提问作者Momo