Python による Fortran の手続きの呼び出し

2016年8月31日

はじめに

Python による Fortran の手続きの呼び出しについて。

環境

• Windows 7 64 bit/MSYS2 64 bit (MinGW w64)
• Ubuntu Linux 12.04 LTS
• Anaconda3 4.1.1

方針

Fortran の手続きを呼び出すため、C でラッパー関数を書いて、CFFI で呼び出す。

コード

サンプルコードは以下の通り。

test.f

      subroutine test(n, a)
        implicit none
        integer :: n
        double precision :: dx
        double precision, parameter::pi=3.1415926535897932
        double precision, dimension(2, n) :: a
        double precision :: x, f
        integer i

        dx = 2*pi/(n-1)
        do i = 1, n
          x = i*dx
          f = sin(x)
          a(1, i) = x
          a(2, i) = f
        enddo
      end

test_f.h

#ifndef TEST_H
#define TEST_H

void test_(int *n, double a[][2]);

#endif /* TEST_H */

test.c

#include "test_f.h"

void test(int n, double a[][2])
{
	test_(&n, a);
}

Makefile (Windows)

all:
	gfortran -Wall -O2 -o test_f.o -c test.f
	gcc -Wall -O2 -c test.c
	gcc -shared -o test.so test_f.o test.o

clean:
	rm -f *.o

Makefile (Linux)

all:
	gfortran -fPIC -Wall -O2 -o test_f.o -c test.f
	gcc -fPIC -Wall -O2 -c test.c
	gcc -shared -o test.so test_f.o test.o

clean:
	rm -f *.o

x = [0, 2π] の範囲で sin 関数を計算した結果をリストで返すものである。

注意点

• Fortran の手続きの名前には、最後にアンダースコアがつく。
• 引数はすべて参照渡しである。
• 配列の添字は、C は 0 から始まるが Fortran は 1 から始まる。
• 多次元配列のメモリ内での並び方が C と Frotran で異なる。たとえば 2x3 の配列について、C では (0, 0)、(0, 1)、(1, 0)、(1, 1) ...　となっているが、Fortran では (1, 1)、(2, 1)、(1, 2)、(2, 1)... という順番である。

Python から C の関数を呼び出す

Python から CFFI を用いて C の関数を呼び出す。

call_c_func_cffi.py

from matplotlib import pyplot as plt
import numpy as np
from cffi import FFI


ffi = FFI()

ffi.cdef(r'''
void test(int n, double a[][2]);
''')

lib = ffi.dlopen('./test.so')

def test(n):
    _r = ffi.new('double[%d][2]' % n)

    lib.test(n, _r)

    r = np.array([[elem[i] for i in (0, 1)] for elem in _r])

    return r

r = test(20)
for i in range(20):
    print('%g %g' % (r[i, 0], r[i, 1]))

plt.plot(r[:, 0], r[:, 1])
plt.show()

Python から Fortran の手続きを呼び出す

実は CFFI で直接 Fortran の手続きを呼び出すこともできる (Fortran のコードを C でラッピングできているのだから、当然)。

Makefile を次のように書き換える。

Makefile (Windows)

all:
	gfortran -Wall -O2 -o test_f.o -c test.f
	gcc -shared -o test.so test_f.o

clean:
	rm -f *.o

Makefile (Linux)

all:
	gfortran -fPIC -Wall -O2 -o test_f.o -c test.f
	gcc -shared -o test.so test_f.o

clean:
	rm -f *.o

Fortran の手続きを呼び出すコードは以下の通り。

call_f_func_cffi.py

from matplotlib import pyplot as plt
import numpy as np
from cffi import FFI


ffi = FFI()

ffi.cdef(r'''
void test_(int *n, double a[][2]);
''')

lib = ffi.dlopen('./test.so')

def test(n):
    _n = ffi.new('int[]', [n])
    _r = ffi.new('double[%d][2]' % n)

    lib.test_(_n, _r)

    r = np.array([[elem[i] for i in (0, 1)] for elem in _r])

    return r

r = test(20)
for i in range(20):
    print('%g %g' % (r[i, 0], r[i, 1]))

plt.plot(r[:, 0], r[:, 1])
plt.show()

関数名にアンダースコアを付ける必要があること、引数をすべて配列の形にしなければならないところが少しややこしいかもしれない。