Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
33
Help us understand the problem. What is going on with this article?
@nqomorita

pythonからfortranを呼ぶ

More than 1 year has passed since last update.

はじめに

本資料は、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)してからポインタを渡すことが考えられます。
33
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ricos
FEMによる構造解析、機械学習の専門家集団。計算資源のクラウド提供もしています。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
33
Help us understand the problem. What is going on with this article?