为何MKL SGEMM性能低于自研分块OpenMP实现与Eigen？

我正在测试三种float类型n×n方阵的矩阵乘法实现性能：朴素分块OpenMP实现、Eigen库、MKL 2021.4.0的SGEMM。所有矩阵按64字节对齐，编译环境为GCC 8.3.1，编译参数-msse4.2 -O3 -fopenmp，操作系统CentOS 7。

但测试结果显示MKL SGEMM性能最差，甚至不如朴素OpenMP实现，完全不符合预期，求解释。

分块OpenMP实现（BS = n / 64） #

#pragma omp for collapse(2)
for(int i=0; i<n; i++)
  for(int j=0; j<n; j++)
    C[i*n+j] *= beta;

#pragma omp parallel for schedule(dynamic)
for (int i = 0; i < n; i+=BS) 
for (int k = 0; k < n; k+=BS) 
for (int j = 0; j < n; j+=BS) 
    for (int ii = i; ii < i+BS; ii++)
    for (int kk = k; kk < k+BS; kk++)
    for (int jj = j; jj < j+BS; jj++)
        C[ii*n+jj] += alpha*A[ii*n+kk]*B[kk*n+jj];

Eigen实现 #

Eigen::Map<const Eigen::MatrixXf> AM(A, n, n);
Eigen::Map<const Eigen::MatrixXf> BM(B, n, n);
Eigen::Map<Eigen::MatrixXf> CM(C, n, n);
CM.noalias() = beta*CM + alpha*(BM * AM); // fortran order!

MKL SGEMM实现 #

cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
    n, n, n, alpha, A, n, B, n, beta, C, n);

测试环境 #

Intel Xeon Silver 4114（2路插槽、2个NUMA节点）

测试结果 #

Google Benchmark测试结果 #

Benchmark                                     Time             CPU   Iterations
-------------------------------------------------------------------------------
MatMul/OmpBlk/4096/64/real_time            1132 ms         1038 ms            1
MatMul/OmpBlk/16384/64/real_time          83668 ms        80612 ms            1
MatMul/OmpBlk/32768/64/real_time        1562980 ms      1492184 ms            1
MatMul/Eigen/4096/real_time                 878 ms          867 ms            1
MatMul/Eigen/16384/real_time              36140 ms        31629 ms            1
MatMul/Eigen/32768/real_time             259762 ms       246788 ms            1
MatMul/Blas/4096/real_time                 4091 ms         3719 ms            1
MatMul/Blas/16384/real_time              219940 ms       219581 ms            1
MatMul/Blas/32768/real_time             1773874 ms      1750015 ms            1

三次运行平均时间（无预热，未用Google Benchmark） #

OmpBlk/4096:    1452 ms
OmpBlk/16384:  87494 ms
Eigen/4096:      818 ms
Eigen/16384:   34719 ms
Blas/4096:      4060 ms
Blas/16384:   225647 ms

依赖库片段（ldd输出） #

libmkl_intel_ilp64.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_intel_ilp64.so.1 
libmkl_core.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_core.so.1 
libmkl_intel_thread.so.1 => /opt/intel/oneapi/mkl/2021.4.0/lib/intel64/libmkl_intel_thread.so.1 
libiomp5.so => /opt/intel/oneapi/compiler/latest/linux/compiler/lib/intel64/libiomp5.so 

可能的性能异常原因分析 #

1. ILP64接口与数据类型不匹配
   你使用了MKL的ILP64版本库，但代码中矩阵维度n是32位int类型。ILP64要求所有整数参数（维度、步长等）为64位long，传入32位int可能导致参数解析错误，触发低效路径或额外类型转换开销。建议切换到LP64版本的MKL库，或把相关参数改为long类型。

2. 线程库冲突
   编译时用了GCC的-fopenmp（绑定libgomp），但MKL链接的是Intel OpenMP库libiomp5.so。两种OpenMP库共存会导致线程调度混乱，比如MKL线程池与应用层OpenMP线程竞争CPU资源，严重影响并行效率。解决方法：统一使用一种OpenMP库，要么编译时添加-liomp5并去掉-fopenmp，要么设置环境变量MKL_THREADING_LAYER=GNU让MKL使用GNU OpenMP。

3. NUMA节点亲和性问题
   双NUMA节点架构下，若矩阵数据未绑定到对应NUMA节点，MKL会出现跨NUMA访问的高延迟。而你的OpenMP实现可能通过动态调度让线程更贴近本地内存。建议用numactl --cpunodebind=0 --membind=0绑定到单个NUMA节点测试，或在代码中显式设置线程NUMA亲和性。

4. MKL参数与矩阵布局不匹配
   你给cblas_sgemm指定了CblasRowMajor（行优先），但Eigen代码注释了fortran order（列优先），需确认矩阵A/B的实际存储布局是否与MKL参数匹配。若布局不匹配，MKL会进行隐式转置，带来巨大性能开销。

5. MKL运行时环境未优化
   默认情况下MKL可能未充分利用CPU特性，建议设置以下环境变量：

   • MKL_ENABLE_INSTRUCTIONS=AVX2（Xeon Silver 4114支持AVX2）
   • MKL_NUM_THREADS：设置为CPU核心数（双路共40线程），避免线程过载。

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