Create interface to fortran modules with ccall

New to Julia
,
1

Hi all,

the julia documentation provides quite some example on how to call C functions. However, I want to accress some fortran modules and currently I’m a bit lost. Say I have a fortran module like:

module foo
    use :: iso_fortran_env
    implicit none
contains

    subroutine foo1(n,x,y)
        integer(int64), intent(in) :: n
        integer(int64), intent(in),  dimension(n) :: x
        integer(int64), intent(out), dimension(n) :: y

        y = 2*x
    end subroutine foo1

    subroutine foo2(x,y)
        integer(int64), intent(in),  dimension(:) :: x
        integer(int64), intent(out), dimension(size(x)) :: y

        y = 2*x
    end subroutine foo2

    function bar1(n,x) result(y)
        integer(int64), intent(in) :: n
        integer(int64), intent(in), dimension(n) :: x

        integer(int64), dimension(n) :: y

        y = 2*x
    end function bar1

end module foo

which I compile (on Ubuntu) with

gfortran -Wall -Wextra -fPIC -c mod_foo.F08
gfortran -shared -Wl,-soname,libfoo.so -o libfoo.so mod_foo.o

then, within julia the following works:

julia> t = ones(Int64,2); ccall((:__foo_MOD_foo1, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ptr{Int64}), 2, [3, -1], t); t
2-element Array{Int64,1}:
  6
 -2

but any attempt to call foo2 or bar1 has ended in segmentation faults. Can anybody tell me how to access these?

Regards,
Laurent

1 Like
2

You don’t state how you’re trying to do that.

I’ve never used Fortran… (or the C API of Julia…) just curious how this can work. I’m not sure foo2 can work. With foo1 you’re sending the length n, as you would do with C. I’m not familiar with how size(x) works in Fortran, but it seems to me Julia must send the length somehow and that wouldn’t work.

function bar1 is superficially similar to soubroutine foo1, still not sure how an array is returned.

3

I tried the following to call bar1:

q = ccall((:__foo_MOD_bar1, "libfoo.so"), Ptr{Int64}, (Ref{Int64}, Ref{Int64}), 2, [3, 5])

which results in a segmentation fault. Similarly:

q = ones(Int64,2); q = ccall((:__foo_MOD_foo2, "libfoo.so"), Void, (Ref{Int64},Ref{Int64}), [3, 5], q); q

doesn’t crash, but won’t contain anything. On the other hand:

q = ones(Int64,2); q = ccall((:__foo_MOD_foo2, "libfoo.so"), Void, (Ref{Int64},Ptr{Int64}), [3, 5], q); q

yields a segmentation fault.

4

Julia doesn’t support “assumed shape” Fortran arrays. Theoretically you could match the ABI gfort uses for array descriptors by hand (it’s just a struct, AFAICT) but probably much easier to write size-specifed wrappers.

Note also that you should declare the subroutine bind(C).

1 Like
5

Fortran is tricky. This is in particular due to a Fortran compiler being permitted to do more or less whatever it wants to speed up computation. The only reliable way is indeed to go through the Iso C Binding features provided by newer Fortran standards.

However, the following two links might give some rather detailed picture of how to handle this: Maclaren (2014) and Pletzer et. al. (2008).

Essentially, you’ll have to either change the original Fortran code or write wrappers in Fortran to do the translation of types. Something like:

use, intrinsic :: iso_c_binding

subroutine foo_wrapper(n,x,y)    bind(C, name=foo)
    integer(kind=c_long), intent(in)  :: x(n), n
    integer(kind=c_long), intent(out) :: y(n)
    call foo1(n, x, y)
end subroutine foo_wrapper

Then it is guaranteed to be C-compliant, portable etc.

And you can directly use the specified name= as the function name, as in

ccall( (:foo, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ref{Int64}), ...)
#      ^^^^^^

All other versions using arrays where size and or shape is not explicitly given as a function argument can never work, since the Fortran subroutine is not provided with that information and couldn’t possibly infer it.

By the way, this also accounts for Fortran derived types (i.e. classes). Where the compiler is in addition permitted to change the layout / order of the members.

Hope that helps. Above code is untested.

4 Likes
6

Thanks for this helpful information. I guess I’ll write some simple wrapper functions then. I can confirm that the following Fortran code:

module foo
    use, intrinsic :: iso_fortran_env, only: int64
    implicit none

contains

    subroutine foo1(n, x, y)
        integer(int64), intent(in) :: n
        integer(int64), intent(in),  dimension(n) :: x
        integer(int64), intent(out), dimension(n) :: y

        y = 2*x
    end subroutine foo1

    subroutine foo2(x, y)
        integer(int64), intent(in),  dimension(:) :: x
        integer(int64), intent(out), dimension(size(x)) :: y

        y = 2*x
    end subroutine foo2

    function bar1(n, x) result(y)
        integer(int64), intent(in) :: n
        integer(int64), intent(in), dimension(n) :: x

        integer(int64), dimension(n) :: y

        y = 2*x
    end function bar1

    function bar2(x) result(y)
        integer(int64), intent(in), dimension(:) :: x

        integer(int64), dimension(size(x)) :: y

        y = 2*x
    end function bar2

end module foo

module foo_wrapper
    use, intrinsic :: iso_c_binding, only: c_long
    use :: foo
    implicit none

contains

    subroutine foo1_wrapper(n, x, y) bind(C, name="foo1")
        integer(c_long), intent(in) :: n
        integer(c_long), intent(in),  dimension(n) :: x
        integer(c_long), intent(out), dimension(n) :: y

        call foo1(n, x, y)
    end subroutine foo1_wrapper

    subroutine foo2_wrapper(n, x, y) bind(C, name="foo2")
        integer(c_long), intent(in) :: n
        integer(c_long), intent(in),  dimension(n) :: x
        integer(c_long), intent(out), dimension(n) :: y

        call foo2(x, y)
    end subroutine foo2_wrapper

    subroutine bar1_wrapper(n, x, y) bind(C, name="bar1")
        integer(c_long), intent(in) :: n
        integer(c_long), intent(in),  dimension(n) :: x
        integer(c_long), intent(out), dimension(n) :: y

        y = bar1(n, x)
    end subroutine bar1_wrapper

    subroutine bar2_wrapper(n, x, y) bind(C, name="bar2")
        integer(c_long), intent(in) :: n
        integer(c_long), intent(in),  dimension(n) :: x
        integer(c_long), intent(out), dimension(n) :: y

        y = bar2(x)
    end subroutine bar2_wrapper

end module foo_wrapper

gives the following output in julia:

julia> t = ones(Int64,2); ccall((:foo1, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ref{Int64}), 2, [3, -1], t); t
2-element Array{Int64,1}:
  6
 -2

julia> t = ones(Int64,2); ccall((:foo2, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ref{Int64}), 2, [3, -1], t); t
2-element Array{Int64,1}:
  6
 -2

julia> t = ones(Int64,2); ccall((:bar1, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ref{Int64}), 2, [3, -1], t); t
2-element Array{Int64,1}:
  6
 -2

julia> t = ones(Int64,2); ccall((:bar2, "./libfoo.so"), Void, (Ref{Int64}, Ref{Int64}, Ref{Int64}), 2, [3, -1], t); t
2-element Array{Int64,1}:
  6
 -2

So it seems to work :slight_smile:

Regards,
Laurent

3 Likes
7

Glad to hear it helped. In retrospect, I forgot to give one important remark regarding types: You might consider using c_int64_t rather than c_long to be sure to always have the 64 bit variant and get automatically a Fortran compiler warning or error if it doesn’t match. See GFortran types for a full overview.

Then, I stumbled today over the this summary of the IBM Fortran compiler, which is incredibly brief and precise in what to do in which situation. And on the bottom of this section it is explained how to change the Fortran code directly.
Still, the ISO C Bindings is the portable way to go…