1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@rino_yume(rino yume)

Fortran 2018をPythonから使う ― NumPy 2.5.4のf2pyとMesonで数値計算を組み込む

1
Posted at

Fortran 2018をPythonから使う ― NumPy 2.5.4のf2pyとMesonで数値計算を組み込む

はじめに

Pythonはデータの読み込みや加工、グラフの作成、Web APIの開発に便利です。一方、河川の流れや氾濫、降雨流出のような数値計算では、Fortranで書かれたプログラムを使う機会が今でも多くあります。

既存のFortranプログラムをPythonで書き直す必要はありません。計算部分はFortranのまま使い、入力データの準備や計算結果の整理をPythonに任せる方法があります。その橋渡しをするのが、NumPyに付属する**f2py（Fortran to Python interface generator）**です。

今回は、2026年10月15日時点で公開されているNumPy 2.5.4を対象に、WSL2のUbuntu 24.04で使う手順をまとめます。簡単な配列計算から始め、2次元メッシュの計算、インターフェースの定義、Fortran 2018との付き合い方、最適化、Docker環境まで取り上げます。

f2pyで気を付けたいのは、最近のバージョンでビルド方法が変わっていることです。以前の記事にあるnumpy.distutils前提の手順は、現在のNumPyではそのまま使えない場合があります。

1. f2pyでできること

f2pyはFortranのサブルーチンや関数を解析し、Pythonから呼び出すための拡張モジュールを生成します。NumPy配列をFortranの配列引数に渡せるため、数値計算で使いやすい仕組みです。

Python+f2py+Fortran_workflow.png

例えば、時間積分やセルごとの更新をFortranで計算し、Python側ではNumPy配列の作成やGeoTIFFの読み書き、FastAPIとの連携を担当する構成にできます。

ただし、f2pyを使うだけで処理が高速化するわけではありません。単純な配列加算ならNumPyだけで十分に高速です。複雑なループや繰り返し計算をFortranにまとめて渡すと、PythonとFortranを往復する回数を減らせます。

NumPy 2.5で変わった点

NumPy 1.26系でMesonを使うビルド方法が導入され、NumPy 2.5では従来のdistutilsバックエンドが削除されました。現在はMesonを使ってビルドします。

項目 現在の扱い
f2py本体 NumPyに同梱されています
ビルド Mesonを使用します
Fortranコンパイラ gfortranなどを使用します
Python側のデータ 主にNumPy配列を使用します
インターフェース ソースから自動生成するか、.pyfで定義します
Fortran 2018 計算内部で使用できますが、すべての構文を自動公開できるわけではありません

基本的なビルドコマンドは次の形です。

python -m numpy.f2py -c -m fastcalc add_vectors.f90

-cがコンパイル、-m fastcalcがPythonに読み込むモジュール名です。python -m numpy.f2pyと書いておくと、使用中のPython環境にインストールされたNumPyを明示的に使えます。

2. WSL2 Ubuntu 24.04に環境を作る

想定する環境はWindows 11、WSL2 Ubuntu 24.04、Python 3.12、gfortranです。Windows側のVS CodeからWSL2のフォルダーを開いて編集しても構いません。

まず、システム側のツールをインストールします。

sudo apt update
sudo apt install -y \
    python3 python3-venv python3-dev \
    build-essential gfortran ninja-build pkg-config

次に、Pythonの仮想環境を作成します。

mkdir -p ~/f2py_project
cd ~/f2py_project

python3 -m venv .venv
source .venv/bin/activate

python -m pip install --upgrade pip
python -m pip install numpy==2.5.4 meson pytest

確認用のコマンドです。

python --version
python -c "import numpy; print(numpy.__version__)"
python -m numpy.f2py -v
gfortran --version
meson --version
ninja --version

ここではNumPyのバージョンを固定しています。ほかのNumPy 2.xを利用する場合も、実際に使うバージョンの公式マニュアルを確認してください。

3. 最初の例：2つの配列を加算する

まずは、2つの配列を足し合わせるサブルーチンをPythonから呼び出します。

作業ディレクトリは次の構成です。

f2py_project/
├── .venv/
├── add_vectors.f90
└── test_add.py

3.1 Fortranのコード

add_vectors.f90を作成します。

subroutine add_vectors(n, a, b, c)
    implicit none

    integer, intent(in) :: n
    double precision, intent(in) :: a(n), b(n)
    double precision, intent(out) :: c(n)
    integer :: i

    do concurrent (i = 1:n)
        c(i) = a(i) + b(i)
    end do
end subroutine add_vectors

do concurrentはFortran 2008で導入された構文で、Fortran 2018でも利用できます。ここでは各要素の計算が独立していることを示しています。ただし、これを書いただけで必ずSIMD化やマルチスレッド化が行われるわけではありません。

intent(in)は入力、intent(out)は出力です。f2pyはこれらの宣言を読み取り、Python向けの引数や戻り値を用意します。

3.2 拡張モジュールを作る

同じフォルダーで次のコマンドを実行します。

python -m numpy.f2py -c -m fastcalc add_vectors.f90

成功すると、Linuxでは次のような名前のファイルができます。PythonのバージョンやCPU環境によって実際のファイル名は変わります。

fastcalc.cpython-312-x86_64-linux-gnu.so

3.3 Pythonから呼び出す

test_add.pyを作成します。

"""Fortranの配列加算をPythonから呼び出す例。"""

import numpy as np
import fastcalc


def main() -> None:
    """配列を渡して、計算結果をNumPyで検証する。"""
    a = np.array([1.0, 2.0, 3.0], dtype=np.float64)
    b = np.array([10.0, 20.0, 30.0], dtype=np.float64)

    c = fastcalc.add_vectors(a, b)

    print("a =", a)
    print("b =", b)
    print("c =", c)

    np.testing.assert_allclose(c, a + b)


if __name__ == "__main__":
    main()

実行します。

python test_add.py

期待する結果は次のとおりです。

a = [1. 2. 3.]
b = [10. 20. 30.]
c = [11. 22. 33.]

Fortran側の引数には配列長nがありますが、Pythonでは指定していません。f2pyが配列の大きさを取得して渡せるためです。

生成された引数の説明は、Pythonから確認できます。

python -c "import fastcalc; print(fastcalc.add_vectors.__doc__)"

f2pyがどのような引数を公開したか迷ったときは、まずこの説明を見ると分かりやすいです。

4. intentとNumPy配列の受け渡し

Fortranの引数宣言は、Pythonからの使い方にも影響します。

Fortranの宣言 Python側での扱い
intent(in) 入力として渡します
intent(out) 多くの場合、戻り値として受け取ります
intent(inout) 入出力用の配列として渡します
dimension(n) 1次元配列として扱います
dimension(nx,ny) 2次元配列として扱います
double precision 通常はnp.float64に対応します

はじめのうちは、入力と出力を別々にしておく方が実装しやすいと思います。intent(inout)の場合、型やメモリ配置によってコピーや書き戻しが関係するためです。

特に2次元以上の配列では、NumPyとFortranのメモリ配置の違いを意識する必要があります。

5. 2次元メッシュをFortranで計算する

もう少し実際の数値計算に近い例として、2次元の拡散方程式を1ステップだけ計算します。

\frac{\partial h}{\partial t}
=D\left(\frac{\partial^2 h}{\partial x^2}
+\frac{\partial^2 h}{\partial y^2}\right)

ここでは格子間隔を等しくした5点差分を使います。これはf2pyの配列受け渡しを確かめるための例であり、洪水氾濫の浅水流方程式ではありません。

5.1 Fortran側の計算処理

diffusion.f90を作成します。

diffusion.f90（クリックして開く） 
subroutine diffusion_step(nx, ny, dt, dx, diffusivity, h, hnew)
    implicit none

    ! メッシュ数と計算条件
    integer, intent(in) :: nx, ny
    double precision, intent(in) :: dt, dx, diffusivity

    ! 現在時刻と次時刻の配列
    double precision, intent(in) :: h(nx, ny)
    double precision, intent(out) :: hnew(nx, ny)

    integer :: i, j

    ! 境界セルは現在の値を維持する
    hnew = h

    ! 内部セルを5点差分で更新する
    do concurrent (j = 2:ny-1, i = 2:nx-1)
        hnew(i, j) = h(i, j) &
            + dt * diffusivity / (dx * dx) * ( &
                h(i-1, j) + h(i+1, j) &
              + h(i, j-1) + h(i, j+1) &
              - 4.0d0 * h(i, j) )
    end do
end subroutine diffusion_step

拡散係数をD、時間刻みをdt、格子間隔をdxとすると、この陽解法の安定条件は、等間隔格子の場合、おおむね次のようになります。

\frac{D\,\Delta t}{\Delta x^2}\leq\frac{1}{4}

以下の例ではD=0.1、dt=1.0、dx=1.0なので、この条件を満たします。

5.2 2つのFortranソースをまとめてビルド

python -m numpy.f2py -c -m fastcalc \
    add_vectors.f90 diffusion.f90

5.3 Pythonから計算する

test_diffusion.pyを作成します。

test_diffusion.py（クリックして開く） 
"""2次元拡散方程式の1ステップ計算を検証する。"""

import numpy as np
import fastcalc


def main() -> None:
    """初期値を設定し、Fortranの計算結果を確認する。"""
    nx, ny = 100, 100
    dt = 1.0
    dx = 1.0
    diffusivity = 0.1

    # Fortranと同じ列優先のメモリ配置にする
    h = np.zeros((nx, ny), dtype=np.float64, order="F")
    h[nx // 2, ny // 2] = 1.0

    hnew = fastcalc.diffusion_step(
        dt, dx, diffusivity, h
    )

    print("配列の形状:", hnew.shape)
    print("中央セル:", hnew[nx // 2, ny // 2])
    print("総和:", hnew.sum())

    assert hnew.shape == (nx, ny)
    np.testing.assert_allclose(
        hnew[nx // 2, ny // 2], 0.6, atol=1.0e-12
    )
    np.testing.assert_allclose(
        hnew.sum(), 1.0, atol=1.0e-12
    )


if __name__ == "__main__":
    main()

python test_diffusion.py

今回の初期条件では、中央セルは0.6、配列全体の総和は1.0になります。中央セルの周囲4セルにそれぞれ0.1ずつ移るためです。

なお、境界セルを固定するこの単純な処理が、どのような条件でも総量を保存するわけではありません。保存性を要求する実際の水理計算では、フラックスや境界条件を別途検証する必要があります。

6. C-orderとF-orderの違い

NumPyで通常作る2次元配列はC-order（行優先）です。一方、Fortranの配列はF-order（列優先）になっています。

# 一般的なNumPy配列
h_c = np.zeros((100, 100), dtype=np.float64)

# Fortranへの受け渡しを意識した配列
h_f = np.zeros((100, 100), dtype=np.float64, order="F")

# 既存配列を必要に応じて変換する
h_f2 = np.asfortranarray(h_c)

print(h_c.flags.f_contiguous)
print(h_f.flags.f_contiguous)

Fortranでは先頭の添字が連続するため、通常は次のループ順序がメモリアクセスに適しています。

do j = 1, ny
    do i = 1, nx
        h(i, j) = h(i, j) + 1.0d0
    end do
end do

f2pyは必要に応じてメモリ配置を変換してくれますが、変換のためにコピーが発生することがあります。例えば、1000×1000の倍精度配列は約8 MBです。これを時間ステップごとに何枚もコピーすると、計算以外の処理時間が増えてしまいます。

コピーが発生するかを調べたいときは、f2pyの配列コピーレポート機能も使えます。詳細は公式のUsing F2PYにあります。

また、GeoTIFFなどのラスターデータでは、配列の(行, 列)とモデルの(x, y)の対応にも注意します。F-orderに変更することと、座標の向きを入れ替えることは別です。

7. .pyfファイルで公開する関数を選ぶ

実務のFortranプログラムには、内部計算用のサブルーチンがたくさんあります。それらをすべてPythonに公開する必要はありません。

f2pyのシグネチャファイル.pyfを使うと、公開する関数や引数を明示できます。

まず、インターフェースのひな型を生成します。

python -m numpy.f2py -h fastcalc.pyf -m fastcalc \
    add_vectors.f90 diffusion.f90

例えば、配列加算のみを公開する場合のfastcalc.pyfは次のように書けます。

! -*- f90 -*-
python module fastcalc
    interface
        subroutine add_vectors(n, a, b, c)
            integer intent(hide), depend(a) :: n = shape(a, 0)
            double precision dimension(n), intent(in) :: a
            double precision dimension(n), intent(in) :: b
            double precision dimension(n), intent(out) :: c
        end subroutine add_vectors
    end interface
end python module fastcalc

ビルドします。

python -m numpy.f2py -c fastcalc.pyf add_vectors.f90

.pyfにモジュール名が定義されているので、この場合は-mを付けません。既存のソルバーにPython用の入口を設けるときには、この方法が便利です。

8. Fortran 2018の機能を使うときの注意点

Fortran 2018で書いたコードをgfortranがコンパイルできることと、f2pyがその構文を理解してPython向けのインターフェースを自動生成できることは別です。

機能 考え方
通常のサブルーチン・関数 f2pyで扱いやすい構成です
数値配列・intent 主な用途として対応しています
do concurrent Fortran側の計算内部で使えます
real(real64) KINDの型マッピングを確認します
allocatable配列 公開の仕方に制約があります
派生型typeや多態性class そのまま自動公開できない場合があります
submodule、Coarray 計算内部で使い、境界は単純なAPIにする方が確実です

real(real64)の型指定

Fortran 2018では、次のように精度を明示する書き方を使えます。

use, intrinsic :: iso_fortran_env, only: real64
real(real64) :: water_level

f2py側でこのKINDを意図どおりに認識できない場合は、.f2py_f2cmapファイルを用意します。

dict(real=dict(real64="double"))

この対応付けはNumPyの公式ドキュメントにも掲載されています。複雑な派生型まで自動で扱えるようになる仕組みではないため、公開APIには単純な実数・整数・配列を使う方が管理しやすいです。

ソルバー本体と公開APIを分離する

例えば河川計算ソルバーなら、次のような構成です。

hydraulic_solver/
├── src/
│   ├── continuity.f90
│   ├── momentum.f90
│   ├── boundary.f90
│   └── hydraulic_api.f90
├── interface/
│   └── hydraulic_api.pyf
└── python/
    ├── run_solver.py
    └── test_solver.py

Fortran内部ではモジュールや派生型を利用しても構いません。Pythonとの境界にあるhydraulic_api.f90では、配列と基本データ型だけを受け渡すようにします。既存の計算コアを大きく変更せずに済み、Fortran単独での試験も続けられます。

9. 高速化するときの考え方

f2pyはPythonとFortranを接続するツールであり、数値計算そのものの高速化はFortran側の実装やコンパイルオプションが担当します。

コンパイラの最適化

例えば、gfortranで最適化を有効にしてビルドするなら、次のように指定できます。

FC=gfortran \
FFLAGS="-O3 -march=native -std=f2018" \
python -m numpy.f2py -c -m fastcalc \
    add_vectors.f90 diffusion.f90

-O3は積極的な最適化、-march=nativeはビルドするCPUに合わせた命令の利用です。-march=nativeで作ったバイナリを別のCPUへ配布する場合は注意してください。

ベクトル化の状況を確認したいときは、-fopt-info-vec-optimizedによるコンパイラの診断メッセージも参考になります。

水理計算では、計算速度だけでなく水収支や再現性が重要です。浮動小数点演算の意味が変わり得る-ffast-mathは、検証が済むまで安易に使わない方がよいと思います。

SIMDとOpenMP

配列の連続領域を更新するループはSIMD化の候補になります。また、セル同士が独立している処理はOpenMPによる並列化も検討できます。

!$omp parallel do private(i,j) schedule(static)
do j = 1, ny
    do i = 1, nx
        hnew(i,j) = h(i,j) + source(i,j)
    end do
end do
!$omp end parallel do

OpenMPを使う場合は、コンパイル時だけでなくリンク時にも必要な設定を行います。複数のFortranファイルや外部ライブラリがある場合は、Meson側で依存関係を明示した構成にすると管理しやすくなります。

ただし、共有変数やsave変数を扱うコードは、並列実行時のデータ競合に注意が必要です。

PythonとFortranを何度も往復しない

大量のセルを計算する場合、例えば次のような違いがあります。

避けたい例
Python: 1セル分をFortranに渡す → 戻す → 次のセルを渡す → …

改善例
Python: メッシュ全体を渡す
Fortran: 全セル・複数時間ステップをまとめて計算
Python: 必要な結果を受け取る

高速化を考えるときは、Fortranのループだけでなく、配列のコピー量とPythonからの呼び出し回数も確認する必要があります。

10. Dockerで実行環境を揃える

同じ構成をDockerでも利用できます。以下は配列加算と2次元拡散をコンテナ内でビルドする例です。

Dockerfile（クリックして開く） 
FROM python:3.12-slim

# C / FortranコンパイラとMesonの実行に必要なツール
RUN apt-get update && \
    apt-get install -y --no-install-recommends \
        build-essential \
        gfortran \
        ninja-build \
        pkg-config && \
    rm -rf /var/lib/apt/lists/*

# f2py本体はNumPyに含まれる
RUN python -m pip install --no-cache-dir \
    numpy==2.5.4 \
    meson \
    pytest

WORKDIR /work

COPY add_vectors.f90 diffusion.f90 test_diffusion.py ./

RUN python -m numpy.f2py -c -m fastcalc \
    add_vectors.f90 diffusion.f90

CMD ["python", "test_diffusion.py"]

docker build -t f2py-demo .
docker run --rm f2py-demo

Dockerfileは実装例です。Pythonパッケージのバージョンだけでなく、コンパイラや依存ライブラリも固定したい場合は、ベースイメージのダイジェストやパッケージのバージョンも管理します。

11. エラーが出たときに確認すること

症状 確認すること
No module named numpy 仮想環境が有効か、NumPyが入っているか
No module named mesonbuild python -m pip install mesonを実行したか
Fortranコンパイラが見つからない gfortran --versionとFCの設定
No module named fastcalc .soの生成場所と実行するPython環境
配列の型・次元に関するエラー dtype、shape、order、引数の順序
undefined symbol リンク対象のオブジェクトやライブラリ
Fortranの構文解析に失敗する .pyfでAPIを限定できないか

f2pyの公開関数を確認するには、次のコマンドが便利です。

python -c "import fastcalc; print(fastcalc.__doc__)"

また、NumPy 2.xでは以前のnumpy.f2py.compile()が削除されています。Pythonスクリプトからビルドを呼ぶ場合は、現在のコマンドラインツールをsubprocess.run()などで起動する形に変更します。

12. 動作確認で大切にしたいこと

今回の2つのFortranプログラムは、gfortran 14.2.0でFortran 2018としてコンパイルし、配列加算と2次元拡散の計算値を確認しています。また、NumPy 2.3.5のf2pyで.pyfを生成・解析できることも確認しています。

一方、NumPy 2.5.4とMesonを使った拡張モジュールの最終ビルド、およびPythonからの実行までは、この検証環境では実施していません。掲載したDockerfileについても、コンテナをビルドした結果を示すものではありません。環境を構築したら、記事中のテストを実行して確認してください。

河川・氾濫計算のようなプログラムでは、結果が返ってくることだけでなく、次の項目も確認しておくべきです。

  • Fortran単独実行とf2py経由で、同じ入力に対して計算結果が一致するか
  • 配列の形状、座標軸、メモリ配置が想定どおりか
  • 水位、流量、水収支などの物理量に不自然な差が生じていないか
  • 配列コピーやPythonとの往復に、どの程度の時間を使っているか
  • OpenMPや最適化オプションを変えても、許容できる計算精度が保たれるか

おわりに

f2pyを使うと、Fortranで作った数値計算プログラムをPythonから比較的少ない変更で利用できます。

特に既存の計算ソルバーがある場合、まずは簡単なサブルーチンを公開して配列の受け渡しを確認し、次に時間積分などの計算処理をまとまった単位で呼び出す形にすると進めやすくなります。

Fortran 2018の機能を全面的にf2pyへ見せる必要はありません。計算コアとPython向けAPIを分けておけば、Fortran側の実装を維持しながら、Pythonのデータ処理、可視化、Web APIなどと組み合わせることができます。

次の段階では、1次元河道や2次元メッシュの時間積分をFortran側にまとめ、Fortran単独実行とf2py経由で計算結果・水収支・処理時間を比較してみたいと思います。

参考資料

1
0
0

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
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?