pythonからfortranを呼ぶ

# はじめに

本資料は、pythonからfortranのsubroutineを呼び出す手順と、呼び出す際に気をつけておくべき事項をまとめたものです。また、pythonはpython3.6.3、fortranはfortran90（gcc7.2.0）で動作を確認しています。さて、昨今様々なプログラム言語がありますが、pythonから呼ばれるプログラムがfortranなのは、筆者がfortran過激派だからです。

pythonは（fortranに比べ）データ生成、データ入出力、可視化などが容易で、大変頼もしいです。また、筆者は普段、行列計算のプログラムを開発していますが、pythonと各種ライブラリを用いて、fortanで作成した自前のライブラリをテストするために利用しています。

# 参考資料

[PythonからFortranのサブルーチンを呼ぶ。](http://d.hatena.ne.jp/ignisan/20121017/p1)
[pythonのfortranによる高速化](http://computational-sediment-hyd.hatenablog.jp/entry/2016/04/22/200141)
[16.16 ctypes Pythonのための外部関数ライブラリ](https://docs.python.jp/3/library/ctypes.html)

# pythonの話

はじめに、python側の関数を作成します。ここでは、`main.py`を以下のように作成します。

```python
import ctypes
import numpy as np

def call_fortran(n, A):
    f = np.ctypeslib.load_library("libfort.so", ".")
    f.test_.argtypes = [
        ctypes.POINTER(ctypes.c_int32),
        np.ctypeslib.ndpointer(dtype=np.float64)
        ]
    f.test_.restype = ctypes.c_void_p

    fn = ctypes.byref(ctypes.c_int32(n))

    f.test_(fn, A)

print("** calling fortran from python")

n = 3
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.float64)

print("** write from python")
print(n)
print(A[0,0], A[0,1], A[0,2])
print(A[1,0], A[1,1], A[1,2])
print(A[2,0], A[2,1], A[2,2])

call_fortran(n, A)
```

`f.{関数名}_.argtypes`で、引数の型を指定します。
`f.{関数名}_.restype`で、関数の戻り値の型を指定します。本稿の状況の場合、void型です。
`fa = ctypes.byref(ctypes.c_int32(a))`で、定数をポインタとして渡します。

numpyのdtypeのデータ型の指定は、以下を参照のこと。

```python
a = np.array(a, dtype=np.int32) # 32 bit integer
a = np.array(a, dtype=np.int64) # 64 bit integer
a = np.array(a, dtype=np.float32) # 32 bit real
a = np.array(a, dtype=np.float64) # 64 bit real
```

ctypesのデータ型の指定は、以下を参照のこと。

```python
a = ctypes.byref(ctypes.c_int32(a)) # 32 bit integer
a = ctypes.byref(ctypes.c_int64(a)) # 64 bit integer
a = ctypes.byref(ctypes.c_float(a)) # 32 bit real
a = ctypes.byref(ctypes.c_double(a)) # 64 bit real
```

# fortranの話

次に、pythonから呼ばれるfortranの関数を作成します。ここでは、`test.f90`を以下のように作成します。

```fortran
subroutine test(n, A)
  implicit none
  integer(4), intent(in) :: n
  real(8),    intent(in) :: A(n, n)

  write(*,"(a)")"** write from fortran"
  write(*,"(i8)")n
  write(*,"(1p3e12.4)")A(1,1), A(1,2), A(1,3)
  write(*,"(1p3e12.4)")A(2,1), A(2,2), A(2,3)
  write(*,"(1p3e12.4)")A(3,1), A(3,2), A(3,3)
end subroutine test
```

作成したプログラムは、以下のコマンドでコンパイルを実行します。

```bash
gfortran -shared -o libfort.so test.f90
```

また、intel compilerでコンパイルしたプログラムを呼ぶ場合には、以下のようなコンパイラ指示行が必要です。

```fortran
subroutine test(n, A)
  !DEC$ ATTRIBUTES DLLEXPORT :: test
  implicit none
  integer(4), intent(in) :: n
  real(8),    intent(in) :: A(n, n)

  write(*,"(a)")"** write from fortran"
  write(*,"(i8)")n
  write(*,"(1p3e12.4)")A(1,1), A(1,2), A(1,3)
  write(*,"(1p3e12.4)")A(2,1), A(2,2), A(2,3)
  write(*,"(1p3e12.4)")A(3,1), A(3,2), A(3,3)
end subroutine test
```

# 実行結果

fortranプログラムのコンパイルが実行できたことを確認して、以下のコマンドを実行します。

```bash
python3 test.py 
```

実行結果は以下のとおりです。

```
** calling fortran from python
** write from python
3
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 9.0
** write from fortran
       3
  1.0000E+00  4.0000E+00  7.0000E+00
  2.0000E+00  5.0000E+00  8.0000E+00
  3.0000E+00  6.0000E+00  9.0000E+00
```

さて、この結果を踏まえ、以下の2点に注意する必要があります。

- pythonは0-origin、fortranは1-origin
    - 行列のインデックスが変わるので注意しなければなりません。

対処法としては、以下のようにfortran側で、配列の範囲を定めることが挙げられます。
ただし、既存のfortranコードが1-originで記載されていると、バグが生まれやすい状況になるかもしれません（あしからず）。

```fortran
  real(8),    intent(in) :: A(0:n-1, 0:n-1)
```

- pythonはrow major、fortranはcolumn major
    - python：メモリへ連続的にアクセスするためには、行方向に添字をインクリメントします
    - fortran：メモリへ連続的にアクセスするためには、列方向に添字をインクリメントします
    - この違いによって、プログラムいかんによっては、fortranでも大幅に演算性能が低下する場合があります
    - 対処法として、あらかじめpython側で転置（`AT = A.T`）してからポインタを渡すことが考えられます。