Как я могу улучшить производительность этого кода Fortran, взаимодействующего с Python?

Я изучаю возможность интерфейса кода Fortran для использования в Python. Я знаю о f2py, но, поскольку мне не удалось использовать его с внешними библиотеками (такими как lapack), я вернулся к использованию ctypes. У меня все работает, но я заметил такое поведение: для одного и того же алгоритма (простое поэлементное умножение двух матриц) код на Фортране медленнее numpy, а я ожидал, что будет быстрее. Я пробовал с двумя разными ПК, один с Windows 10, а другой с Windows 11. Я использовал gfortran из MinGW-64, установленный через winlibs.com версии gcc-12.2.0-llvm-15.0.7-mingw-w64ucrt-10.0.0-r4, моя версия Python 3.9.13

Измеренные времена со временем это (для матриц 1000x1000)

С Numpy: 2,77 мс ± 110 мкс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, по 100 циклов в каждом)

С Fortran: 12,4 мс ± 314 мкс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, по 100 циклов в каждом)

Код Fortran находится в файле "libfortran3.f90"

subroutine prodotto(a, b, n) bind(c, name='prodotto')
  ! declarations
  use iso_c_binding, only: c_double, c_int

  integer(c_int), intent(in) :: n
  real(c_double), dimension(n,n), intent(in) :: b
  real(c_double), dimension(n,n), intent(inout) :: a
  
  do i = 1,n
    do j = 1,n
        a(i,j) = a(i,j)*b(i,j)
    end do
  end do
  
end subroutine prodotto

и он скомпилирован с помощью gfortran -O3 -funroll-loops -ffast-math -fPIC -shared -o libfortran3.so libfortran3.f90. Обратите внимание, что я использовал одну и ту же матрицу для ввода и вывода, поэтому мне не нужно передавать другую матрицу для результата (так как я ожидаю, что это усугубит проблему).

Код Python (запущенный в блокноте)

from ctypes import CDLL, POINTER, c_int, c_double
import numpy as np
import time
import sys

fortran = CDLL('./libfortran3.so')
fortran.prodotto.argtypes = [ POINTER(c_double), 
                                 POINTER(c_double), 
                                 POINTER(c_int)]
fortran.prodotto.restype = None
#
N = 10
A = np.random.rand(N,N)
B = np.random.rand(N,N)
#
print('With Numpy')
%timeit A*B
#
print('With Fortran')
A = np.asfortranarray(A, dtype=c_double)
B = np.asfortranarray(B, dtype=c_double)
Act = A.ctypes.data_as(POINTER(c_double))
Bct = B.ctypes.data_as(POINTER(c_double))
%timeit fortran.prodotto( Act, Bct, c_int(N) )

Конечно, мой код Fortran не так оптимизирован, как Numpy, но я хотел бы знать, есть ли что-то, что я мог бы изменить, чтобы добиться лучшей производительности.

Редактировать 3 апреля 2023 г.

Благодаря mkrieger1 и janneb, я действительно принял неправильное решение во вложенном цикле! Переключение i с j в "libfortran3.f90" дает мне это время (с матрицами 1000x1000):

С Numpy: 2,67 мс ± 148 мкс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, по 100 циклов в каждом)

С Fortran: 1,01 мс ± 34,9 мкс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1000 циклов в каждом)

Я также пробовал с f2py со следующим кодом в файле «libfortran1.f90»

subroutine prodotto(a, b, c, n) 
  ! declarations

  integer, intent(in) :: n
  real(kind=8), dimension(n,n), intent(in) :: a, b
  real(kind=8), dimension(n,n), intent(out) :: c
  
  do j = 1,n
    do i = 1,n
        c(i,j) = a(i,j)*b(i,j)
    end do
  end do
  
end subroutine prodotto

скомпилировано с помощью python -m numpy.f2py --compiler=mingw32 --fcompiler=gfortran -c libfortran1.f90 -m libfortran1, запуск его в блокноте дает мне:

С f2py: 12,3 мс ± 260 мкс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 100 циклов в каждом)

что намного медленнее, чем чистый numpy.

Просто для полноты, если я попытаюсь скомпилировать этот код на фортране, который использует Lapack (который я создал в своей системе) в «libfortran0.f90»,

subroutine risolvi(A,b,x,n) bind(c, name='risolvi')
  ! declarations
  use iso_c_binding, only: c_double, c_int

  integer(c_int), intent(in) :: n
  real(c_double), dimension(n,n), intent(inout) :: A
  real(c_double), dimension(n), intent(inout) :: b
  real(c_double), dimension(n), intent(out) :: x
  integer(c_int), dimension(n) :: pivot(n)
  integer(c_int) :: ok
  ! call LAPACK subroutine
  call DGESV(n, 1, A, n, pivot, b, n, ok) 
  x = b
  
end subroutine risolvi

скомпилировано с помощью gfortran -O3 -funroll-loops -ffast-math -fPIC -shared -L. -lblas -llapack -o libfortran0.so libfortran0.f90, и я сравниваю с numpy

from ctypes import CDLL, POINTER, c_int, c_double
import numpy as np
import time
import sys
import libfortran1

fortran = CDLL('./libfortran0.so')
fortran.risolvi.argtypes = [ POINTER(c_double), 
                                 POINTER(c_double), 
                                 POINTER(c_double), 
                                 POINTER(c_int)]
fortran.risolvi.restype = None
#
N = 1000
A = np.random.rand(N,N)
b = np.random.rand(N)
#
print('With Numpy')
%timeit np.linalg.solve(A,b)
#
print('With Fortran')
A = np.asfortranarray(A, dtype=c_double)
b = np.asfortranarray(b, dtype=c_double)
x = np.asfortranarray(np.zeros(N), dtype=c_double)
Act = A.ctypes.data_as(POINTER(c_double))
bct = b.ctypes.data_as(POINTER(c_double))
xct = x.ctypes.data_as(POINTER(c_double))
fortran.risolvi( Act, bct, xct, c_int(N) )
%timeit fortran.risolvi( Act, bct, xct, c_int(N) )

Я получаю следующие сроки (с N = 1000)

С Numpy: 13,6 мс ± 596 мкс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, по 100 циклов в каждом)

С Fortran: 221 мс ± 9,74 мс на цикл (среднее значение ± стандартное отклонение для 7 запусков, по 1 циклу в каждом)

В любом случае, при N=30 использование Фортрана происходит быстрее.

С Numpy: 15,6 мкс ± 244 нс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, 100000 циклов в каждом)

С Fortran: 11,4 мкс ± 237 нс на цикл (среднее значение ± стандартное отклонение для 7 прогонов, 100000 циклов в каждом)

Возможно ли, что в моем коде что-то не так при передаче больших массивов? Или просто, как я читал, "ctypes сделан для совместимости, а не для скорости"?