借助C扩展Python模块功能和Numpy的C-API接口实现Python调用Fortran代码

虫儿飞

Python调用Fortran代码，最简单方便的方法是f2py，论坛里就有f2py使用简介，通过在已有的Fortran源代码，添加形如‘!f2py intent(inout/in/out)’的语句再编译成动态库，即可以实现在Python中import。目前f2py是包含在Numpy项目中，但是支持的功能有点少，而且很长时间不更新，所以我要重新找一条技术路线，可以实现Python调用Fortran，目标就是将Fortran代码编译成动态库。
经过调研与初步的技术实践，我选定的技术路线是用C扩展Python模块功能，并借助C&Fortran混编，给Fortran代码编写C接口，最终实现Python调用Fortran。C扩展Python模块功能是Python自带的功能，并且不同于ctypes（ctypes是Python自带的标准库，C扩展Python模块功能是借助Python.h/PyObject实现的一个功能，所以后者具有一定的实现难度），但是只支持简单数据类型作为参数进行传递，如果要支持数组Array功能，就必须借助Numpy的C-API接口。这条技术路线的实现对C、Fortran和Python都有一定要求（特别是Fortran支持默认的位域设定），单纯通过文字将问题讲清楚需要大量的篇幅，而且学语言就是要通过练习来达到融汇贯通，因此通过一个例子来说明。
这个例子是我想实现eemd时剥离出来的。
首先是Fortran代码：
    subroutine array(x,y,n,m)
        implicit none

        integer*4 :: n,m
        real*8 :: x(n),y(m,n)
        integer*8 :: i,j
        write(*,*) n,m
        do i=1,n
            write(*,*) i,x(i)
        enddo
        do i=1,n
            do j=1,m
                y(j,i)=x(i)+j*1.0
            enddo
        enddo
    end

这个代码极其简单，输入参数有四个，一维数组x和二维数组y，以及数组维度n和m，一位数组是输入，二维数组是输出，通过对一维数组加上常数赋值给二维数组。但是有一个问题需要注意：位域需要清晰给定，并且要清楚Python和C中不同类型的位域分别是多少。
然后给出C的接口代码：
#include <Python.h>
#include <numpy/arrayobject.h>
#include <stdio.h>

#define in_N(i) (in_array->dimensions[(i)])
/*#define in_Value(i,j) (*(npy_float64*)((PyArray_DATA(in_array) +                \
                                    (i) * PyArray_STRIDES(in_array)[0] +   \
                                    (j) * PyArray_STRIDES(in_array)[1])))*/
#define in_Value(i) (*(npy_float64*)((PyArray_DATA(in_array)+(i)*PyArray_STRIDES(in_array)[0])))

#define out_N(i) (out_array->dimensions[(i)])
#define out_Value(i,j) (*(npy_float64*)((PyArray_DATA(out_array) +                \
                                    (i) * PyArray_STRIDES(out_array)[0] +   \
                                    (j) * PyArray_STRIDES(out_array)[1])))

static PyObject* eemd(PyObject* self, PyObject* args)
{
    PyArrayObject *in_array;
    PyObject      *out_array;
    int i,j;
    int N;
    int nd = 2;
    int M=10;

    if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &in_array))
        return NULL;

    printf("eemd %i\n",PyArray_NDIM(in_array));

    N = in_N(0);

    for (i=0; i < N; ++i){
        printf("%i %f\n",i,in_Value(i));
    }
    npy_intp dims[2] = {N,M};
    out_array = PyArray_SimpleNew(nd, dims, NPY_FLOAT64);
    if (out_array == NULL)
        return NULL;
    array(PyArray_DATA(in_array),PyArray_DATA(out_array),&N,&M);
    for (i=0; i < N; ++i){
        for(j=0; j < M; ++j){
            printf("%i %i %f %f\n",i+1,j+1,in_Value(i),out_Value(i,j));
        }
    }
    return out_array;
}
C代码部分，首先需要注意的是头文件部分，包含三个Python.h和numpy/arrayobject.h都是实现功能不可缺少的，stdio.h是因为代码中有标准输出。然后是宏定义部分，可以简化代码。接下来就是一个函数，PyObject*对象（有一句话叫做Python的一切皆对象，读一读Python.h就能理解），技术核心。熟悉C语言的应该知道，static PyObject* eemd(PyObject* self, PyObject* args)的括号之前部分表明的是函数名和函数返回类型，表明这个函数返回的是一个Python对象，括号内self不用管（在Python内实现过class都可以理解self），args是参数列表，是从Python接收的参数列表，参数只有一个，是一个一维数组in_array，是一个PyArrayObject，通过PyArg_ParseTuple识别，返回是一个二维数组out_array，是一个PyObject，本质上也是一个PyArrayObject，但是需要注意的是，这些数组与C和Fortran中的数组不是一回事，但是肯定包含有我们需要的数组类型，位域，大小以及数据（这才是C和Fortran的数组）这些信息。细细体会这些C代码，其实也很简单，就是Python传递过来一个一维数组，通过PyArg_ParseTuple识别转换给C语言使用，再创建一个二维数组，将一位数组的数据部分作为输入，二维数组的数据部分作为输出，以及数组维度信息（正好四个参数）作为参数，调用Fortran的subroutine（名称是array），再将输出的二维数组返回Python程序。另外，C&Fortran混编，要注意参数传递是地址传递，另外二维数组的存储方向是相反的（比较一下参数N，M，n，m）。

接下来是封装的Python接口：
static PyMethodDef EemdMethods[] =
{
     {"eemd", eemd, METH_VARARGS,NULL},
     {NULL, NULL, 0, NULL}
};

/* module initialization */
PyMODINIT_FUNC
initEemd(void)
{
     (void) Py_InitModule("Eemd", EemdMethods);
     /* IMPORTANT: this must be called */
     import_array();
}

两个部分，第一部分是函数方法的列表，包含信息有名字，对象，参数以及说明，第二部分是模块初始化部分，init<模块名>，其中import_array()是使用C-API所特有的，否则不能将数组作为参数返回。
接下来还有编译部分：
gcc -c -fpic -fno-underscoring -lpthread CFpy_array.f90
gcc -c -fpic -I/g1/app/apps/python/2.7.14/gnu/include/python2.7 -I/g1/app/apps/python/2.7.14/gnu/lib/python2.7/site-packages/numpy/core/include -lpthread -ldl -lutil -lm -lpython2.7 eemd.c
gcc -o Eemd.so CFpy_array.o eemd.o -shared -lgfortran -lpthread

全部都是gcc编译的（虽然我很想通过intel编译器编译，但是真的没有能够实现），简单说明一下，Fortran代码的编译-fno-underscoring可以实现函数名没有下划线，C代码吧部分，需要有python和numpy的include，链接部分-lgfortran是为了支持Fortran代码。
感兴趣的童鞋评论和私信都可以，网上这些资料太尼玛少了，没有人讨论很痛苦。
C扩展Python模块的小例子，帮助各位理解。

言深深

辛苦，赞！

又是那隻貓

这不如f2py好用呀

931404656

这不如f2py好用呀