15 Julia 与 Fortran 的互操作

15.1 在 Julia 中调用 Fortran

在 Julia 中调用 Fortran 与调用 C 几乎完全相同。作为示例，我们在 interoperability/fortran 目录下提供了一个简单的动态链接库 libJuliaFortran，其中包含了几个简单的 Fortran 函数。

interoperability/fortran/src/lib.f90

module Julia
    integer :: last = 0;

    contains

subroutine add(a,b,c)
  implicit none
  integer, intent(in) :: a,b
  integer, intent(out) :: c
  c = a + b
  last = c
end subroutine add

function fadd(a,b) result(c)
  implicit none
  integer, intent(in) :: a,b
  integer :: c
  c = a + b
end function fadd

subroutine vadd(a,b,c,n)
  implicit none
  integer, intent(in) :: n, a(n), b(n)
  integer, intent(out) :: c(n)
  c = a + b
end subroutine vadd

end module

接下来我们尝试在 Julia 中调用这些函数。

using Libdl: dlopen, dlsym, dlclose

dlopen("libJuliaFortran") do libJuliaFortran
    # 使用 `ccall` 调用 Fortran 中的函数
    add = dlsym(libJuliaFortran, :__julia_MOD_add)
    ans = Ref(Int32(0))
    ccall(add, Cvoid, (Ref{Cint}, Ref{Cint}, Ref{Cint}), 1, 2, ans)
    @show ans[]

    # 使用 `cglobal` 获取 Fortran 中的全局变量
    last_sym = dlsym(libJuliaFortran, :__julia_MOD_last)
    last_ptr = cglobal(last_sym, Cint)
    last = unsafe_load(last_ptr)
    @show last

    dlclose(libJuliaFortran)
    return nothing
end

ans[] = 3
last = 3

说明

调用 Fortran 函数时，传入的标量类型必须是引用类型，而不是值类型。例如，Int32 不能被传入，而必须使用 Ref{Int32}。所以，在上面的例子中，我们将 ans 定义为 Ref(Int32(0)) 而不是直接定义为 Int32(0)，否则我们无法得到正确的结果。

读者可能会注意到我们并没有对 1 或 2 进行包装，这是因为 Julia 会自动对它们进行包装，因此我们可以直接将它们传入 Fortran 函数。但对于我们需要获取的 ans 来说，如果由 Julia 自动包装，Fortran 将会对包装后的指针进行修改，而 Julia 这边 ans 的值将不会发生变化，就像下面这样：

dlopen("libJuliaFortran") do libJuliaFortran
    add = dlsym(libJuliaFortran, :__julia_MOD_add)
    ans = 0
    ccall(add, Cvoid, (Ref{Cint}, Ref{Cint}, Ref{Cint}), 1, 2, ans)
    @show ans
    dlclose(libJuliaFortran)

    return nothing
end

ans = 0

如果调用的是 Fortran 中的 function，则不需要将函数返回值定义为引用类型，就像下面这样：

dlopen("libJuliaFortran") do libJuliaFortran
    add = dlsym(libJuliaFortran, :__julia_MOD_fadd)
    ans = ccall(add, Cint, (Ref{Cint}, Ref{Cint}), 1, 2)
    @show ans
    dlclose(libJuliaFortran)

    return nothing
end

ans = 3

如果传入的参数是数组，则应当将对应的参数类型声明为 Ptr{T}，就像下面这样：

dlopen("libJuliaFortran") do libJuliaFortran
    add = dlsym(libJuliaFortran, :__julia_MOD_vadd)
    a = Int32[1, 2, 3, 4, 5]
    b = Int32[5, 4, 3, 2, 1]
    ans = zeros(Int32, 5)
    ccall(add, Cvoid, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ref{Cint}), a, b, ans, length(a))
    @show ans
    dlclose(libJuliaFortran)

    return nothing
end

ans = Int32[6, 6, 6, 6, 6]

说明

目前，在 Julia 中调用 Fortran 时，还存在一些缺陷：

函数调用时必须显式传入数组的长度，而不能使用 a(:) 这样的方式动态获取数组长度。
必须传入 Fortran 中定义的所有参数，即使某些参数在 Fortran 中被声明为 optional。

15.2 在 Fortran 中调用 Julia

利用 Fortran 的 ISO_C_BINDING 模块，我们可以在 Fortran 中调用 Julia 的 C API，进而调用 Julia。在 interoperability/fortran 目录下，我们提供了一个在 Fortran 中调用 Julia 函数的例子。项目仓库中已经使用 CMake 配置好了编译和链接的有关选项，直接构建即可生成可执行程序 FortranJulia。

interoperability/fortran/src/main.f90

module JuliaWrapper
   use, intrinsic :: iso_c_binding, only : c_ptr
   type(c_ptr), bind(c, name="jl_main_module") :: jl_main_module

   interface
      subroutine jl_init() bind(c, name="jl_init")
      end subroutine jl_init

      type(c_ptr) function jl_eval_string(str) bind(c, name="jl_eval_string")
         use, intrinsic :: iso_c_binding, only : c_char, c_ptr
         character(kind=c_char), dimension(*), intent(in) :: str
      end function jl_eval_string

      type(c_ptr) function jl_symbol(str) bind(c)
         use, intrinsic :: iso_c_binding, only : c_char, c_ptr
         character(kind=c_char), dimension(*), intent(in) :: str
      end function jl_symbol

      type(c_ptr) function jl_get_global(m, var) bind(c, name="jl_get_global")
         use, intrinsic :: iso_c_binding, only : c_ptr
         type(c_ptr), value, intent(in) :: m, var
      end function jl_get_global

      type(c_ptr) function jl_box_float64(x) bind(c, name="jl_box_float64")
         use, intrinsic :: iso_c_binding, only : c_double, c_ptr
         real(c_double), value, intent(in) :: x
      end function jl_box_float64

      real(c_double) function jl_unbox_float64(x) bind(c, name="jl_unbox_float64")
         use, intrinsic :: iso_c_binding, only : c_double, c_ptr
         type(c_ptr), value, intent(in) :: x
      end function jl_unbox_float64

      type(c_ptr) function jl_call1(f, arg1) bind(c, name="jl_call1")
         use, intrinsic :: iso_c_binding, only : c_ptr
         type(c_ptr), value, intent(in) :: f, arg1
      end function jl_call1

      subroutine jl_atexit_hook(code) bind(c, name="jl_atexit_hook")
         use, intrinsic :: iso_c_binding, only : c_int
         integer(c_int), value :: code
      end subroutine jl_atexit_hook
   end interface
end module JuliaWrapper

program FortranJulia
   use JuliaWrapper
   use iso_c_binding, only : c_char, c_null_char
   implicit none

   type(c_ptr) :: res, sym, jl_sqrt

   call jl_init()
   sym = jl_symbol(c_char_"sqrt"//c_null_char)
   jl_sqrt = jl_get_global(jl_main_module, sym)
   res = jl_call1(jl_sqrt, jl_box_float64(2.0d0))
   print *, jl_unbox_float64(res)
   call jl_atexit_hook(0)
end program FortranJulia

运行 FortranJulia，输出结果应当为：

   1.4142135623730951

说明

在上面的例子中，需要特别注意的有以下两点：

JuliaWrapper 模块中定义的辅助函数中，需要为非数组的输入参数声明 value 属性。
因为 C 中的字符串是以 NULL 结尾的，所以在对应函数需要字符串时，我们需要将字符串写成 c_char_"sqrt"//c_null_char 这样的形式。