从 Julia 调用 Fortran 函数，返回数组：未知函数，段错误？

我想从 Julia 调用 Fortran 库中的函数。在这种情况下，我有一个函数eye它接受一个 Integer，并返回一个二维整数数组。

Fortran 模块被编译成共享库，使用

$ gfortran -shared -fPIC -o matrix_routines.so matrix_routines.f90

此后我尝试从交互式 Julia 解释器中调用它（名称从nm):

julia> n=5
5

julia> ccall( (:__matrix_routines_MOD_eye, "/path/to/library/matrix_routines.so"), Array{Int64,2} , (Ptr{Int64},), &n )

然而，这立即导致 Julia 向我抛出一个段错误：

signal (11): Segmentation fault
__matrix_routines_MOD_eye at /path/to/library/matrix_routines.so (unknown line)
anonymous at no file:0
unknown function (ip: -1137818532)
jl_f_top_eval at /usr/bin/../lib/julia/libjulia.so (unknown line)
eval_user_input at REPL.jl:53
jlcall_eval_user_input_19998 at  (unknown line)
jl_apply_generic at /usr/bin/../lib/julia/libjulia.so (unknown line)
anonymous at task.jl:95
jl_handle_stack_switch at /usr/bin/../lib/julia/libjulia.so (unknown line)
julia_trampoline at /usr/bin/../lib/julia/libjulia.so (unknown line)
unknown function (ip: 4199613)
__libc_start_main at /usr/bin/../lib/libc.so.6 (unknown line)
unknown function (ip: 4199667)
unknown function (ip: 0)
zsh: segmentation fault (core dumped)  julia

我调用该函数的方式是否错误？该函数的正确名称是什么？ （看来不仅仅是eye，因为这也不起作用。）

另外一个问题：Julia 对结果数组的内存方向做了什么处理吗？ Fortran 和 Julia 都是列优先的，但我想知道由于 ccall() Julia 是否可能认为应该转置它们？

module matrix_routines
    implicit none

    private

    public :: eye

    contains

        pure function eye(n,offset) result(um) !{{{
            integer, intent(in) :: n
            integer, intent(in), optional :: offset

            integer, dimension(n,n) :: um

            integer :: i, l, u, os

            um = 0

            l = 1
            u = n
            os = 0

            if (present(offset)) then
                os = offset
            end if

            if (abs(os) < n) then
                if (os > 0) then
                    u = n - os
                else if (os < 0) then
                    l = 1 - os
                end if

                do i=l, u
                    um(i, i+os) = 1
                end do
            end if

        end function eye !}}}
end module matrix_routines

您的方法有几个问题。直接将数组返回给 julia 是有问题的，因为除非满足特定条件，否则 Fortran 数组不能与 C 互操作。当您使数组可互操作时（添加bind(C)给你的程序并给数组一个C类型）编译器（gfortran）会抱怨：

Error: Return type of BIND(C) function 'um' at (1) cannot be an array

为了解决这个问题，我们可以通过虚拟参数返回数组。你会想把它变成一个intent(inout)参数并在 julia 中构造数组，以避免在 Fortran 中创建数组时出现任何内存/范围问题。

其次，可选参数是有问题的，浏览 Julia 文档我不确定它是否受支持。请注意，即使 Fortran 也无法在没有显式接口的情况下使用可选参数调用 Fortran，并且 Julia 不与.mod文件并且似乎期望 C 的处理方式，它可能不会工作（Fortran 2008 15.3.7 p2.6 似乎说它不受支持）。不过，也有解决方法——您可以创建多个具有不同数量参数的 Fortran 过程，然后使用其中的可选参数调用该过程。

首先，考虑这个 Fortran 模块，它以您的示例开始，但被简化为演示互操作所需的内容：

module matrix_routines
  implicit none

  private
  public :: eye

contains

  pure subroutine eye(n,um) bind(C,name="eye") !{{{
    use, intrinsic :: iso_c_binding, only: c_int
    implicit none
    integer(c_int), intent(in) :: n
    integer(c_int), intent(inout), dimension(n,n) :: um

    integer :: i, j

    do j=1,n
       do i=1,n
          um(i,j) = i+j
       end do
    end do

  end subroutine eye !}}}
end module matrix_routines

请注意，我已经搬家了um成为一个inout虚拟参数，由于我们不返回值，因此将过程更改为子例程。我还删除了可选参数。我还使用了 C 互操作类型并将 C 名称绑定到过程。您可以按照您的问题进行编译。

在 Julia 中，您现在可以执行以下操作：

julia> n = 2
2

julia> um = zeros(Int32, n, n)
2x2 Array{Int32,2}:
 0  0
 0  0

julia> ccall((:eye, "matrix_routines.so"), Void, (Ptr{Int32}, Ptr{Array{Int32,2}}), &n, um)

julia> um
2x2 Array{Int32,2}:
 2  3
 3  4

julia> n = 4
4

julia> um = zeros(Int32, n, n)
4x4 Array{Int32,2}:
 0  0  0  0
 0  0  0  0
 0  0  0  0
 0  0  0  0

julia> ccall((:eye, "matrix_routines.so"), Void, (Ptr{Int32}, Ptr{Array{Int32,2}}), &n, um)

julia> um
4x4 Array{Int32,2}:
 2  3  4  5
 3  4  5  6
 4  5  6  7
 5  6  7  8

请注意，我们可以将该函数称为:eye因为我们使用了bind(C,name="eye")Fortran 中的 C 互操作。

最后，如果我们将 Fortran 示例中的 do 循环更改为um(i,j) = i*10+j，我们可以看到数组中没有发生转置：

julia> um
3x3 Array{Int32,2}:
 11  12  13
 21  22  23
 31  32  33

段错误的具体原因可能有很多原因——数据类型不匹配、返回类型问题、可选参数问题或实际调用和传递的变量不匹配。