Низкая производительность функции, вызывающей Fortran Intel MKL dgemm против numpy, но также и против matmul

Я открываю для себя fortran и ctypes в python, поскольку я планирую предоставить python библиотеку fortran с использованием intel MKL, библиотеку, которая на данный момент вызывается из самой библиотеки c, вызываемой из библиотеки c++... наконец, используемой c#...

Для начала мне удалось написать простой код на Фортране, определяющий различные функции умножения матриц: наивную, одну с использованием fortran matmul и одну с использованием Intel MKL dgemm:

        subroutine matmultorig(M1, M2, M3, M, N, K) bind(c, name='matmultorig')
            !DEC$ ATTRIBUTES DLLEXPORT :: matmultorig
            use iso_c_binding, only: c_double, c_int
        
            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)

            M3 = matmul(M1,M2)
        
        end subroutine

        subroutine matmultmy(M1, M2, M3, M, N, K) bind(c, name='matmultmy')
            !DEC$ ATTRIBUTES DLLEXPORT :: matmultmy
            use iso_c_binding, only: c_double, c_int
        
            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)
        
            integer :: i,j,l
        
            do i = 1,M
                do j = 1,K
                    M3(i,j) = 0.
                    do l = 1,N
                        M3(i,j) = M3(i,j) + M1(i,l) * M2(l,j)
                    end do
                end do
            end do
        
        end subroutine

        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
        
            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)

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

Я компилирую фортран, используя файл .bat (я под окнами):

@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 test.obj test.f
ifort -dll -o mylib.dll test.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, который я выполняю из кода Visual Studio:

from ctypes import *
import time

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 numpy as np

mylib = CDLL(r"C:/path/to/the/fortran/mylib.dll")

mylib.matmultmy.argtypes = [ POINTER(c_double), 
                                POINTER(c_double), 
                                POINTER(c_double),
                                POINTER(c_int),
                                POINTER(c_int),
                                POINTER(c_int) ]

mylib.matmultorig.argtypes = [ POINTER(c_double), 
                                POINTER(c_double), 
                                POINTER(c_double),
                                POINTER(c_int),
                                POINTER(c_int),
                                POINTER(c_int) ]

mylib.matmultmkl.argtypes = [ POINTER(c_double), 
                                POINTER(c_double), 
                                POINTER(c_double),
                                POINTER(c_int),
                                POINTER(c_int),
                                POINTER(c_int) ]

# Setup    
M=500
N=500
K=500

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

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


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

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

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

# Numpy
start = time.time()
c = a.dot(b)
stop = time.time()
print(f"Numpy \t\t {stop - start}s")

Результат:

Fortran my       0.11234903335571289s
Fortran matmul   0.023325443267822266s
Fortran mkl      0.5279343128204346s
Numpy            0.001001596450805664s

Я играл с теми же идеями, но в контексте python С++ pybind11 здесь (те же размеры матрицы):

pybind11 и numpy для матричного продукта

где цифры сейчас несопоставимы с моим случаем, но, по крайней мере, numpy и intel mkl были примерно на одном уровне производительности. Здесь функция, вызывающая dgemm, требует в 500 раз больше матричного произведения numpy. Я подозреваю, что это из-за второстепенной сортировки и в основном из-за привязки c. Однако я обнаружил все это как два дня назад, так что если у профи есть идея ...