Weird behavior calling FORTRAN function

I have a Julia function that calls a FORTRAN subroutine that I have written. Everything was working fine until recently I tried to change the return type of the Julia function - an array of the FORTRAN returns instead of a tuple. i.e.

return x, y, z -> return hcat(x, y, z).

Having changed this, and nothing else, the code now crashes within the FORTRAN routine from an error that would mean the input to the subroutine has somehow changed to bad values. Can someone help me make sense of this?

Edit: The function in question:

function ComputeDerivatives_Pressure( D::Array{Float64, 1}, E::Array{Float64, 1}, Ne::Array{Float64,1};
    Units_Option::Bool=true )

    # Initialize arrays to hold derivatives
    nx     :: Int64             = length( D )
    dPdD   :: Array{Float64, 1} = zeros( nx );
    dPdT   :: Array{Float64, 1} = zeros( nx );
    dPdY   :: Array{Float64, 1} = zeros( nx );
    dPdE   :: Array{Float64, 1} = zeros( nx );
    dPdDe  :: Array{Float64, 1} = zeros( nx ); 
    dPdTau :: Array{Float64, 1} = zeros( nx ); 

    # =============================================================
    # This calls the FORTRAN function :computederivatives_pressure_
    # FORTRAN compilation mangles the name. Cvoid is the return 
    # type, the next parameters are input types, followed 
    # by the arguements. 
    # =============================================================
    ccall( (:computederivatives_pressure_, "./EoS_jl.so"), Cvoid, 
    ( Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, 
    Ref{Float64}, Ref{Bool} ), 
    D, E, Ne, nx, dPdD, dPdT, dPdY, dPdE, dPdDe, dPdTau, Units_Option )    

    # return dPdD, dPdT, dPdY, dPdE, dPdDe, dPdTau
    return hcat( copy(dPdD), copy(dPdT), copy(dPdY), copy(dPdE), copy(dPdDe), copy(dPdTau) )

end   

Update: Even without changing the return type, adding a line such as y = hcat( zeros(nx), zeros(nx), zeros(nx), zeros(nx) ) after the FORTRAN call causes a crash in the FORTRAN.

I think all of your arrays should be passed to Fortran with a Vector{Float64} signature, not Ref{Float64}. As far as I can tell, the only Ref items in the list should be the nx and Units_option variables, which should also have explicit Ref wrappers (assuming it’s still using traditional Fortran calling convention and hasn’t been setup with ISO C bindings instead):

function ComputeDerivatives_Pressure(D::Array{Float64, 1},
                                     E::Array{Float64, 1},
                                     Ne::Array{Float64,1};
                                     Units_Option::Bool=true)

    # Initialize arrays to hold derivatives
    nx     :: Int64             = length( D )
    dPdD   :: Array{Float64, 1} = zeros( nx );
    dPdT   :: Array{Float64, 1} = zeros( nx );
    dPdY   :: Array{Float64, 1} = zeros( nx );
    dPdE   :: Array{Float64, 1} = zeros( nx );
    dPdDe  :: Array{Float64, 1} = zeros( nx ); 
    dPdTau :: Array{Float64, 1} = zeros( nx ); 

    ccall((:computederivatives_pressure_, "./EoS_jl.so"), Cvoid, 
        (Vector{Float64}, Vector{Float64}, Vector{Float64}, Ref{Int64},
         Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64},
         Vector{Float64}, Vector{Float64}, Ref{Bool}), 
        D, E, Ne, Ref(nx), dPdD, dPdT, dPdY, dPdE, dPdDe, dPdTau, Ref(Units_Option))

    # return dPdD, dPdT, dPdY, dPdE, dPdDe, dPdTau
    return hcat(dPdD, dPdT, dPdY, dPdE, dPdDe, dPdTau)
end

No, they should be Ref (or Ptr, which is equivalent), see Calling fortran function.

2 Likes

Double check that you are using the correct types for nx and Units_Option: if you’re using gfortran and the default KINDs for INTEGER and LOGICAL, then I think these should both be Int32: https://gcc.gnu.org/onlinedocs/gfortran/KIND-Type-Parameters.html

1 Like

This is also explained in Julia’s documentation: https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/index.html#man-bits-types-1, see the “System Independent Types” table

1 Like

Thanks for correcting me — I misremembered what I’d read [apparently too] long ago in the manual: https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/#When-to-use-T,-Ptr{T}-and-Ref{T}-1