调用Fortran Intel MKL dgemm的函数性能远低于numpy及Fortran matmul的问题排查

算力阿峰

2026-4-29

分析与解决方案

首先，你的问题核心在于Fortran与numpy的内存布局不匹配，导致MKL的DGEMM无法高效利用缓存，性能暴跌。另外，MKL多线程的配置也可能需要调整，我们一步步来解决：

1. 内存布局的关键问题

Fortran默认使用列优先（Column-Major）存储数组，而numpy默认是行优先（Row-Major）。当你把numpy的行优先数组传给Fortran子程序时，Fortran会按照列优先的逻辑去解读内存——这意味着数组的实际存储顺序和Fortran声明的维度完全不匹配。

举个例子：你的M1(M,N)在Fortran里应该是按列存储（先存第1列所有元素，再第2列...），但numpy的行优先数组是按行存储（先存第1行所有元素，再第2行...）。DGEMM是高度缓存优化的函数，需要连续的内存访问，这种错位会导致DGEMM每次访问内存都跳过大块区域，缓存命中率几乎为0，直接让性能下降几百倍。

而numpy的dot函数内部会自动处理内存布局（要么转置数组，要么调用适配行优先的MKL接口），所以能保持高性能。

解决方法：让numpy数组使用列优先存储

修改Python里数组的创建代码，指定order='F'：

# Setup
M=500
N=500
K=500
# 改为列优先存储，匹配Fortran
a = np.empty((M,N), dtype=c_double, order='F')
b = np.empty((N,K), dtype=c_double, order='F')
c = np.empty((M,K), dtype=c_double, order='F')
a[:] = np.random.rand(M,N)
b[:] = np.random.rand(N,K)

2. 确保MKL多线程生效

你的编译脚本已经链接了MKL的线程库，但还需要在编译时启用OpenMP，让MKL能利用多核CPU：

修改bat脚本中的IFORT_OPTIMIZATION_FLAGS，添加-qopenmp：

SET "IFORT_OPTIMIZATION_FLAGS=/O3 -qopenmp"

另外，你可以在Python脚本里显式设置MKL的线程数，确保充分利用CPU资源：

import os
# 设置为你的CPU核心数，比如8
os.environ['MKL_NUM_THREADS'] = str(os.cpu_count())

或者在Fortran代码中调用MKL的线程设置函数（需要添加use mkl_service）：

subroutine matmultmkl(M1, M2, M3, M, N, K) bind(c, name='matmultmkl')
!DEC$ ATTRIBUTES DLLEXPORT :: matmultmkl
use iso_c_binding, only: c_double, c_int
use mkl_service, only: mkl_set_num_threads
integer(c_int),intent(in) :: M, N, K
real(c_double), intent(in) :: M1(M, N), M2(N, K)
real(c_double), intent(inout):: M3(M, K)
! 设置线程数（根据你的CPU调整）
call mkl_set_num_threads(8)
CALL DGEMM('N','N',M,K,N,1.,M1,M,M2,N,0.,M3,M)
end subroutine

3. 验证结果正确性

之前的内存布局错误不仅影响性能，还会导致计算结果错误！调整完布局后，你可以对比Fortran MKL版本和numpy的计算结果，确认是否一致：

# 计算后验证
c_mkl = c.copy()
c_numpy = a.dot(b)
print("结果是否一致：", np.allclose(c_mkl, c_numpy))

预期效果

调整后，Fortran MKL版本的性能应该和numpy处于同一量级（甚至可能因为少了numpy的封装开销略快一点），你会看到时间从0.5s左右降到0.001~0.002s的水平。

额外提示

如果你确实需要处理行优先的numpy数组，也可以在Fortran里调整DGEMM的TRANSA和TRANSB参数为'T'，相当于告诉MKL数组是转置后的，但这种方法不如直接用列优先数组高效。
编译时可以添加-xHost选项，让编译器针对你的CPU架构做深度优化，进一步提升性能：
```
SET "IFORT_OPTIMIZATION_FLAGS=/O3 -qopenmp -xHost"
```

内容的提问来源于stack exchange，提问作者Olórin

调用Fortran Intel MKL dgemm的函数性能远低于numpy及Fortran matmul的问题排查

1. 内存布局的关键问题 #

解决方法：让numpy数组使用列优先存储 #

2. 确保MKL多线程生效 #

3. 验证结果正确性 #

预期效果 #