LoginSignup
50
42

More than 5 years have passed since last update.

pythonからfortranを呼ぶ

Last updated at Posted at 2018-04-09

はじめに

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

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
50
42