ctypeslibを使ってFortranを呼び出すときにnumpy.ndarrayのviewを渡してはいけない

# はじめに

標題の通りで、ctypeslibを使ってFortranのサブルーチン(あるいはCの関数)を呼び出すとき、numpy.ndarrayのviewを渡すと正しく動作せず、場合によってはSegmentation faultでプログラムが落ちることがあります。

numpy.ndarrayを操作したときにcopy(コピー)が作成されるのかview(ビュー)が作成されるのかの区別と、メモリに確保された配列をFortran(あるいはC)がどのようにアクセスしているかを考えると当たり前ですが、うっかりハマって時間を溶かしたのでまとめました。

# サンプルコード

## Fortran側の準備

まず、以下のように3次元配列を受け取ってその要素を標準出力するサブルーチンを作成します。

```Fortran
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
```

このプログラムを共有ライブラリとしてコンパイルします。

```shell
gfortran -cpp -fPIC -shared mod.f90 -o mod.so
```

## Python

実際にPythonで試してみます。Jupyter Notebook上で実行すると便利です。

まず、必要なパッケージをインポートしておきます。

```python
from ctypes import POINTER, c_int32, c_float, c_void_p
import numpy as np
```

次にテスト用に3次元配列を用意します。

```python
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)

共有ライブラリの関数を呼び出す準備として、ライブラリのロードと引数と返値の型の設定を行います。

```python
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`のポインタとして渡す点が注意が必要です。

```python
_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側で配列のスライスを使って順番を入れ換えたいと思います。

```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の順に表示されることを期待したくなります。実際にはどうなるでしょうか。


```python
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のポインタを渡してしまうと、正しく読めないようです。

これを回避するためには、配列のコピーを生成してから渡せばよいです。

```python
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です。