Calling a Fortran function from Julia, returning an array: unknown function, segfault?

Question

I want to call functions in my Fortran library from Julia. In this case, I have a function eye that takes an Integer, and returns a two-dimensional array of integers.

The Fortran module is compiled into a shared library using

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

And thereafter I am attempting to call it from the interactive Julia interpreter like that (name obtained from nm):

julia> n=5
5

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

This, however, immediately results in Julia throwing a segfault at me:

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

Am I calling the function the wrong way? What is the correct name of the function? (It doesn't appear to be just eye, as that doesn't work either.)

As an additional question: does Julia do anything with the memory-orientation of the resulting arrays? Fortran and Julia are both column-major, but I wonder if due to ccall() Julia might think it should tranpose them?

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

Answer 1

There a couple issues with your approach. Returning an array directly to julia is problematic because Fortran arrays are not interoperable with C unless specific conditions are met. When you make the array interoperable (add bind(C) to your procedure and give the array a C type) the compiler (gfortran) will complain:

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

To get around this we can return the array through a dummy argument. You'll want to make this an intent(inout) argument and construct the array in julia to avoid any memory/scope issues with creating the array in Fortran.

Secondly, the optional argument is problematic and skimming the Julia docs I'm not sure it is even supported. Note that not even Fortran can call Fortran with optional arguments without an explicit interface and as Julia doesn't interact with the .mod file and seems to expect the C way of doing things, it probably won't work (and Fortran 2008 15.3.7 p2.6 seems to say it isn't supported). There are workarounds though -- you can create multiple Fortran procedures with varying numbers of arguments and then call the procedure with optional arguments from them.

First, consider this Fortran module, which started with your example but is pared down to just what is necessary to demonstrate the interop:

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

Note that I've moved um to being an inout dummy argument and since we're not returning a value changed the procedure to a subroutine. I have also removed the optional argument. I have also used C interop types and bound a C name to the procedure. You can compile this as in your question.

In Julia you can now do the following:

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

Note that we can call the function as just :eye since we used the bind(C,name="eye") C interop in our Fortran.

And lastly, if we change the do loop in my Fortran example to be um(i,j) = i*10+j, we can see that no transposition happens in the array:

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

---

The particular reason for your segfault could have been a number of things -- mismatched data types, issues with the return type, issues with the optional argument or mismatch in the actual call and variables passed.