Fortranの数値計算プログラムに、PythonのインタプリタへのAPIを埋め込んで、 計算結果の二次元配列をmatplotlibで描画し、マンデルブロ集合を書き出します。
想定するシチュエーションとして、Fortranで開発された既存の数値計算コードに対して数値計算以外の機能を追加したり、デバッグや調査のためにちょっとしたコードを差し込みたい時のことを考えます。
Fortranの数値計算コードはそこそこの規模があって、成熟しており(あるいはぐちゃぐちゃになっていてカプセル化が困難になっており...) Fortranの数値計算コード自体をいじることはしないこととします。
さて、差し込みたいコードも数値計算コードなら、そのままFortranのソースを追加すればいいのですが、たとえば指定された形式(xmlとかjsonとか...)へのファイル出力をFortranで書くのはあまり良い方法だと僕は思えません。Fortranで書きにくいからといってなんでも.csvだとか.datで保存するのはどうかと思いますよ。たいていのファイル形式にはPythonのライブラリがあるのでそれを使いましょう。
Fortranでコマンドライン引数のパーサを作るのは、地味なのに面倒です。根本的な原因として、Fortranでは「長さが異なる文字列の配列」をそのまま表現することができないせいでargvを自然に定義できず、argcを得るget_command_count argv[i]の長さおよび内容を取得するget_command_argumentを組み合わせて順にコマンドライン引数を読み込んで、そこからパースして... みたいな地味で面倒なコーディングが延々と続きます。そんなわけでFortranのプログラムはinput.txtに指定された順でオプションを書いてそれを読むことで妥協するケースが多いように思います。Pythonなら標準ライブラリargparseにまかせておけば、ヘルプメッセージまでよい感じになります。
以上のように、Fortranのプログラムで数値計算した多次元の配列をFotran以外の高水準言語に渡して処理するようなコードを書きたいというケースがあり、それに関するデモになります。今回はタイトルの通りmatplotlibを使いたいので、Python3を選びました。
デモ用なのでFortranによるマンデルブロ集合の計算自体は非常に簡単で、そもそもNumpyで書けばよいだけのお話なのですが、実際に想定するのはFortranの数値計算部分がもっと膨大だと考えてください。
このプログラムには大きく分けて以下の３つの機能があります。
- 入力(コマンドライン引数のパース)
- 数値計算(mandelbrot集合、２次元複素数配列計算)
- 出力(画面描画)
前段の説明の通り、このうち入力と出力はPythonで書いて、cythonを用いてC言語のインターフェースを作成し、Fortranのprogram mainで既存の数値計算コードと組み合わせるような形にします。
F複素平面上の領域と、実部虚部の分割数を引数に受け取り、分割数だけの整数の配列を返します。戻り値の整数は発散が確定するまでの反復回数で、反復回数の上限に達しても発散しない場合には計算を打ち切り、-1を返すこととします。 Fortran95の配列計算に関する機能を組み合わせて、並行性のある部分はOpenMPの指示行を入れて効率アップを図りました。 これくらいのコードでもPythonではできないことができてよいですね。
matplotlib.pyplot.imshowをそのまま使っています。発散しない場合(二次元配列の要素が-1のとき)黒で表示し、それ以外はmatplotlib組み込みのカラーマップでそれっぽく表示するようにしたいので、整数の二次元配列を[実部解像度,虚部解像度,RGB]の三次元配列に変換して、imshowしています。
matplotlibのfigureのデフォルトで、画像内の拡大縮小やファイル保存などの機能が揃っていて使いやすいですね。同じ機能をFortranのみで書くことを考えたら... 冷や汗が出ます。