PythonからFortranのサブルーチンを呼ぶ。

Fortran Python

本題の前に

どうやら日記放置しすぎたようですが、それはまあ気にしない。
今回のお題はPythonです。
最近Science with PythonというコミュニティがGoogle groupにできましてですね、Twitterでのハッシュタグは#sciwpyなんですが、せっかく参加したので科学技術計算で残念ながら現役バリバリのFortranをPythonから利用する方法を書いてみよーかと。

PythonからFortranを使う利点

Pythonにもnumpyとかscipyとか便利な数値計算のためのライブラリがあるので、軽い数値計算ならPythonだけでも十分いけます。PyPyとかCythonという選択肢もありますし。
ところが、かなり重い計算とか、すでにFortranでコードがある状態だとPythonで書くor書き直す気にはならないのですよ。なぜか？

  • Pythonで速く書けても、CやFortranならもっと速くなるのではないかという疑念が憑き纏う
  • Fortranで書かれたコードがあるならそれ使えばいい。コードがどれだけ良くなっても結果が変わらないなら成果にならない。

だからヘビーな計算だと、Pythonのみってのはやる気にならない。やるならアルゴリズムのテストとかちょっとした実験とか。
とは言ってもですよ、Fortranとか色々不便なのでPythonを使いたいわけですよ。手軽だし。というわけで、

  • 既存のFortranコードをライブラリ化してPythonから利用する
  • 処理の重い部分だけをFortranで書く

という方針がたちます、これによって、

  • 今までのコードを大きく書きなおさずに利用可能
  • Fortranコードのプリプロセス/ポストプロセス(例えば、計算のバッチ処理化とか計算結果のグラフ出力とかいろいろ)をPythonで手軽に記述でウマー

というメリットがあるわけですね！

どうやって呼び出すのか？

過去にはf2pyとかいうものもあったみたいですが、なんだか更新されてないようです。これについては今回はスルーして、Pythonの標準ライブラリctypesを使います。あと、以下コンパイラ依存が強くなるので話をgcc/gfortranに限定します。今回はPython2.7でやります。
ctypesはPython2.5から追加された、C言語の関数を呼び出すためのライブラリです。C言語とFortranを組み合わせたプログラムを書いたことのある人は知ってると思いますが、実はC言語からFortranの関数(サブルーチン)を呼び出すのは結構簡単です。C言語からFortranのサブルーチンを呼ぶときの注意点は

  • サブルーチン名の末尾にアンスコ"_"を付けたものがC言語から呼び出す時の関数名となる。
  • サブルーチンの返値はvoid型
  • 引数は参照(ポインタ)渡し

といったところです。C言語とFortranは簡単に連携でき、C言語とPythonは標準ライブラリctypesで連携できる！
ならば、C言語でFortranのラッパーを書いてやればいける！



と思うかもしれませんが(最初私はそう思いましたが)、C言語でラッパーを書く必要はありません。ctypesでFortran直接呼べます。少なくともgcc/gfortranであれば。

Example 1 基本編

具体的に、次のようなFortranの関数をPythonから呼び出すことを考えます。

! add.f90
subroutine add(a,b)
    implicit none
    integer(8),intent(in) :: a
    integer(8),intent(inout):: b
    c = a + b
end subroutine

やってることは簡単です。bにaを加えるだけです。intentとか無くても動きますがつけます。注意事項ですが、integerのサイズ(32bitか64bitか)が環境依存なので、きちんと指定したほうが良いかと思います。この記事では整数・浮動小数いずれも64bitです。このソースを共有ライブラリとしてビルドします。

$ gfortran -shared -fPIC -o add.so add.f90

これはLinuxの場合の例です。Windowsの場合はfPICを消して、soではなくdllにして下さい。

次にPythonでのコードです。

#!/usr/bin/env python 
#vim:fileencoding=utf8  
from ctypes import * 

addmodule = cdll.LoadLibrary("add.so")

# 引数の型を指定する。今回はintのポインタ型
# Fortranのサブルーチン名"add"にアンスコ"_"をつける(gccでは)
addmodule.add_.argtypes = [ POINTER(c_int64),POINTER(c_int64)]
      
# 戻り値の型を指定する。Fortranのサブルーチンはvoidしか返せない
addmodule.add_.restype = c_void_p  

# 呼び出しに使う引数はctypesの型でラップする
a,b = 10,8 
a = c_int64(a)
b = c_int64(b)
 
# byrefでポインタにして渡す         
addmodule.add_(byref(a),byref(b)) 
print b.value # 18

詳しい説明はソースコード中に書いたので省きます。実はC言語(というかポインタ)が不慣れなウチの後輩がハマっていたので、それを紹介するためにa,bの初期化で若干回りくどい書き方をしています。というのは、この例ではc_intでラップした変数をbyrefを利用して参照で渡していますが、もし

a,b = 10,8  
addmodule.add_(byref(c_int64(a)),byref(c_int64(b)))
print b

と呼び出した場合はどうなるか？ これ結果は8になります。わかっている人には何でもないことですが、c_int64でラップすると別のオブジェクトが作られるので、その参照をとって渡してしまうと意図した結果が得られません。注意。
あと、よくわかりませんが実はbyref使わなくても動いてしまいました。謎。
byrefではなくpointerを使う方法もありますが、速さで言えばbyrefの方が速いそうです。詳しくはドキュメントへGoってことで。

とまあ、こんな感じで直接Fortanのサブルーチンを呼び出すことができます。
まとめると、

  1. Fortranサブルーチンを共有ライブラリとしてコンパイル
  2. Pythonではライブラリのロード
  3. 関数プロトタイプ(引数、返り値の型)を指定
  4. ctypesのデータ型でラップして呼び出し。Fortranは参照渡しな点に注意。

です！

Example 2 実践編

Fortranのサブルーチンを呼べるようになったわけですが、実際のところ、int型の変数1つとか渡さないわけですよ、数値計算やってる人は。渡すのは主に配列。しかも、Pythonとの連携をするならPython側はnumpyの配列使ってることでしょう。だから、numpyの配列を渡すとこまでやらないといけませんよね？

というわけで、基本編のnumpy配列版。まずはFortranソースコードから

! add_np.f90
#define MAXSIZE 1024 
subroutine add_array(a,b,N)
    implicit none
    integer(8),intent(in) :: N  
    real(8),dimension(0:MAXSIZE),intent(in) :: a
    real(8),dimension(0:MAXSIZE),intent(inout):: b
#ifndef NDEBUG
    if (N>MAXSIZE) then
        print *,"MAXSIZE is too small."
    endif 
#endif 
    b(0:N-1) = b(0:N-1) + a(0:N-1)
end subroutine

Fortranではほとんど使われない(というかFortranじゃない)プリプロセッサを使って配列のサイズチェックを行なってます。gfortranではコンパイルオプションに-cppを追加することで、C言語のプリプロセッサが使えるようになります。
でもBOOST.PPは使えなかったはず...

なんでこんなことするのかというと、Pythonから渡された配列は型としてサイズがわからないからですね。こんな風に大きめに配列をとって、もし小さければ警告を出すなどの工夫が必要になります。プリプロセッサならデバッグが終わったら外せるので便利です。ただし、コンパイラ依存？
渡された配列サイズがわからないので、引数でサイズをとるようにしています。C言語と同じです。
この辺り、もっとマシな方法があればいいのですが・・誰か教えて下さい。

次にPythonの方ですね。numpyにはctypeslibがあるのでそれを使うことにします。使い方はほとんどctypesと同じです。

#!/usr/bin/env python
#vim:fileencoding=utf8
from ctypes import *
import numpy as np

add_np = np.ctypeslib.load_library("add_np.so",".")
add_np.add_array_.argtypes = [
               np.ctypeslib.ndpointer(dtype=np.float64),
               np.ctypeslib.ndpointer(dtype=np.float64),
               POINTER(c_int64),]
add_np.add_array_.restype = c_void_p   
                    
a = np.arange(0.,10.,dtype=np.float64) # 0,1,2,3,4,5,..,9 
b = a*2                                # 0,2,4,6,8,...,18
size = byref(c_int64(b.size)) 
add_np.add_array_(a,b,size) # ndarrayはそのまま渡して良い 
print b # [ 0 3 6 9 12 .. 27 ]

ほぼExample1と同じですね。numpyの配列は大抵ndarrayとして扱うことになりますので、Fortranにはそのポインタを渡してやります。それがnp.ctypeslib.ndpointer(dtype=np.float64)の部分。私はかなり変数の型に注意を払ってます。
コメントにも書きましたが、ndarrayな配列はそのままFortranのサブルーチンに渡してもOK。
整数型はNG。面倒な場合は整数もndarrayで渡してしまうという手も・・。


numpy配列の場合もあまり説明することないですね。(久しぶりに書いたから疲れたわけではないですよ決して)
こんな感じで、Python+numpy/scipy+Fortran/C,C++が簡単にできるのでScienceな分野だと結構Python強いんじゃないかなーという気がします。

みんなでScience with Pythonしましょう！