Correct way to ccall a Fortran subroutine requiring pointers

General Usage
1

Given the following Fortran code:

MODULE PTR
   IMPLICIT NONE
   REAL, SAVE, TARGET, ALLOCATABLE :: A(:, :), B(:, :)

CONTAINS

   SUBROUTINE INIT(AP, BP)
      REAL, POINTER :: AP(:, :), BP(:, :)

      ALLOCATE(A(6, 6), B(6, 6))
      AP => A(1:4, 2:5)
      BP => B(2:3, 3:6)

      AP = 3.2
      PRINT *, AP

      BP = 4.5
      PRINT *, BP
   END SUBROUTINE
END MODULE

Currently, I can only ccall the subroutine with the following hacky method:

using Libdl

lib = dlopen("minimal.so")

none = () -> nothing
ap = Ref(@cfunction(none, Cvoid, ()))
bp = Ref(@cfunction(none, Cvoid, ()))

ccall(dlsym(lib, :__ptr_MOD_init), Cvoid, (Ref{Ptr{Cvoid}}, Ref{Ptr{Cvoid}}), ap, bp)

a = unsafe_wrap(Array{Float32,2}, convert(Ptr{Float32}, ap[]), (6, 6))
@show a[1:4, 1:4]

b = unsafe_wrap(Array{Float32,2}, convert(Ptr{Float32}, bp[]), (6, 6))
@show b[1:2, 1:4]

Output is:

   3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005    
   4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000    
a[1:4, 1:4] = Float32[3.2 3.2 3.2 3.2; 3.2 3.2 3.2 3.2; 3.2 3.2 3.2 3.2; 3.2 3.2 3.2 3.2]
b[1:2, 1:4] = Float32[4.5 4.5 4.5 4.5; 4.5 4.5 4.5 4.5]

I used @cfunction to get a function pointer here and pass them to the Fortran subroutine.

All other attempts failed. For example, if I use ap = Ref{Ptr{Float32}}(C_NULL) and change the ccall arguments accordingly,

ap = Ref{Ptr{Float32}}(C_NULL)
bp = Ref{Ptr{Float32}}(C_NULL)

ccall(dlsym(lib, :__ptr_MOD_init), Cvoid, (Ref{Ptr{Float32}}, Ref{Ptr{Float32}}), ap, bp)

The output will be:

   4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000    
a[1:4, 1:4] = Float32[0.0 1.0f-44 7.0f-45 4.3228914f-38; 0.0 0.0 0.0 0.0; 9.0f-44 0.0 1.1f-44 0.0; 0.0 0.0 0.0 1.0f-45]

signal (11): Segmentation fault
...error messages omitted...

Using ap = zeros(Float32, 6, 6), i.e.

ap = zeros(Float32, 6, 6)
bp = zeros(Float32, 6, 6)

ccall(dlsym(lib, :__ptr_MOD_init), Cvoid, (Ptr{Float32}, Ptr{Float32}), ap, bp)

@show ap
@show bp

will not cause a segmentation fault, however, it is obvious that we cannot get the correct data:

   3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005       3.20000005    
   4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000       4.50000000    
ap = Float32[4.7818537f-38 0.0 1.0f-45 1.0f-45 0.0 0.0; 0.0 1.079f-42 0.0 0.0 0.0 0.0; NaN 6.0f-45 6.0f-45 6.0f-45 0.0 0.0; NaN 0.0 0.0 0.0 0.0 0.0; 6.0f-45 1.0f-45 8.0f-45 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]
bp = Float32[1.7421598f-38 0.0 1.0f-45 1.0f-45 0.0 0.0; 0.0 1.079f-42 0.0 0.0 0.0 0.0; NaN 6.0f-45 3.0f-45 6.0f-45 0.0 0.0; NaN 0.0 0.0 0.0 0.0 0.0; 6.0f-45 1.0f-45 8.0f-45 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0]

Can anyone explain the internal mechanism? Is there a better way to handle Fortran pointers other than using function pointers?

Another question

Another question is about the indices. Note that in the first example, b[1:2, 1:4] in Julia corresponds to B(2:3, 3:6) in Fortran, is it always true that indices will start from 1 at the Julia side?

2

I think the answer to your general question is this one: Calling Fortran variables from Julia - #3 by anowacki

Note that this is not specific of julia, but of the interface. Since you pass a reference, one side does not know about how the indexing was at the other side:

program main
    integer :: i, a(3), b
    do i = 1,3
        a(i) = i
    end do
    print *, a
    call test(a(2:3),b)
    print *, b
end

subroutine test(a,b)
    integer :: a(2), b
    b = a(1)
end

Will print:

           1           2           3
           2

which means that b assumed the value of what was a(2) in the main program.

You can, in Julia, use OffsetArrays to use arbitrary indexing (but then you need one additional layer of care in passing the array to fortran, since an OffsetArray is an array contained in a struct that carries the information about the offsets.

1 Like
3

Thanks! What about the first question? How should we pass Fortran pointers correctly (if the Fortran code cannot be modified)?

4

I would probably write, in Fortran, an interfacing module, following those directives, and call that module from Julia instead of the original one.

1 Like