Python与Fortran混合编程

简介：Fortran程序在地球物理学中并不少见，90年代Fortran十分流行，386时代计算机时很珍贵，Fortran作为当时的高级语言，具有运行速度快，语法简单等诸多优点。现在虽然Fortran已经成为小众游戏，但过去曾经的那些基础计算程序还是有价值的，如果全要重新改写，不但耗时耗力还容易出错。因此，今天我们讨论一下如何在Python里面调用Fortran的函数。

1、认识f2py

混合编程方法有很多，但是f2py可以作为众多方法的首选。f2py是NumPy的一部分，不需要独立安装配置。f2py简单来说就是 Fortran to Python的接口生成器。使用它可以实现Python和Fortran语言之间提供的连接。

在安装了NumPy之后，它还可以作为独立的命令行工具f2py使用，它可以方便地创建/构建Python C/API扩展模块。支持调用Fortran 77/90/95外部子程序和Fortran 90/95模块子程序以及C函数; 支持访问Fortran 77通用块和Fortran 90/95模块数据，包括可分配数组。

f2py通过创建一个可在Python中使用import关键字导入的扩展模块来工作。该模块包含可以从Python调用的自动生成的包装器函数，充当Python和编译后的Fortran例程之间的接口。

编译后的fortran程序依据平台的不同，在磁盘上会产生一个以pyd为扩展名的文件，可以在python程序中被import。

MATLAB和Python这样的脚本语言虽然容易上手，但是对于大规模的循环和数值计算速度还是逊色于C/ c++和Fortran这样的编译语言。要实现两全其美，混合编程是一个好的选择，而Numpy这种大型的python包，其底层也是要通过混合编程来加速的。

2、Fortran的函数接口设计

我们从一个简单的例子开始，Fortran代码保存在fortran.f95中。一般作为接口函数的代码都放到一个subroutine中，这里面在fortran代码里面通过intent关键字来定义输入/输出变量和类型。

下面是fortran代码中的subroutine的典型写法：

subroutine logical_to_integer(prime_numbers, is_prime, num_primes, n)
! =====================================================
! Translates the logical array from sieve to an array
! of size num_primes of prime numbers.
! =====================================================
    integer                 :: i, j=0
    integer, intent(in)     :: n
    logical, intent(in)     :: is_prime(n)
    integer, intent(in)     :: num_primes
    integer, intent(out)    :: prime_numbers(num_primes)
    do i = 1, size(is_prime)
        if (is_prime(i)) then
            j = j + 1
            prime_numbers(j) = i
        end if
    end do
end subroutine

需要注意的是在接口函数中不要使用动态数组，与python交换数据一定要指定数组大小。

3、Python的函数调用方法

Python要调用fortran的函数，必须对fortran函数进行编译。
可以在虚拟环境中，使用下面的命令：

f2py -c fortran.f95 -m fortran

或者：

python -m numpy.f2py -c fortran.f95 -m fortran

在Python中，可以直接import上面的函数名

import fortran
print(fortran.__doc__)
print(fortran.logical_to_integer.__doc__)

注意上面的__doc __是f2py自动生成的，可以看到fortran模块里面包含几个函数，每个函数里面还可以再调用doc看到接口参数类型。

在实践中更常见的做法是在python中在设计一个与fortran函数同名的wrapper函数，并设计新的docstring，这样用起来更像一个python函数。

4、Setup.py中的编译

如果您在程序包里面用到了混合编程方法，还需要编写setup.py，这样软件包安装后，自动会完成fortran代码的编译。方法如下：

# Fortran extension
from numpy.distutils.core import Extension 
from numpy.distutils.core import setup 
ext_for1 = Extension(name='geoist.name.for1',
                sources=['geoist/name/for1.f90'])
setup(
  ext_modules         = [ext_for1],
)

上面代码在setup.py中添加后，当geoist这个软件包安装时，for1.f90就会被自动的根据平台类型进行编译。

5、什么时候使用混合编程

如果想问什么时候要用到混合编程，我想这永远不会有没有明确的答案。然而，一个好的经验法则是，当在(嵌套的)循环中执行多个操作/-计算时，可以使用f2py或编译语言。最典型的例子可能是对多维矩阵中的元素进行操作。也就是一般的线性代数。其他好的例子可以是计算积分的程序或进行蒙特卡罗模拟。此时，您可能想知道是否已经有人为您的特定问题制作了免费的模块。

答案很可能是肯定的!

NumPy和SciPy中的大多数函数和例程实际上都是用Fortran(或C/ c++)例程编译的，它们为解决多个问题提供了高效和快速的解决方案。

一句话总结：如果您是从fortran阵营转到python来的，并且希望使用更多python生态工具，或者您期望通过fortran来加速python程序运行速度，我想这篇文章是有帮助的。对于一些经典的数值计算算法，我们建议您经常检查这Numpy和Scipy这两个包中是否已经提供了解决方案。如果没有，您可以尝试通过混合编程的方法，继承原来的程序框架，而f2py可能会为您的问题提供最佳解决方案。