はじめに
標題の通りで、ctypeslibを使ってFortranのサブルーチン(あるいはCの関数)を呼び出すとき、numpy.ndarrayのviewを渡すと正しく動作せず、場合によってはSegmentation faultでプログラムが落ちることがあります。
numpy.ndarrayを操作したときにcopy(コピー)が作成されるのかview(ビュー)が作成されるのかの区別と、メモリに確保された配列をFortran(あるいはC)がどのようにアクセスしているかを考えると当たり前ですが、うっかりハマって時間を溶かしたのでまとめました。
サンプルコード
Fortran側の準備
まず、以下のように3次元配列を受け取ってその要素を標準出力するサブルーチンを作成します。
module mod
contains
subroutine print_array3d(array3d, imax, jmax, kmax) bind(C)
use, intrinsic :: ISO_C_BINDING
integer(c_int32_t), intent(in) :: imax, jmax, kmax
real(c_float), intent(in) :: array3d(imax, jmax, kmax)
integer(c_int32_t) :: i, j, k
do k = 1, kmax
do j = 1, jmax
do i = 1, imax
write(*,*) i, j, k, array3d(i, j, k)
end do
end do
end do
end subroutine print_array3d
end module mod
このプログラムを共有ライブラリとしてコンパイルします。
gfortran -cpp -fPIC -shared mod.f90 -o mod.so
Python
実際にPythonで試してみます。Jupyter Notebook上で実行すると便利です。
まず、必要なパッケージをインポートしておきます。
from ctypes import POINTER, c_int32, c_float, c_void_p
import numpy as np
次にテスト用に3次元配列を用意します。
imax = 2
jmax = 2
kmax = 2
array3d = np.empty((kmax, jmax, imax), dtype=np.float32)
count = 0
for k in range(kmax):
for j in range(jmax):
for i in range(imax):
array3d[k, j, i] = count
count += 1
用意したarray3dの中身は当然以下のようになっています。
array([[[0., 1.],
[2., 3.]],
[[4., 5.],
[6., 7.]]], dtype=float32)
共有ライブラリの関数を呼び出す準備として、ライブラリのロードと引数と返値の型の設定を行います。
mod = np.ctypeslib.load_library("mod.so", ".")
mod.print_array3d.argtypes = [
np.ctypeslib.ndpointer(dtype=np.float32),
POINTER(c_int32),
POINTER(c_int32),
POINTER(c_int32),
]
mod.print_array3d.restype = c_void_p
実際にmod.print_array3dにarray3dを渡してみます。配列サイズを表す変数は、c_int32のポインタとして渡す点が注意が必要です。
_kmax = POINTER((c_int32))(c_int32(kmax))
_jmax = POINTER((c_int32))(c_int32(jmax))
_imax = POINTER((c_int32))(c_int32(imax))
res = mod.print_array3d(array3d, _imax, _jmax, _kmax)
1 1 1 0.00000000
2 1 1 1.00000000
1 2 1 2.00000000
2 2 1 3.00000000
1 1 2 4.00000000
2 1 2 5.00000000
1 2 2 6.00000000
2 2 2 7.00000000
期待通り、0から7までが表示されました。
次に、Python側で配列のスライスを使って順番を入れ換えたいと思います。
array3d[::-1,:,:]
array([[[4., 5.],
[6., 7.]],
[[0., 1.],
[2., 3.]]], dtype=float32)
これをmod.print_array3dに渡したときに、4,5,6,7,0,1,2,3の順に表示されることを期待したくなります。実際にはどうなるでしょうか。
res = mod.print_array3d(array3d[::-1,:,:], _imax, _jmax, _kmax)
1 1 1 4.00000000
2 1 1 5.00000000
1 2 1 6.00000000
2 2 1 7.00000000
1 1 2 -1.57098575E-34
2 1 2 5.96014276E-41
1 2 2 6.00097418
2 2 2 7.00000000
k = 2の場合(Python側ではk = 1)に、値がおかしくなりました。これは配列外参照を起こしていて、場合によってはSegmentation faultでプログラムが落ちます。配列のスライスでは配列のcopyではなくviewが作成されるため、Fortranのサブルーチンにviewのポインタを渡してしまうと、正しく読めないようです。
これを回避するためには、配列のコピーを生成してから渡せばよいです。
res = mod.print_array3d(array3d[::-1,:,:].copy(), _imax, _jmax, _kmax)
1 1 1 4.00000000
2 1 1 5.00000000
1 2 1 6.00000000
2 2 1 7.00000000
1 1 2 0.00000000
2 1 2 1.00000000
1 2 2 2.00000000
2 2 2 3.00000000
なお、配列がviewであるかどうかはbase属性を調べることで分かります。array3d.base is NoneはTrueなので、array3dはviewではありません。一方array3d[::-1,:,:].base is NoneはFalseのため、こちらはviewです。