用f2py在Python中调用Fortran矩阵乘法会变慢

Question

我一直在尝试使用f2py将优化的Fortran代码与Python接口化，用于向量和矩阵乘法。 为了获得对比性能的数据，我在循环中执行相同的乘积100000次。 使用完整的Fortran代码，乘积需要2.4秒（ifort），而使用f2py需要大约11秒。仅供参考，使用numpy需要约20秒。 我要求Fortran和Python部分在循环前后写入时间差异，使用f2py时它们都写入11秒，因此代码没有在传递数组时浪费时间。我试图理解是否是numpy数组存储方式引起的问题，但我无法理解问题所在。 你有任何想法吗？ 提前致谢

program Main
    implicit none
    save

    integer :: seed, i, j, k
    integer, parameter :: states =15
    integer, parameter :: tessere = 400
    real, dimension(tessere,states,states) :: matrix
    real, dimension(states) :: vector
    real :: start, finish
    real  :: prod(tessere)

    do i=1,tessere
       do j=1,states
          do k=1,states
              matrix(i,j,k) = i+j+k
          end do
       enddo
    end do
    do i=1,states
        vector(i) = i
    enddo
    call doubleSum(vector,vector,matrix,states,tessere,prod)

end program

Fortran子程序：

subroutine doubleSum(ket, bra, M , states, tessere,prod)
    integer :: its, j, k,t
    integer :: states
    integer :: tessere
    real, dimension(tessere,states,states) :: M
    real, dimension(states) :: ket
    real, dimension(states) :: bra
    real, dimension(tessere) :: prod
    real,dimension(tessere,states) :: ctmp

    call cpu_time(start)
    do t=1,100000
        ctmp=0.d0
        do k=1,states
             do j=1,states
                do its=1,tessere
                   ctmp(its,k)=ctmp(its,k)+ M(its,k,j)*ket(j)
                enddo
             enddo
        enddo
        do its=1,tessere
            prod(its)=dot_product(bra,ctmp(its,:))
        enddo
    enddo
    call cpu_time(finish)
    print '("Time = ",f6.3," seconds.")',finish-start
end subroutine

Python脚本

import numpy as np
import time
import cicloS


M= np.random.rand(400,15,15)
ket=np.random.rand(15)

#M=np.asfortranarray(M)
#ket=np.asfortranarray(ket)

import time


start=time.time()  
prod=cicloS.doublesum(ket,ket,M)
end=time.time()
print(end-start)

使用f2py生成并编辑的.pyf文件

!    -*- f90 -*-
! Note: the context of this file is case sensitive.

python module cicloS 
    interface  
        subroutine doublesum(ket,bra,m,states,tessere,prod) 
            real dimension(states),intent(in) :: ket
            real dimension(states),depend(states),intent(in) :: bra
            real dimension(tessere,states,states),depend(states,states),intent(in) :: m
            integer, optional,check(len(ket)>=states),depend(ket) :: states=len(ket)
            integer, optional,check(shape(m,0)==tessere),depend(m) :: tessere=shape(m,0)
            real dimension(tessere),intent(out) :: prod
        end subroutine doublesum
    end interface
end python module cicloS

- MartaR

---

您提供的Fortran子程序似乎存在错误，integer，dimension（tessere）:: prod 应该是 real，dimension（tessere）:: prodotto，并且子程序中所有出现的 prod 都应替换为 prodotto。您使用什么编译标志来编译ifort版本的代码，并且是否对f2py使用相同的标志？（使用 gfortran -O3 sub.f90 main.f90 进行编译，我得到了与纯Fortran相同的运行时间，与f2py版本大约为1.3秒） - jbdv

Fortran例程中的错误是问题中的拼写错误，我已经进行了更正。实际上，在检查编译标志时，我发现在我的实验室集群中调用f2py的标准方式是gfortran -O3，而如果加载ifort模块，则变为ifort -O1，因此这就是性能差异的原因。您是否知道是否可以告诉f2py使用哪个编译器和标志？ - MartaR

我本来想说“错误”，但更应该说“笔误”或“不一致” :) 不同的编译器标志（以及 ifort vs. gfortran！）肯定可以解释您观察到的时间差异。 您可以通过在f2py命令中添加例如 --opt='-O1' 来指定（附加）编译器标志，并搜索 --help-fcompiler 以了解如何指定f2py应使用哪个Fortran编译器。 - jbdv

f2py标志完美地工作！我确定它们存在，但我无法在网上找到它们！非常感谢您的帮助！ - MartaR

我刚刚发布了一个回答，总结了一切（并附有适当文档的链接），以备将来参考。 - jbdv

Answer 1

该OP指出，独立版本和F2PY编译版本的代码执行时间差异是由于使用了不同的编译器和编译器标志。

为了获得一致的结果，并回答问题，有必要确保F2PY使用所需的1）编译器和2）编译器标志。

第一部分：指定F2PY应使用哪个Fortran编译器

可以通过执行例如python -m numpy.f2py -c --help-fcompiler来显示目标系统上F2PY可用的Fortran编译器列表。在我的系统上，这将产生以下结果（缩短）：

Fortran compilers found:
  --fcompiler=gnu95    GNU Fortran 95 compiler (7)
  --fcompiler=intelem  Intel Fortran Compiler for 64-bit apps (19.0.1.144)

您可以通过在编译命令中添加相应的 --fcompiler 标志来指示 F2PY 使用可用的 Fortran 编译器。例如，要使用 ifort（假设您已经创建并编辑了一个 cicloS.pyf 文件）：

python -m numpy.f2py --fcompiler=intelem -c cicloS.pyf sub.f90

第二部分：指定供F2PY使用的附加编译器标志

请注意，在上一步中，--help-fcompiler输出还显示了F2PY为每种可用编译器定义的默认编译器标志（例如compiler_f90）。在我的系统上，这是（裁剪并简化为最相关的标志）：

• gnu95: -O3 -funroll-loops
• intelem: -O3 -xSSE4.2 -axCORE-AVX2,COMMON-AVX512

您可以通过在编译命令中使用--opt标志（也可参见文档中的--f90flags）来指定F2PY所需的首选优化标志，示例如下：

python -m numpy.f2py --fcompiler=intelem --opt='-O1' -c cicloS.pyf sub.f90

比较独立和F2PY版本的运行时间：

使用ifort -O1 sub.f90 main.f90 -o main编译独立可执行文件和从第二部分编译的F2PY版本，我得到以下输出：

./main
Time =  5.359 seconds.

python test.py
Time =  5.297 seconds.
5.316878795623779

然后，使用ifort -O3 sub.f90 main.f90 -o main编译一个独立的可执行文件，以及从第一部分得到的（默认的）F2PY编译版本，我得到了以下结果：

./main
Time =  1.297 seconds.

python test.py
Time =  1.219 seconds.
1.209657907485962

因此，独立版本和F2PY版本显示类似的结果，并且编译器标志对结果有影响。

关于临时数组的注释

尽管不是您观察到的减速的原因，但请注意，在您的Python示例中，F2PY被迫复制数组M（和ket）的临时副本，原因如下：

• 您传递给的3D数组M是默认的NumPy数组，具有C排序（行主要）。由于Fortran使用列主要排序，因此F2PY将制作数组副本。注释掉的应该解决此问题的一部分。
• 用于数组副本的下一个原因（也适用于ket）是Python（默认64位，双精度浮点数）和Fortran（real提供默认精度，可能为32位浮点数）方面的实际类型不匹配。因此，再次进行复制以解决此问题。

您可以通过在F2PY编译命令中添加-DF2PY_REPORT_ON_ARRAY_COPY = 1 标志（也在文档中）来获得有关数组副本的通知。在您的情况下，可以通过在Python中更改M和ket矩阵的dtype（即，M = np.asfortranarray（M，dtype = np.float32）和ket = np.asfortranarray（ket，dtype = np.float32）），或者通过在Fortran代码中使用适当的kind定义real变量（例如，在您的子程序和主程序中添加use，intrinsic :: iso_fortran_env，only：real64并使用real（kind = real64）定义实数）。

- jbdv