如何提升Fortran与Python接口代码的性能？

背景与初始问题 #

我正在研究Fortran代码与Python对接的可行性。了解f2py工具，但因无法成功调用LAPACK等外部库，转而使用ctypes实现对接。目前功能正常，但发现相同的矩阵逐元素乘法算法，Fortran代码运行速度慢于NumPy，这与预期不符。测试在Windows 10和Windows 11两台电脑上进行，使用winlibs.com提供的gcc-12.2.0-llvm-15.0.7-mingw-w64ucrt-10.0.0-r4版本gfortran，Python版本为3.9.13。

初始性能测试结果（1000×1000矩阵） #

• NumPy：2.77 ms ± 110 µs per loop（7次运行，每次100循环的均值±标准差）
• Fortran：12.4 ms ± 314 µs per loop（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调用代码（Notebook中运行） #

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，但想了解是否有可修改之处以提升性能。

2023年4月3日更新 #

根据建议调整嵌套循环顺序后，性能得到大幅提升！将libfortran3.f90中的i与j交换后，1000×1000矩阵的测试结果如下：

• NumPy：2.67 ms ± 148 µs per loop（7次运行，每次100循环的均值±标准差）
• Fortran：1.01 ms ± 34.9 µs per loop（7次运行，每次1000循环的均值±标准差）

f2py测试尝试 #

我还尝试了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

Notebook中运行结果：

• f2py：12.3 ms ± 260 µs per loop（7次运行，每次100循环的均值±标准差）
  性能远低于纯NumPy。

LAPACK调用性能测试补充 #

我尝试编译调用LAPACK的Fortran代码（文件libfortran0.f90，LAPACK为自行编译）：

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

Python调用代码 #

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 ms ± 596 µs per loop（7次运行，每次100循环的均值±标准差）
  • Fortran：221 ms ± 9.74 ms per loop（7次运行，每次1循环的均值±标准差）
• N=30时，Fortran更快：
  • NumPy：15.6 µs ± 244 ns per loop（7次运行，每次100000循环的均值±标准差）
  • Fortran：11.4 µs ± 237 ns per loop（7次运行，每次100000循环的均值±标准差）

想确认是大数组传递代码存在问题，还是ctypes本身侧重兼容性而非性能？

内容的提问来源于stack exchange，提问作者Federico Cluni