Error calling Fortran from Julia

General Usage
1

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!

2

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
3

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))
4

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

5

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
6

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(c) 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

7

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