Error calling Fortran from Julia

I’ve been trying to make a very simple example of Julia calling fortran to work but it seems I’m missing something, right now what I’ve done is:

Fortran code:

module ipmf
  use, intrinsic :: iso_c_binding
  implicit none
  private
  public :: wavg_f
subroutine wavg_f(x, w, n, wavg) bind(C, name="wavg_f_")

  integer(kind = c_int), intent(in), value        :: n    !Length of x and w
  real(kind = c_double), intent(in), dimension(n) :: x    !Vector of values
  real(kind = c_double), intent(in), dimension(n) :: w    !Vector of weights
  real(kind = c_double), intent(out)              :: wavg !output
  real(kind = c_double)                           :: wsum !Sum of weights
  integer                                         :: i    !Internal count

  wavg = 0.0
  wsum = 0.0

  do i = 1, n
      wsum = wsum + w(i)
  end do

  do i = 1, n
      wavg = wavg + x(i)*w(i)/wsum
  end do

end subroutine wavg_f
end module ipmf


Compiled the code using:

> gfortran -shared -fPIC ipmf.f95 -o ipmf.so

and finally, my Julia code looks like this:

using Libdl
lib = Libdl.dlopen("libs/ipmf.so")
sym = Libdl.dlsym(lib, :wavg_f_)
x = Float64[1,2,3]
w = Float64[1,1,1]
n = 3
wavg = 0.0
ccall(sym, Cvoid, (Ref{Float64}, Ref{Float64}, Ref{Int64}, Ref{Float64}), x, w, Ref(n), Ref(wavg))
    

but instead of getting wavg = 2.0, I’m getting this:

julia> wavg
0.0

I’m sure I’m doing something wrong, I do not have that much experience using Julia… I’d appreciate any help. Thanks!

Your Fortran code declares this parameter as integer(kind = c_int), which would correspond to Ref{Cint} (Cint == Int32, a 32-bit integer) not Int64.

Also, the array arguments x and w should probably be declared as Ptr{Float64} (or equivalently Ptr{Cdouble}) rather than Ref (which is used to wrap scalars in pointers).

4 Likes

Thank you for your answer, I tried some variatios of what you suggest, but I keep getting the same value or (sometimes) I get this error:

julia> ccall(sym, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ref{Cint}, Ref{Cdouble}), x, w, Ref(n), Ref(wavg))
ERROR: ReadOnlyMemoryError()
Stacktrace:
 [1] top-level scope at ./REPL[28]:

My code right now looks like this:

using Libdl
lib = Libdl.dlopen("libs/ipmf.so")
sym = Libdl.dlsym(lib, :wavg_f_)
x = Cdouble[1,2,3]
w = Cdouble[1,1,1]
n = convert(Cint, 3)
wavg = 0.0::Cdouble
ccall(sym, Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ref{Cint}, Ref{Cdouble}), x, w, Ref(n), Ref(wavg))

I also created an smaller example avoiding vectors… but somehow I’m making the same mistake:

Fortran code:

subroutine sum(a, b, c) bind(C, name = "sum_")
  real(kind = c_double), intent(in)    :: a, b 
  real(kind = c_double), intent(out)   :: c
  c = 0.0
  c = a+b
end subroutine sum

and my Julia Code:

lib = Libdl.dlopen("libs/ipmf.so")
sym_sum = Libdl.dlsym(lib, :sum_)
a=1.0
b=2.0
c=0.0
ccall(sym_sum, Cvoid, (Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), Ref{Cdouble}(a), Ref{Cdouble}(b), Ref{Cdouble}(c))

And I’m getting

julia> c
0.0

Instead of 3.0

In this example, c is not changed because it is never passed to the Fortran function (you are passing a Ref that is initialized to the same value as c, but isn’t related to c in any other way). Define c = Ref{Cdouble}(0.0) from the start and use that in the ccall.

1 Like

Thank you @sostock, with your suggestion I finally got some numbers out of fortran. Somehow I kept having problems reading integers. Everytime I wanted to print them I got trash. I had to change the way I was sending them in the Julia code and now is working

using Libdl
lib = Libdl.dlopen("libs/ipmf.so")
sym = Libdl.dlsym(lib, :wavg_f_)
x = Cdouble[1.0,2.0,3.0]
w = Cdouble[1.0,1.0,1.0]
n = 3
wavg = Ref{Cdouble}(50.0)
ccall(sym, Cvoid, (Ref{Cdouble}, Ref{Cdouble}, Cint, Ref{Cdouble}), x, w, n, wavg)

For reasons I don’t yet understand Ref{Cint} / Ref© wasn’t passed through properly and because I’m using that variable to define my arrays it make everything crash.

Thank you very much for your help

Since you are using bind(C, ...) then scalar in values are passed through using C calling conventions, not as pointers as would happen for Fortran functions… So n should be passed as Cint, not as Ref{Cint}.

4 Likes