首页 > 解决方案 > 当由 python 使用 f2py 调用时,fortran 矩阵乘积变慢

问题描述

我一直在尝试使用 f2py 将优化的 fortran 代码用于向量和矩阵乘法与 python 的接口。为了获得对我的目的有用的性能比较,我在一个周期内执行相同的产品 100000 次。使用完整的 fortran 代码,该产品需要 2.4 秒(ifort),而使用 f2py 大约需要 11 秒。仅供参考,使用 numpy 大约需要 20 秒。我要求 fortran 和 python 部分都写出循环前后的时间差,使用 f2py 他们都写 11 秒,所以代码在传递数组时不会浪费时间。我试图了解它是否是存储 numpy 数组的方式,但我无法理解问题所在。你有什么主意吗?提前致谢

fortran 主要

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

蟒蛇脚本

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

标签: pythonf2py

解决方案


OP 表示,在代码的独立编译版本和 F2PY 编译版本之间观察到的执行时间差异是由于使用了不同的编译器和编译器标志。

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

第 1 部分:指定F2PY 应使用哪个Fortran 编译器

通过执行例如,可以显示目标系统上 F2PY 可用的 Fortran 编译器列表python -m numpy.f2py -c --help-fcompiler。在我的系统上,这会产生(截断):

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

第 2 部分:指定F2PY 使用的其他编译器标志

请注意,上一步的输出还显示了F2PY 为每个可用编译器定义--help-fcompiler的默认编译器标志(参见例如)。compiler_f90再次在我的系统上,这是(截断并简化为最相关的标志):

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

您可以使用编译命令中的标志指定 F2PY 的首选优化标志--opt(另请参见文档),现在变为例如:--f90flags

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

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

使用 编译一个独立的可执行文件,以及来自第 2 部分ifort -O1 sub.f90 main.f90 -o main的 F2PY 编译版本,我得到以下输出:

./main
Time =  5.359 seconds.

python test.py
Time =  5.297 seconds.
5.316878795623779

然后,使用第 1 部分ifort -O3 sub.f90 main.f90 -o main中的(默认)F2PY 编译版本编译一个独立的可执行文件,我得到以下结果:

./main
Time =  1.297 seconds.

python test.py
Time =  1.219 seconds.
1.209657907485962

因此显示了独立版本和 F2PY 版本的类似结果,以及编译器标志的影响。

评论临时数组

虽然不是您观察到的减速的原因,但请注意 F2PY 在您的 Python 示例中被迫制作数组M(和ket)的临时副本,原因有两个:

  • M您传递给的 3D 数组cicloS.doublesum()是默认的 NumPy 数组,具有 C 排序(行优先)。由于 Fortran 使用列优先排序,F2PY 将制作数组副本。注释掉的np.asfortranarray()应该纠正这部分问题。
  • 数组副本(也用于 )的下一个原因是示例的 Python(默认 64 位,双精度浮点数)和 Fortran(给出默认精度,可能是 32 位浮点数)ket上的真实类型之间存在不匹配。real因此,再次制作副本以说明这一点。

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


推荐阅读