Fortran, C言語 との連携
Python からC言語、Fortranを呼び出す¶
C言語やFortranで書かれたライブラリは、直接Pythonから呼び出すことが出来ます。
まず、C言語やFortranをコンパイルして共有ライブラリファイル(*.so (Linux) や *.dll)にして、 そのライブラリをPythonで読み込みます。
!cat libsobol.c
!gcc -cpp -fPIC -shared libsobol.c -lm -o libsobol.so
ライブラリが出来ました。では早速Pythonから呼び出してみましょう。 ライブラリ ctypes と numpy を読み込みます。
from ctypes import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
ついでに、動いていることを確認するためにMatplotlibも読み込みました。
np.ctypeslib.load_library で読み込む共有ライブラリファイルを指定します。
libsobol = np.ctypeslib.load_library("libsobol.so",".")
libsobol.c で定義した関数がすべて読み込まれました。
次に、関数の引数・戻り値の型を教えてあげます。
# C の関数宣言式:
# void sobol(double *points, int N);
# 引数の型を指定
# v 関数の名前
libsobol.sobol.argtypes = [np.ctypeslib.ndpointer(dtype=np.float64), c_int32]
# 戻り値の型を指定
libsobol.sobol.restype = c_void_p
代表的な型は以下のように指定します。
| |型(C言語) |型 (C99 stdint.h) |Python での型指定 | |---|----------------|------------------|---------------------------------------| |倍精度浮動小数点数の配列|double * | |np.ctypeslib.ndpointer(dtype=np.float64)| |単精度浮動小数点数の配列|float * | |np.ctypeslib.ndpointer(dtype=np.float32)| |64bit整数の配列|long long int *|int64_t * |np.ctypeslib.ndpointer(dtype=np.int64) | |32bit整数の配列|int * |int32_t * |np.ctypeslib.ndpointer(dtype=np.int32) | |倍精度浮動小数点数|double | |c_float64| |単精度浮動小数点数|float | |c_float32| |64bit整数|long long int |int64_t |c_int64 | |32bit整数|int |int32_t |c_int32 | |void |void | | c_void_p |
# C言語の関数 sobol に渡す numpy の配列を用意します
a = np.zeros((1024, 2))
# 関数を呼び出します。
libsobol.sobol(a, 1024)
C言語の関数に引数として渡された多次元配列は、C言語からは1次元配列として見えます。(a[0, 0], a[0, 1], a[1, 0], a[1, 1] ...)
では結果を確認してみましょう。
a
plt.scatter(a[:,0], a[:,1])
うまくできました。
Python から Fortran の関数を呼び出す¶
Fortran 2003 の機能により、 Fortran のサブルーチン・関数は C 言語の関数と同様に扱うことができます。Fortran 2003 に対応していないコンパイラでも、場合によっては動かすことができます。
まず、Python から呼び出したいサブルーチン、関数の宣言の最後に bind(C) を書き加えます。
!cat RK4.f90
これを共有ライブラリにコンパイルします。
!gfortran -cpp -fPIC -shared RK4.f90 -o RK4.so
Python から呼び出してみます。Fortran は参照渡しであるため、関数の引数の型を、すべてポインタとして教えていることに注意してください。
libRK = np.ctypeslib.load_library("RK4.so",".")
libRK.rungekutta.argtypes = [POINTER(c_double),
POINTER(c_int32),
np.ctypeslib.ndpointer(dtype=np.float64)]
libRK.rungekutta.restype = c_void_p
a = np.zeros((500,))
# まず、Python の float 型変数 0.01 を c_double 型に変換し、
# つぎにその値の入った場所へのポインタに変換している。
dt = POINTER(c_double)(c_double(0.01))
len = POINTER(c_int32)(c_int(500))
libRK.rungekutta(dt, len, a)
plt.plot(np.arange(0,5,0.01), a)
Fortran 2003 の bind(C) に対応していないコンパイラの場合は、C言語やPythonから呼び出す時に、Fortranで宣言された関数名の末尾にアンダースコア( _ )をつけると呼び出せることが多いです。
f2py を使う¶
scipy に同梱されている f2py を使えば簡単に Python から Fortran を呼び出せます。
6スペース空けてFortranプログラムを書き、入出力の変数を Cf2py intent(in) :: (変数名) で指示します。
!cat RK4_f2py.f90
コンパイルします。-mオプションで出来上がりのモジュール名を指定します。
!f2py -c -m rk RK4_f2py.f90
モジュールができたので、使ってみます。
import rk
intent(in) で指定した変数のみを渡します。
x = rk.rungekutta(0.01, 500)
plt.plot(np.arange(0,5,0.01), x)
うまくできました。