FortranプログラムをPythonから活用する

第6回：ctypesとbind(C)によるFortranコードのラッピング

前回はf2pyを使い、FortranコードをPythonから呼び出す手法を探りました。f2pyはラッパーを自動生成できる便利なツールですが、環境構築の複雑さや、Python/NumPyのバージョンアップに伴う再ビルドが頻繁に必要になるという実運用上の考慮すべき点があることが分かりました。

そこで今回は、外部ツールへの依存をなくし、より堅牢な連携を実現するため、Fortran 2003標準のbind(C)とPython標準モジュールctypesを組み合わせるアプローチを検証します。この手法はFortranコンパイラとPython本体の標準機能だけで実現できます。

bind(C)とctypesによる連携の仕組み

このアプローチでは、両言語が標準で備えるC言語連携機能を架け橋として利用します。

  • bind(C) (Fortran側): Fortran 2003の標準機能。Fortranのサブルーチンや関数を、C言語から呼び出せる形式で共有ライブラリ（.dllや.so）へ公開します。
  • ctypes (Python側): Pythonの標準モジュール。共有ライブラリを読み込み、内部のC言語互換関数をPythonから直接呼び出します。

この2つにより、f2pyのような追加ツールなしでFortranとPythonを連携させます。

実装プロセス

連携は以下の3ステップで実現します。

  1. Fortranラッパーの作成: bind(C)を使い、C言語互換のFortran関数を作成する。
  2. 共有ライブラリの生成: 元のコードとラッパーをコンパイルし、共有ライブラリを生成する。
  3. Pythonラッパーの作成: ctypesを使い、共有ライブラリの関数を呼び出せるようにする。

ステップ1：bind(C)によるFortranラッパーの作成

元のdsvdcサブルーチンを変更せず、C言語から直接呼び出せるようにbind(C)でラッパーを作成します。このラッパーで、FortranとC言語間にある以下の非互換性を吸収します。

  • 呼び出し規約: 引数の渡し方やシンボル名の違い。
  • データ型: INTEGERとintなどがメモリ上で同一とは限らない点。
  • API設計の差異を吸収: Fortran（引数で結果を返す）とC（戻り値でステータスを返す）の慣習の違い。

具体的に以下の実装を行いました。

  • C互換インターフェース: function ... bind(c, name=...)で、関数をC言語から直接呼出し可能な形で公開。
  • C互換のデータ型: iso_c_bindingモジュールのc_intやc_doubleを使い、C言語の型と正確に対応。
  • 戻り値の設計: 元のサブルーチン（戻り値なし）とは異なり、C言語と連携しやすいようにinfoコードを返す関数(FUNCTION)として定義。
  • 元のサブルーチンの呼び出し: ラッパー内部で、オリジナルのdsvdcの呼び出し。

以下に、今回作成したFortranラッパーコードを示します。

Fortranラッパーコード (dsvdc_wrapper.f90)

module dsvdc_wrapper_mod
    use svd, only: dsvdc, dp
    use, intrinsic :: iso_c_binding, only: c_int, c_double
    implicit none
    private
    public :: dsvdc_c

contains
    ! C言語互換のラッパー関数（bind(c, ...)を指定する）
    function dsvdc_c(x, n, p, s, e, u, v, job) result(info) bind(c, name='dsvdc_c')

        ! --- 引数定義 (C言語の型に対応させる) ---
        integer(c_int), value, intent(in) :: n, p, job
        real(c_double), intent(inout) :: x(n, p)
        real(c_double), intent(out) :: s(min(n+1,p)), e(p), u(n, n), v(p, p)
        integer(c_int) :: info

        ! --- 元のFortranサブルーチン呼び出し ---
        integer :: info_fortran
        call dsvdc(x, n, p, s, e, u, v, job, info_fortran)

        ! --- 結果を関数の戻り値として設定 ---
        info = info_fortran
    end function dsvdc_c
end module dsvdc_wrapper_mod

ステップ2：共有ライブラリの生成

次に、dsvdc.f90とdsvdc_wrapper.f90をコンパイルし、単一の共有ライブラリ（Windowsでは.dll、Linux/Macでは.so）を作成します。

今回指定したgfortranでのコマンドは以下の通りです。

# 1. 元のdsvdc.f90をオブジェクトファイルにコンパイル
gfortran -O3 -Wall -c fortran_src/ctypes/dsvdc.f90 -o build/dsvdc_ctypes.o

# 2. ラッパーとオブジェクトファイルをリンクし、共有ライブラリ(.dll)を生成
gfortran -O3 -Wall -shared fortran_src/ctypes/dsvdc_wrapper.f90 build/dsvdc_ctypes.o -o bin/svd_ctypes.dll

これにより、Pythonからロード可能な共有ライブラリ svd_ctypes.dll が生成されます。

ステップ3：ctypesによるPythonラッパーの作成

最後に、共有ライブラリを利用するためのPythonラッパー(svd_ctypes.py)を作成します。このラッパーの主目的は、共有ライブラリのロード、あるいは関数の引数や戻り値の型を指定するといった、ctypesを直接扱う上で必要となる煩雑な手続きを隠蔽することです。これにより、利用者はこれらの低レベルな準備作業を意識する必要がなくなり、手軽にFortranの機能を利用できるようになります。

具体的には、このラッパーで以下の設定を行います：

  • ライブラリのロード: ctypes.CDLLで共有ライブラリをロードする。
  • 戻り値の型の指定: restypeで、Fortran関数の戻り値がC言語の整数型であることを指定する。
  • 引数の型の指定: argtypesに、Fortranで定義した順序でC言語互換の型リストを指定する。
    • NumPy配列はnp.ctypeslib.ndpointerでC言語のポインタとして渡す。
    • Fortranは列優先(Column-major)配列を扱うため、flags='F_CONTIGUOUS'の指定が必須となる。
  • Fortran互換配列の準備: Fortranに渡す配列はnp.asfortranarray()で列優先に揃える。
  • 出力用配列の確保: Fortran側で書き込む配列は、呼び出し側（Python）でnp.zerosなどを使って事前にメモリ確保する。
  • dsvdc関数の実装: C互換関数を呼び出す際に、以下の処理を行う。
    • 入力配列xが列優先 (Fortran Order) のメモリレイアウトを持つよう、必要に応じてnp.asfortranarray()で変換する。
    • 関数呼び出し後、出力用配列uとvがPython側の元の配列と異なる場合（メモリレイアウトの変換が行われた場合）、データをコピーして戻す。

以下に、今回作成したPythonラッパーコードを示します。

Pythonラッパーコード (svd_ctypes.py) 抜粋

import numpy as np
import ctypes
from pathlib import Path

# 共有ライブラリをロード
dll_path = Path(__file__).parent / 'svd_ctypes.dll'
fortran_lib = ctypes.CDLL(str(dll_path))

# Fortran関数への参照を取得
dsvdc_c_func = fortran_lib.dsvdc_c

# 関数の戻り値の型を指定 (C言語のint)
dsvdc_c_func.restype = ctypes.c_int

# 関数の引数の型をリストで指定
dsvdc_c_func.argtypes = [
    np.ctypeslib.ndpointer(dtype=np.float64, flags='F_CONTIGUOUS'), # x
    ctypes.c_int,                                                   # n
    ctypes.c_int,                                                   # p
    np.ctypeslib.ndpointer(dtype=np.float64),                       # s
    np.ctypeslib.ndpointer(dtype=np.float64),                       # e
    np.ctypeslib.ndpointer(dtype=np.float64, flags='F_CONTIGUOUS'), # u
    np.ctypeslib.ndpointer(dtype=np.float64, flags='F_CONTIGUOUS'), # v
    ctypes.c_int,                                                   # job
]

def dsvdc(x, n, p, s, e, u, v, job):
    """Fortranのdsvdc_c関数を呼び出すPythonインターフェース"""
    # 入力配列をFortran互換の列優先に変換
    x_f = np.asfortranarray(x)
    u_f = np.asfortranarray(u)
    v_f = np.asfortranarray(v)
    
    # C互換のFortran関数を呼び出す
    info = dsvdc_c_func(x_f, n, p, s, e, u_f, v_f, job)
    
    # 必要に応じて、出力用配列u, vの内容を元の配列にコピーし直す
    # (asfortranarrayがコピーを生成した場合に備える)
    if u is not u_f:
        u[...] = u_f
    if v is not v_f:
        v[...] = v_f
        
    return info

実行と検証

前述の手順で実現したFortranとPythonの連携が正しく機能するかを検証します。ここでは、同じ入力データ（1600x1600の行列）を用いて特異値分解を実行し、結果の正確性とパフォーマンスをオリジナルのFortranプログラムと比較します。

今回の検証では以下の検証用Pythonコード（抜粋）を使用しました。このコードは、ファイルから行列を読み込み、Fortran側で必要となる出力用配列を事前に確保した上で、作成したラッパー関数を実行します。

呼び出しの主要部分は以下の通りです。

# ファイルから行列'x'を列優先(order='F')で読み込む
x, n = get_matrix_from_file(data_file)
# ...

# --- Fortranサブルーチン呼び出しのための配列を準備 ---
p = n
job = 11

# 結果を格納する配列をあらかじめ確保する
s = np.zeros(p, dtype=np.float64)
e = np.zeros(p, dtype=np.float64)
u = np.zeros((n, n), dtype=np.float64, order='F')
v = np.zeros((p, p), dtype=np.float64, order='F')

# --- 時間を計測してSVD計算を実行 ---
start_time = time.perf_counter()

# 作成したctypesラッパー関数を呼び出す
info = svd_ctypes.dsvdc(x, n, p, s, e, u, v, job)

end_time = time.perf_counter()

# --- 結果を出力 ---

実行結果と考察

上記コードと、比較対象であるオリジナルのFortranプログラムを実行した結果は以下の通りです。

1. bind(c) + ctypes経由でのPythonからの実行結果

> python benchmark_svd_python_ctypes.py matrix_1600x1600_rand_uniform.bin
13.21040380  # 実行時間（秒）
Top 5 singular values (Python ctypes):
   23.007873050107914    22.898887126396772    22.847291879705001    22.725199286295815    22.706511528915520

First 3 components of leading eigenvector (1st column of U):
    0.015510246530786     0.042426665711607     0.010966853963481

2. オリジナルFortranプログラムの実行結果

> benchmark_svd_fortran_original.exe matrix_1600x1600_rand_uniform.bin
 13.39000034  # 実行時間（秒）
Top 5 singular values (Fortran Original):
   23.007873050107914    22.898887126396772    22.847291879705001    22.725199286295815    22.706511528915520

First 3 components of leading eigenvector (1st column of U):
    0.015510246530786     0.042426665711607     0.010966853963481

実行結果が示す通り、bind(C)とctypesを用いたアプローチで得られた特異値および固有ベクトルは、オリジナルのFortranプログラムの結果と数値的に完全に一致しました。また、実行時間もほぼ同等であり、Pythonからの呼び出しによる目立ったオーバーヘッドは見られません。

まとめと次回予告

今回は、Fortran 2003標準のbind(C)とPython標準ライブラリctypesを使い、外部ツールに頼らずにFortranとPythonを連携させる方法を見てきました。この手法は、ラッパーコードを手書きする手間と引き換えに、Pythonのバージョンアップなど環境の変化にも強い、堅固な連携が実現できるのが利点です。

さて、これまで6回にわたり、Fortran資産を現代のPythonエコシステムで活用するための様々なアプローチ、「完全移行」と「ハイブリッド連携」の2つの方針を、具体的なコードと共に探ってきました。

次回、第7回「総括 — あなたの状況に最適なFortran資産活用法は？」では、これまでの各アプローチを「性能」「開発工数」「保守性」などの観点から整理・比較し、ご自身の状況に合わせて最適な手法を検討する際の、判断材料を提供できればと考えています。


関連情報
製品
その他