前言

本文介绍了如何在Python脚本中,通过引入Fortran动态库调用Fortran子程序,加速Python中的数值计算。以矩阵乘法为例,比较了调用Python中NumPy和SciPy包的方法和使用Fortran子程序两种方式在计算效率上的差别。

Python作为一种动态的解释型编译器,尽管具有新手友好、直观易读、适用范围广的特点,但比起静态的编译型语言Fortran/C,在数值计算方面的表现要差许多,这是由于Python没有被预先编译到机器语言层面。目前流行的Python数值计算包NumPySciPy提供了大量数值计算相关的函数和方法,在很大程度上弥补了Python这一缺点,但仍然无法满足数值计算的全部需求,特别是难以链接一些尚未Python模块化的数学库和工具库。

作为老牌数值计算语言,Fortran拥有许多高效的数学库,我们很容易在Fortran程序中使用他们,但要把他们用到Python中则并不是那么容易。一种解决方案是使用F2PY产生Python接口,它是NumPy项目的一部分。基于F2PY,在Python中调用Fortran函数的基本流程是

  1. 编写使用了数学库的Fortran代码
  2. 在恰当的编译选项下使用f2py编译Fortran代码,产生可供引入的动态库
  3. 在Python中通过import引入动态库

这样就能像调用Python包一样使用Fortran代码中定义的子程序和函数了。下面先简单介绍如何使用f2py,再以矩阵乘法作为例子做说明。

使用F2PY产生供Python引入的动态库

聪明的方法: 创建署名文件

f2py安装就不多说了, 直接上用法.

1
f2py mysubr.F90 -m mysubr -h mysubr.pyf

其中产生的.pyf就是所谓的署名文件(signature file)。其中定义了Python模块mysubr,它包含一个接口,mysubr.F90中所有函数和子程序都被声明在该接口中。每个声明中包含函数所需参量的类型和维度具体信息。具体可参照官方文档。产生.pyf文件后,可用下面的命令

1
f2py -c mysub.pyf mysubr.F90

产生动态库。编译时链接外部函数库(MKL或者FFTW)和一般的编译器相同

1
2
# link mysubr against libabc.a in /lib/dir/
f2py -c mysub.pyf mysubr.F90 -L/lib/dir/ -labc

所有f2py命令可以用python -m numpy.f2py等价替换, 这样做的好处是f2py版本总是与python版本兼容.

举例: 矩阵乘法

矩阵生成

在进行矩阵乘法前,首先利用random模块产生ndarray类型的矩阵,这里取100维的方阵进行测试。具体代码如下所示:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#!/usr/bin/env python3

import numpy as np
from random import seed, random

# set the dimensions
m = 100
k = 100
n = 100

# create mat1(m,k) and mat2(k,n) matrices
seed()
mat1 = np.array([[random() for col in range(k)] for row in range(m)], \
order='F', dtype='float64')
seed()
mat2 = np.array([[random() for col in range(n)] for row in range(k)], \
order='F', dtype='float64')

order='F'使用Fortran的数据存储方式。为保证计算效率,必须手动设置该值

order默认为C。如果不覆盖默认值,Fortran代码的运算效率将下降很多。

执行矩阵乘法

我们比较四种不同的矩阵乘法实现

  1. NumPy的matmul函数

    1
    mat3 = np.matmul(mat1, mat2)
  2. SciPy的linalg.blas.dgemm函数

    1
    2
    from scipy import linalg
    mat3 = linalg.blas.dgemm(1.0, mat1, mat2)
  3. Fortran的matmul函数

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    ! file: f_matmul.F90
    subroutine f_matmul(m, n, k, mat1, mat2, mat3)

    integer, intent(in) :: m, n, k
    real(8), intent(in) :: mat1(m,k), mat2(k,n)
    real(8), intent(out) :: mat3(m,n)

    mat3 = matmul(mat1, mat2)

    end subroutine f_matmul

    m,n,k用于确定矩阵mat1,mat2,mat3的维度,这在纯Fortran中是必须声明的参数,但f2py会将其转化为可选参数,矩阵规模由从Python输入的mat1,mat2确定。在对应Python文件中加入

    1
    from f_matmul import f_matmul

    即可调用.

  4. Intel MKL的dgemm子程序

    将上面Fortran代码中的matmul行替换为

    1
    2
    3
    4
    5
    6
    ! file: f_dgemm.F90
    subroutine f_dgemm(m, n, k, mat1, mat2, mat3)
    !...
    call dgemm('N', 'N', m, n, k, 1.0D0, mat1, m, mat2, k, 0.0D0, mat3, m)

    end subroutine f_dgemm

    此时编译需要链接MKL库,编译命令为

    1
    f2py -c f_dgemm.pyf --fcompiler=intelem --compiler=intelem -L$MKLROOT/lib/intel64/ -lmkl_rt f_dgemm.F90

    在对应Python文件中加入

    1
    from f_dgemm import f_dgemm

所有源码和编译用Makefile打包在这个压缩包里了, 方便取用.

测试结果

测试结果如下

Dimension 3000 (o=C) 5000 (o=C) 3000 5000 10000
np.matmul 1.1611 5.2340 1.1611 5.2404 41.8116
SciPy dgemm 1.3578 5.8640 1.1568 5.2255 41.4386
Fortran matmul 1.4907 6.5730 1.3204 5.9488 48.3035
Intel MKL dgemm 1.2586 5.4831 1.0553 4.8380 38.9050

其中o=C表示使用C存储方式, 不注明则是用Fortran。通过比较可以得到下面的一些结论

  • 在Fortran order下,计算效率顺序为Fortran matmul<np.matmul<scipy.linalg.blas.dgemm< Intel MKL dgemm
  • np.array在内存中的存储方式(order='C'|'F')显著影响Fortran matmul, scipy.linalg.blas.dgemm和Intel MKL dgemm的计算效率,order='F'要比'C'快约10%。
  • 不合理的存储方式('C')导致MKL效率低于np.matmul
  • 存储方式对np.matmul没有显著影响。

总结

本文介绍了使用F2PY工具将Fortran代码编译为可供Python调用的动态库,并以矩阵乘法为例进行演示。同时,比较了流行的Python数值计算模块NumPy和SciPy中的实现与Fortran语言下matmul和MKLdgemm实现的计算效率。测试结果发现,在高维情况下,Intel MKL具有最高的效率。产生ndarray时的存储方式(order=F)对Fortran库中函数计算效率有显著影响。

备注

使用MKL的DGEMM时,用f2py链接MKL库编译后执行test_mm.py,报错

1
Intel MKL FATAL ERROR: Cannot load libmkl_avx2.so or libmkl_def.so

解决方法是在编译前preload几个核心的Intel库

1
export LD_PRELOAD="$MKLROOT/lib/intel64/libmkl_def.so:$MKLROOT/lib/intel64/libmkl_sequential.so:$MKLROOT/lib/intel64/libmkl_core.so"

必须按顺序全部预载入,否则ld会报错。

更新

2019-05-02

  1. 时隔一年各类包都有更新, 原来的代码运行各种问题, 包括在import f_matmul时出现

    1
    ImportError: dynamic module does not define module export function (PyInit_f_matmul)

    这一条pytorch issue陈述了类似的问题, 我的理解是由于f2py和python版本不兼容所致. 用python -m numpy.f2py替代f2py即解决了这一问题.

  2. order='Fortran'替换为order='F'.

  3. 为方便测试, 将所有源码和Makefile打包.

修改后, Fedora 27下运行正常, 但macOS 10.14.4下会报错误

1
ImportError: dlopen(f_matmul.cpython-37m-darwin.so, 2): __dyld section not supported in f_matmul.cpython-37m-darwin.so

根据这条18年10月的Intel Forum帖, 这个问题似乎跟Xcode版本有关. 因为要赶ddl所以暂时没有考虑macOS, 如果有人知道如何解决的话麻烦指点我一下 :)

Comments