トップページに戻る

Fortran, C言語 との連携

Python からC言語、Fortranを呼び出す

C言語やFortranで書かれたライブラリは、直接Pythonから呼び出すことが出来ます。

まず、C言語やFortranをコンパイルして共有ライブラリファイル(*.so (Linux) や *.dll)にして、 そのライブラリをPythonで読み込みます。

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

まず、呼び出される側のC言語のライブラリを用意します。

ここではモンテカルロ法のための乱数の生成ライブラリ libsobol を使ってみましょう。

In [1]:
!cat libsobol.c
#include <math.h>
#include <stdint.h>
#include <stdlib.h>
int rightmost_unset_bit(int i) {
  return log2(~i & -~i) + 1; // rightmost unset bit
}
void sobol(double *points, int N) {
  int L = (unsigned)ceil(log2((double)N));
  uint64_t V[64];
  for (int j = 0; j < 2; j++) {
    for (int i = 0; i < L; i++) {
      uint64_t m = 0;
      for (int k = 0; k < 64; k += 30) {
        m = m * (RAND_MAX + 1Lu) + rand(); // 64 bit random integer
      }
      V[i] = m << (uint64_t) (63-i);
    }
    uint64_t X = 0;
    for (int i = 0; i < N; i++) {
      X = X ^ V[rightmost_unset_bit(i)];
      points[2*i+j] = X / 0x1.0p64; /* floating-point hexadecimal literal 2**64 */
    }
  }
}
In [2]:
!gcc -cpp -fPIC -shared libsobol.c -lm -o libsobol.so

ライブラリが出来ました。では早速Pythonから呼び出してみましょう。 ライブラリ ctypesnumpy を読み込みます。

In [3]:
from ctypes import *

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline 

ついでに、動いていることを確認するためにMatplotlibも読み込みました。

np.ctypeslib.load_library で読み込む共有ライブラリファイルを指定します。

In [4]:
libsobol = np.ctypeslib.load_library("libsobol.so",".")

libsobol.c で定義した関数がすべて読み込まれました。

次に、関数の引数・戻り値の型を教えてあげます。

In [5]:
# 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
In [6]:
# C言語の関数 sobol に渡す numpy の配列を用意します
a = np.zeros((1024, 2))

# 関数を呼び出します。
libsobol.sobol(a, 1024)
Out[6]:
1024

C言語の関数に引数として渡された多次元配列は、C言語からは1次元配列として見えます。(a[0, 0], a[0, 1], a[1, 0], a[1, 1] ...)

では結果を確認してみましょう。

In [7]:
a
Out[7]:
array([[  7.50000000e-01,   7.50000000e-01],
       [  2.50000000e-01,   0.00000000e+00],
       [  5.00000000e-01,   7.50000000e-01],
       ..., 
       [  7.50007613e-01,   7.50007613e-01],
       [  7.61304574e-06,   7.61304574e-06],
       [  1.67316778e-08,   1.67316778e-08]])
In [8]:
plt.scatter(a[:,0], a[:,1])
Out[8]:
<matplotlib.collections.PathCollection at 0x7fb9a7d29e10>

うまくできました。

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

Fortran 2003 の機能により、 Fortran のサブルーチン・関数は C 言語の関数と同様に扱うことができます。Fortran 2003 に対応していないコンパイラでも、場合によっては動かすことができます。

まず、Python から呼び出したいサブルーチン、関数の宣言の最後に bind(C) を書き加えます。

In [9]:
!cat RK4.f90
subroutine rungekutta(dt, nend, values) bind(C)
  implicit none
  double precision, intent(in) :: dt
  integer, intent(in)  :: nend
  double precision, intent(out) :: values(nend)
  double precision :: t, p1, p2, p3, p4, u
  integer :: n
  u = 1d0
  do n = 1, nend
    t = dt * real(n)
    p1 = u
    p2 = u + 0.5d0 * dt * p1
    p3 = u + 0.5d0 * dt * p2
    p4 = u + dt * p3
    u = u + dt * 1/6d0 * (p1 + 2 * p2 + 2 * p3 + p4)
    values(n) = u
  end do
end subroutine

これを共有ライブラリにコンパイルします。

In [10]:
!gfortran -cpp -fPIC -shared RK4.f90 -o RK4.so

Python から呼び出してみます。Fortran は参照渡しであるため、関数の引数の型を、すべてポインタとして教えていることに注意してください。

In [11]:
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)
Out[11]:
[<matplotlib.lines.Line2D at 0x7fb9a7ac65f8>]

Fortran 2003 の bind(C) に対応していないコンパイラの場合は、C言語やPythonから呼び出す時に、Fortranで宣言された関数名の末尾にアンダースコア( _ )をつけると呼び出せることが多いです。

f2py を使う

scipy に同梱されている f2py を使えば簡単に Python から Fortran を呼び出せます。

6スペース空けてFortranプログラムを書き、入出力の変数を Cf2py intent(in) :: (変数名) で指示します。

In [12]:
!cat RK4_f2py.f90
      subroutine rungekutta(dt, nend, ret)
        implicit none
        double precision :: dt
        integer :: nend
        double precision :: ret(nend)
Cf2py intent(in) :: dt
Cf2py intent(in)  :: nend
Cf2py intent(out) :: ret
      double precision :: t, p1, p2, p3, p4, u
      integer :: n
      u = 1d0
      do n = 1, nend
        t = dt * real(n)
        p1 = u
        p2 = u + 0.5d0 * dt * p1
        p3 = u + 0.5d0 * dt * p2
        p4 = u + dt * p3
        u = u + dt * 1/6d0 * (p1 + 2 * p2 + 2 * p3 + p4)
        ret(n) = u
      end do
      end subroutine

コンパイルします。-mオプションで出来上がりのモジュール名を指定します。

In [13]:
!f2py -c -m rk RK4_f2py.f90
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "rk" sources
f2py options: []
f2py:> /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.c
creating /tmp/tmpf61pdgjh/src.linux-x86_64-3.6
Reading fortran codes...
	Reading file 'RK4_f2py.f90' (format:fix)
Post-processing...
	Block: rk
			Block: rungekutta
Post-processing (stage 2)...
Building modules...
	Building module "rk"...
		Constructing wrapper function "rungekutta"...
		  ret = rungekutta(dt,nend)
	Wrote C/API module "rk" to file "/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.c"
  adding '/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.c' to sources.
  adding '/tmp/tmpf61pdgjh/src.linux-x86_64-3.6' to include_dirs.
copying /usr/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmpf61pdgjh/src.linux-x86_64-3.6
copying /usr/lib/python3.6/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmpf61pdgjh/src.linux-x86_64-3.6
build_src: building npy-pkg config files
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
customize Gnu95FCompiler
Found executable /usr/bin/gfortran
customize Gnu95FCompiler
customize Gnu95FCompiler using build_ext
building 'rk' extension
compiling C sources
C compiler: gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -march=x86-64 -mtune=generic -O2 -pipe -fstack-protector-strong -march=x86-64 -mtune=generic -O2 -pipe -fstack-protector-strong -march=x86-64 -mtune=generic -O2 -pipe -fstack-protector-strong -fPIC

creating /tmp/tmpf61pdgjh/tmp
creating /tmp/tmpf61pdgjh/tmp/tmpf61pdgjh
creating /tmp/tmpf61pdgjh/tmp/tmpf61pdgjh/src.linux-x86_64-3.6
compile options: '-I/tmp/tmpf61pdgjh/src.linux-x86_64-3.6 -I/usr/lib/python3.6/site-packages/numpy/core/include -I/usr/include/python3.6m -c'
gcc: /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.c
In file included from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1788:0,
                 from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
                 from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.h:13,
                 from /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.c:19:
/usr/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^~~~~~~
/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.c:105:12: warning: ‘f2py_size’ defined but not used [-Wunused-function]
 static int f2py_size(PyArrayObject* var, ...)
            ^~~~~~~~~
gcc: /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.c
In file included from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/ndarraytypes.h:1788:0,
                 from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/ndarrayobject.h:18,
                 from /usr/lib/python3.6/site-packages/numpy/core/include/numpy/arrayobject.h:4,
                 from /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.h:13,
                 from /tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.c:2:
/usr/lib/python3.6/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: #warning "Using deprecated NumPy API, disable it by " "#defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
 #warning "Using deprecated NumPy API, disable it by " \
  ^~~~~~~
/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.c: In function ‘format_def’:
/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.c:139:18: warning: comparison between signed and unsigned integer expressions [-Wsign-compare]
         if (size < sizeof(notalloc)) {
                  ^
compiling Fortran sources
Fortran f77 compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /usr/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran fix compiler: /usr/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops
compile options: '-I/tmp/tmpf61pdgjh/src.linux-x86_64-3.6 -I/usr/lib/python3.6/site-packages/numpy/core/include -I/usr/include/python3.6m -c'
gfortran:fix: RK4_f2py.f90
/usr/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmpf61pdgjh/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/rkmodule.o /tmp/tmpf61pdgjh/tmp/tmpf61pdgjh/src.linux-x86_64-3.6/fortranobject.o /tmp/tmpf61pdgjh/RK4_f2py.o -L/usr/lib -lpython3.6m -lgfortran -o ./rk.cpython-36m-x86_64-linux-gnu.so
Removing build directory /tmp/tmpf61pdgjh

モジュールができたので、使ってみます。

In [14]:
import rk

intent(in) で指定した変数のみを渡します。

In [15]:
x = rk.rungekutta(0.01, 500)
plt.plot(np.arange(0,5,0.01), x)
Out[15]:
[<matplotlib.lines.Line2D at 0x7fb9a7533ef0>]

うまくできました。

誤字やおかしい点などがあったら @zawawahoge (Twitter) にお気軽にご連絡ください。

トップページに戻る