はじめに
本資料は、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）してからポインタを渡すことが考えられます。