使用 fortran intel 的 mkl dgemm 和 python ctypes 时出现分段错误，矩阵大小已经非常小

如何解决使用 fortran intel 的 mkl dgemm 和 python ctypes 时出现分段错误，矩阵大小已经非常小

这是对以下内容的跟进：

Dll dependencies missing when using ctypes with fortran dll compiled with intel fortran compiler and linking to intel's mkl library

我还在窗户下。 mkl_example.f 中的 Fortran 代码：

    subroutine matmultmkl(M1,M2,M3,M,N,K) bind(c,name='matmultmkl')
        !DEC$ ATTRIBUTES DLLEXPORT :: matmultmkl
        use iso_c_binding,only: c_float,c_int
    
        integer(c_int),intent(in) :: M,K
        real(c_float),intent(in) :: M1(M,N),M2(N,K)
        real(c_float),intent(out):: M3(M,K)

        CALL DGEMM('N','N',K,1.,M1,0.,M)
    
    end subroutine 

还是一样。我使用以下批处理文件使用命令行（在我之前运行 compilervars.bat 的 cmd 窗口中）编译它：

@Echo off

setlocal ENABLEDELAYEDEXPANSION
SET "IFORT_INITIAL_FLAGS=-c -fpp"
SET "IFORT_OPTIMIZATION_FLAGS=/O3"

ifort %IFORT_OPTIMIZATION_FLAGS% %IFORT_INITIAL_FLAGS% /I"C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2020.4.311\windows\mkl\include" -o mkl_example.obj mkl_example.f
ifort -dll -o mylib.dll mkl_example.obj /link /LIBPATH:"C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2020.4.311\windows\mkl\lib\intel64_win" mkl_intel_lp64.lib mkl_intel_thread.lib mkl_core.lib libiomp5md.lib 

也没有改变。然后我运行以下 python 脚本 mkl_example.py，我将其更改为指向静态库和导入库：

import os
os.add_dll_directory(r"C:/Program Files (x86)/IntelSWTools/compilers_and_libraries_2020.4.311/windows/redist/intel64_win/mkl")
os.add_dll_directory(r"C:/Program Files (x86)/IntelSWTools/compilers_and_libraries_2020.4.311/windows/redist/intel64_win/compiler")
import time
from ctypes import *
import signal
import numpy as np
mylib = CDLL(r"C:/path/to/the/mylib.dll")

def main():
    mylib.matmultmkl.argtypes = [ POINTER(c_float),POINTER(c_float),POINTER(c_int),POINTER(c_int) ]

    mylib.matmultmkl.restype = None

    M=231
    N=231
    K=231

    a = np.empty((M,dtype=c_float)
    b = np.empty((N,K),dtype=c_float)
    c = np.empty((M,dtype=c_float)

    a[:] = np.random.rand(M,N)
    b[:] = np.random.rand(N,K)
    c[:] = np.full((M,0.0)

    # Fortran mkl call
    start = time.time()
    mylib.matmultmkl( a.ctypes.data_as(POINTER(c_float)),b.ctypes.data_as(POINTER(c_float)),c.ctypes.data_as(POINTER(c_float)),c_int(M),c_int(N),c_int(K) )
    stop = time.time()
    print(f"Fortran mkl \t {stop - start}s")

if __name__ == "__main__":
    main()

dgemm 执行矩阵乘积。在这里，我有两个大小为 231 的方阵，这对于今天的矩阵乘积来说是一个完全愚蠢的小尺寸，但是：对 fortran 函数的调用触发了分段错误。如果大小为 230 或更小，则不会。

不过，我不明白为什么 231 大小会导致问题，所以我考虑了混叠，但没有看到任何特别的东西。