SymPy で記号演算した結果の Fortran からの利用

Fortran から Julia を呼んで、そこで PyCall + SymPy して記号処理します。出てきた結果を Julia の関数として定義して、それを Fortran から呼び出します。

実行結果

まず julia で普通に f(x) = x² と定義した関数が Fortran から呼び出せることを確かめます。次に SymPy で log(x)/x を積分してやると log(x)²/2 と求まるので、これを Julia の関数に直します。そうして Fortran から呼び出します。念のため Fortran 側で log(x)²/2 を計算して結果を比較しています。

hp8@HP8:~/forjulia$ ifort -o test  -L$JULIA_DIR/lib -Wl,-rpath,$JULIA_DIR/lib sym.f90 -ljulia
hp8@HP8:~/forjulia$ ./test
 call SymPy from Julia from Fortran
 f(x)=x^2
 f(2.0d0)=   4.00000000000000

 integrate(log(x) / x, x)= log(x)^2/2
  1   1.0000000000000000   0.0000000000000000   0.0000000000000000
  2   2.0000000000000000   0.2402265069591007   0.2402265069591007
  3   3.0000000000000000   0.6034744804062910   0.6034744804062910
  4   4.0000000000000000   0.9609060278364028   0.9609060278364028
  5   5.0000000000000000   1.2951451969901173   1.2951451969901173
  6   6.0000000000000000   1.6052009977842006   1.6052009977842006
  7   7.0000000000000000   1.8932831540982358   1.8932831540982358
  8   8.0000000000000000   2.1620385626319059   2.1620385626319059
  9   9.0000000000000000   2.4138979216251641   2.4138979216251641
 10  10.0000000000000000   2.6509490552391997   2.6509490552391997

今回は ifort でコンパイルしてみました。gfortran より型チェックにうるさいです。

プログラム

SymPy の出力は直では普通の関数には出来ないようなので、一度文字列化してパースして eval して Fortran から呼び出せる普通の関数にしています。この辺は詳しいことがよく分かっていないので素朴にやっています。

eval(Meta.parse( "f(x)=" * string(integrate(y, x)) ))

module julia_mod
    use, intrinsic :: iso_fortran_env
    use, intrinsic :: iso_c_binding
    implicit none

!
! global variables
!
    type(c_ptr) :: jl_base_module
    bind(c, name = 'jl_base_module') :: jl_base_module

    type(c_ptr) :: jl_int64_type
    bind(c, name = 'jl_int64_type') :: jl_int64_type

    type(c_ptr) :: jl_float64_type
    bind(c, name = 'jl_float64_type') :: jl_float64_type
!
! interfaces
!
    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

        type(c_ptr) function jl_eval_string(pstr) bind(c, name = 'jl_eval_string')
            import
            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

        type(c_ptr) function jl_box_float64(x) bind(c, name = 'jl_box_float64')
            import
            real(real64), value :: x
        end function jl_box_float64

        type(c_funptr) function jl_get_global(modu, psymbol) bind(c, name = 'jl_get_global')
            import
            type(c_ptr), value :: modu
            type(c_ptr), value :: psymbol
        end function jl_get_global

        type(c_ptr) function jl_symbol(ptext) bind(c, name = 'jl_symbol')
            import
            type(c_ptr), value :: ptext
        end function jl_symbol

        type(c_ptr) function jl_call1(pfunc, parg) bind(c, name = 'jl_call1')
            import
            type(c_funptr), value :: pfunc
            type(c_ptr), value :: parg
        end function jl_call1

    end interface

contains

    subroutine eval(text, iret)
        character(*), intent(in) :: text
        type(c_ptr), intent(out), optional :: iret
        character(:), allocatable, target :: buff
        type(c_ptr) :: kret
        buff = text // achar(0)
        kret = jl_eval_string(c_loc(buff))
       if (present(iret)) iret = kret
    end subroutine eval

    function jl_get_function(modu, namen)
        type(c_ptr), value :: modu
        character(*), intent(in) :: namen
        type(c_funptr) :: jl_get_function
        character(:), allocatable, target :: text
        text = namen // achar(0)
        jl_get_function = jl_get_global(modu, jl_symbol(c_loc(text)))
    end function jl_get_function

end module julia_mod

program sympy
    use julia_mod
    implicit none
    type(c_ptr) :: iret, arg1, iatyp
    type(c_funptr) :: ifun
    real(real64) :: x, y
    integer :: i

    call jl_init()

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

    call eval('println(" f(x)=x^2")')
    call eval('f(x)=x^2', iret)
    ifun = transfer(iret, ifun) ! c_ptr -> c_funptr
    arg1 = jl_box_float64(2.0_real64)
    iret = jl_call1(ifun, arg1)
    x = jl_unbox_float64(iret)
    print *, 'f(2.0d0)=', x
    print *

    call eval('using PyCall')
    call eval('using SymPy') ! or SymEngine
    call eval('@vars x y')
    call eval('y = log(x) / x')
    call eval('print(" integrate(log(x) / x, x)= ")')
    call eval('println(integrate(y, x))')
    call eval('eval(Meta.parse( "f(x)=" * string(integrate(y, x)) ))', iret)
    ifun = transfer(iret, ifun) ! c_ptr -> c_funptr

    do i = 1, 10
        x = real(i, real64)
        arg1 = jl_box_float64(x)
        iret = jl_call1(ifun, arg1)
        y = jl_unbox_float64(iret)
        print '(i3, 3(1x, f20.16))', i, x, y, log(x)**2 / 2
    end do

    call jl_atexit_hook(0)
end program sympy