Embedding Julia in Fortran

スタバのソファーでふんぞり返って、タラタラとネットサーフィンｗしていたら C から Julia を呼び出す本家のページにたどり着きましたので眺めやると、それほど面倒くさそうでもなさそうなので Fortran から呼び出すことに挑戦してみました。それにつけてもモバイルオーダーまじ楽ちん。

とりあえず、本家のここのページを参考にして、C 用のルーチンのインターフェースを最小限用意して実行してみました。Windows だと色々めんどそうなので WSL + gfortran でやってみることにします。

docs.julialang.org

Fortran ソース・プログラム

基本的にC 用のサンプルを丸写し的にやっています。

文字列の内容を Julia で実行させています。
Julia 側の演算結果の倍精度実数のスカラー変数を Fortran 側にもってきています。
using LinearAlgebra で Julia で逆行列を求めています。(Fortran 側には持ってきていません）
using Plots で~~グラフ表示は出来ませんでした~~。display 命令で表示できることが分かりました。ただし図を表示したままにするには Fortran の puase 文（規格から廃止になりましたが）などを使う必要があります。~~しかし~~図は savefig でファイルに保存できました。

module julia_mod
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: iso_c_binding
    implicit none
    interface
        subroutine jl_init() bind(c, name = 'jl_init__threading')
        end subroutine jl_init

        subroutine jl_atexit_hook(status) bind(c, name = 'jl_atexit_hook')
            import
            integer(c_int), value :: status
        end subroutine jl_atexit_hook

        function jl_eval_string(pstr) bind(c, name = 'jl_eval_string')
            import
            type(c_ptr) :: jl_eval_string
            type(c_ptr), value :: pstr
        end function jl_eval_string

        real(real64) function jl_unbox_float64(ptr) bind(c, name = 'jl_unbox_float64')
            import
            type(c_ptr), value :: ptr
        end function jl_unbox_float64
    end interface

contains
    subroutine eval(text)
        character(*), intent(in) :: text
        type(c_ptr) :: iret
        character(:), allocatable, target :: buff
        buff = text // achar(0)
        iret = jl_eval_string(c_loc(buff))
    end subroutine eval
end module julia_mod

program test
    use julia_mod
    implicit none
    character(:), allocatable, target :: text
    type(c_ptr) :: iret
    real(real64) :: x

    call jl_init()

    call eval('println("Julia from Fortran")')

    call eval('x=1/6')
    call eval('println(sinpi(x))')

    text = 'sqrt(2.0)'//achar(0)
    iret = jl_eval_string(c_loc(text))
    x = jl_unbox_float64(iret)
    print *, 'SQRT(2.0)=', x, sqrt(2.0_real64)

    call eval('using LinearAlgebra')
    call eval('A = [1 2 3; 4 1 6; 7 8 1]')
    call eval('println(A)')
    call eval('println("det(A)=", det(A))')
    call eval('println("A^-1=", inv(A))')

    call eval('using Plots')
    call eval('gui')
    call eval('p = plot(sin)')
    call eval('plot!(cos)')
    call eval('savefig("test.png")')

    call eval('display(p)')
    pause
    call eval('display(plot(rand(10), st=:pie))')
    pause

    call jl_atexit_hook(0)
end program test

実行結果

gfortran の pause 文はうるさいので ifort でやり直しました。

hp8@HP8:~$ export JULIA_DIR=~/julia-1.5.2
hp8@HP8:~/forjulia$ gfortran -o test -fPIC -L$JULIA_DIR/lib -Wl,-rpath,$JULIA_DIR/lib julia.f90 -ljulia
hp8@HP8:~/forjulia$ ifort    -o test       -L$JULIA_DIR/lib -Wl,-rpath,$JULIA_DIR/lib julia.f90 -ljulia
hp8@HP8:~/forjulia$ ./test
Julia from Fortran
0.5
 SQRT(2.0)=   1.41421356237310        1.41421356237310
[1 2 3; 4 1 6; 7 8 1]
det(A)=104.0
A^-1=[-0.4519230769230769 0.2115384615384615 0.08653846153846155; 0.3653846153846154 -0.1923076923076923 0.057692307692307675; 0.24038461538461536 0.057692307692307696 -0.0673076923076923]
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-hp8'
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-hp8'
FORTRAN PAUSE
PAUSE prompt>
FORTRAN PAUSE
PAUSE prompt>
hp8@HP8:~$ eog test.png

f:id:fortran66:20201109133136p:plain f:id:fortran66:20201107022314p:plain