はじめに

標題の通りで、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です。