はじめに

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

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

参考資料

PythonからFortranのサブルーチンを呼ぶ。
pythonのfortranによる高速化
 16.16 ctypes Pythonのための外部関数ライブラリ

pythonの話

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

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のデータ型の指定は、以下を参照のこと。

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のデータ型の指定は、以下を参照のこと。

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を以下のように作成します。

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

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

gfortran -shared -o libfort.so test.f90

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

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プログラムのコンパイルが実行できたことを確認して、以下のコマンドを実行します。

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で記載されていると、バグが生まれやすい状況になるかもしれません（あしからず）。

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

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