🎃

pybind11を使ってC++/FortranプログラムをPythonで動かす <その4>

に公開
C++
Fortran
pybind11
tech

はじめに

Python ➔ C++ ➔ Fortranの順にバインディングを行うことで、PythonからFortranプログラムを動かすテストの第4回です。入力パラメータがPython ➔ C++ ➔ Fortranの順に流れ、計算結果がFortran ➔ C++ ➔ Pythonの順に返ってくるような処理をイメージしています。

プログラム例 7

こんどは多次元のNumPy配列を転送してみます。
2次正方行列を入力パラメータに設定し、その逆行列をFortranプログラムで計算して求めます。

Fortranソースコード

matinv2.f90
module matinv2_mod
    use, intrinsic :: iso_c_binding, only: c_double
    use, intrinsic :: iso_fortran_env, only: real64
    implicit none
contains
    ! Inverse of a 2x2 matrix
    subroutine matinv2(a, a_inv) bind(C)
        real(c_double), intent(in) :: a(2, 2)
        real(c_double), intent(out) :: a_inv(2, 2)
        real(real64) :: detinv

        detinv = 1. / (a(1, 1) * a(2, 2) - a(1, 2) * a(2, 1))
        a_inv(1, 1) =  detinv * a(2, 2)
        a_inv(2, 1) = -detinv * a(2, 1)
        a_inv(1, 2) = -detinv * a(1, 2)
        a_inv(2, 2) =  detinv * a(1, 1)
    end subroutine
end module

matinv2は関数にしたくなりますが、返り値が配列であるような関数にはbind(C)属性をつけられないので、サブルーチンにして引数で計算結果を渡すことにしました。

C++ソースコード

example.cpp
#include <pybind11/pybind11.h>
#include <pybind11/numpy.h>

// Using Fortran code within C++
extern "C" {
    void matinv2(const double *a, double *a_inv);
}

namespace py = pybind11;

// Matrix inverse
py::array_t<double, py::array::f_style> c_matinv2(const py::array_t<double, py::array::f_style> &a) {
    py::buffer_info buf1 = a.request();

    // Allocate the buffer
    auto a_inv = py::array_t<double, py::array::f_style>(buf1.shape);

    py::buffer_info buf2 = a_inv.request();

    const double *a_ptr = static_cast<const double *>(buf1.ptr);
    double *a_inv_ptr = static_cast<double *>(buf2.ptr);

    matinv2(a_ptr, a_inv_ptr);

    return a_inv;
}

// Bindings with Python
PYBIND11_MODULE(example, m) {
    m.def("matinv2", &c_matinv2);
}

ポイントは、

  • py::array_tの第2テンプレート引数の py::array::f_style

です。これでFortranと同じ”列優先”（Column Major）で配列を解釈するように指示しています。Pythonからテストしてみます。

>>> import example
>>> import numpy as np
>>> a = np.arange(1, 5, dtype=float).reshape(2, 2)
>>> a
array([[1., 2.],
       [3., 4.]])
>>> example.matinv2(a)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])
>>> np.linalg.inv(a)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

Pythonの線形代数ライブラリ（linalg）で計算した結果と同じ結果が出ました。

参考

https://pybind11.readthedocs.io/en/stable/

Discussion