Fortran arrays in Julia callbacks

Question

I am trying to write a wrapper for twpbvpc (ODE BVP with remeshing) Fortran-77 solver. The solver needs the input function with signature

subroutine fsub(ncomp, x, u, f, rpar, ipar)

where

• ncomp is an integer (length of the vector),
• x (in) is a float,
• u (in) is a vector of length ncomp,
• f (out) is the location for the result, vector of length ncomp
• rpar and ipar are arrays of float and integer external parameters; Julia's closure would be more preferred way, but apparently there are difficulties with it (see the blog post here). But for a moment they can be ignored.

In Julia, to write fsub I would normally use the signature

function fsub_julia(x :: Float64, y :: Vector{Float64}, dy :: Vector{Float64})
        dy[1] = ...
        dy[2] = ...
        ...
end

ncomp doesn't seem to be necessary as can get the length via length or size (however, can Julia detect the size of the passed arrays from Fortran? for a test code, I know ncomp explicitly, so for now it's not a problem).

So, to comply with twpbvpc format I wrote a wrapper:

function fsub_par(n :: Int64, x :: Float64, y :: Vector{Float64}, dy :: Vector{Float64}, rpar :: Vector{Float64}, ipar :: Vector{Float64})
        fsub_julia(x, y, dy)
end

Now, to pass this function to Fortran routine, it needs to be converted using cfunction to declare the types. The question is how?

If I put it as:

cf_fsub = cfunction(fsub_par, Void, (Ref{Int64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int64}))

When called from Fortran, I get the error:

ERROR: LoadError: MethodError: no method matching (::TWPBVP.#fsub_par#1{TWPBVP_Test.#f})(::Int64, ::Float64, ::Float64, ::Float64, ::Float64, ::Int64)
Closest candidates are:
  fsub_par(::Int64, ::Float64, !Matched::Array{Float64,1}, !Matched::Array{Float64,1}, !Matched::Array{Float64,1}, !Matched::Array{Float64,1})

So, somehow the signature doesn't match for the array arguments...

If I replace Ref{Float64} for an array arguments with Ref{Array{Float64,1}} (it does look a bit odd though...):

cf_fsub = cfunction(fsub_par, Void, (Ref{Int64}, Ref{Float64}, Ref{Array{Float64,1}}, Ref{Array{Float64,1}}, Ref{Array{Float64,1}}, Ref{Array{Int64,1}}))

I get a segmentation fault at the point when fsub_par (cf_fsub) is called in Fortran code (this located approximately, as the error doesn't give exact location).

Replacing Ref{Float54} with Ptr{Float64} for arrays also doesn't do anything.

One interesting thing I found in the Fortran code is how fsub is called:

call fsub (ncomp, xx(1), u(1,1), fval(1,1),rpar,ipar)

where u and fval are declared as:

dimension xx(nmsh), u(nudim,nmsh), fval(ncomp,nmsh)

So, I guess, it uses the fact that Fortran passes all the arguments by reference and reference to u(1,1) is supposed to be the pointer to the first column of the matrix (as far as I remember Fortran, as Julia, stores matrices in column-first order).

What is the way out of it? Do I need to change the signature of fsub_julia to accept the pointers and convert them to arrays manually (this is how ODEInterface.jl works at the lower level)?

Update

Following how it was done in ODEInterface.jl and combining with the idea of passing Julia's functions in void*-thunk parameters in C I came up with this:

immutable TWPBVPCProblem
    fsub :: Function            # RHS function
    dfsub :: Function           # Jacobian of RHS
    gsub :: Function            # BC function
    dgsub :: Function           # gradients of BC function
end

function unsafe_fsub(rn :: Ref{Int64}, rx :: Ref{Float64}, py :: Ptr{Float64}, pdy :: Ptr{Float64}, rpar :: Ptr{Float64}, ipar :: Ptr{Int64}) :: Void
    x = rx[]
    n = rn[]
    y = unsafe_wrap(Array, py, n)
    dy = unsafe_wrap(Array, pdy, n)
    problem = unsafe_pointer_to_objref(rpar) :: TWPBVPCProblem
    problem.fsub(x, y, dy)
    return nothing
end

const fsub_ptr = cfunction(unsafe_fsub, Void, (Ref{Int64}, Ref{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}))

And when I call the solver I do (it's rather long):

function twpbvpc(nlbc :: Int64,
        aleft :: Float64, aright :: Float64,
        fixpnt :: Nullable{Vector{Float64}},
        ltol :: Vector{Int64}, tol :: Vector{Float64},
        linear :: Bool, givmsh :: Bool, giveu :: Bool, nmsh :: Ref{Int64},
        xx :: Vector{Float64}, u :: Array{Float64, 2}, nmax :: Ref{Int64},
        wrk :: Vector{Float64}, iwrk :: Vector{Int64},
        fsub :: Function, dfsub :: Function,
        gsub :: Function, dgsub :: Function,
        ckappa1 :: Ref{Float64}, gamma1 :: Ref{Float64},
        ckappa :: Ref{Float64},
        # rpar :: Vector{Float64},
        # ipar :: Vector{Int64},
        iflbvp :: Ref{Int64})

    # Keep problem functions in rpar
    # HACK!
    rpar = TWPBVPCProblem(fsub, dfsub, gsub, dgsub)

    # Fake external parameters
    # Can't have it 0-length as it would be Any[0] and not Float64[0]
    # local rpar :: Vector{Float64} = [0.0]
    local ipar :: Vector{Int64} = [0]

    # No need to pass these parameters
    # u is a matrix for the solution only!
    ncomp, nucol = size(u)
    # Get the maximum of xx
    nxxdim = length(xx)
    # max for mesh points must be the same as the number of column points of u
    assert(nucol == nxxdim)

    # Sizes of work arrays
    lwrkfl = length(wrk)
    lwrkin = length(iwrk)

    # Number of fixed mesh points
    if isnull(fixpnt)
        nfxpnt = 0
        fixpnt_v = [0.0]
    else
        fixpnt_v = get(fixpnt)
        nfxpnt = length(fixpnt_v)
    end

    # Size of tolerance vector ≤ ncomp
    ntol = length(ltol)

    ccall((:twpbvpc_, libtwpbvpc), Void,
        (Ref{Int64}, Ref{Int64},                    # ncomp, nlbc,
        Ref{Float64}, Ref{Float64},                 # aleft, aright
        Ref{Int64}, Ptr{Float64},                   # nfxpnt, fixpnt
        Ref{Int64}, Ptr{Int64}, Ptr{Float64},       # ntol, ltol, tol
        Ref{Int64}, Ref{Int64}, Ref{Int64},         # linear, givmsh, giveu
        Ref{Int64}, Ref{Int64},                     # nmsh, nxxdim
        Ptr{Float64}, Ref{Int64},                   # xx, nudim
        Ptr{Float64}, Ref{Int64},                   # u, nmax
        Ref{Int64}, Ptr{Float64},                   # lwrkfl, wrk
        Ref{Int64}, Ptr{Int64},                     # lwrkin, iwrk
        Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, # fsub, dfsub, gsub, dgsub
        Ref{Float64}, Ref{Float64},                 # ckappa1, gamma1
        Ref{Float64}, Any, Ptr{Int64},              # ckappa, rpar, ipar
        Ref{Int64}),                                # iflbvp
        ncomp, nlbc, aleft, aright,
        nfxpnt, fixpnt_v, ntol, ltol, tol,
        linear, givmsh, giveu, nmsh,
        nxxdim, xx, nucol, u, nmax,                 # nudim = nucol
        lwrkfl, wrk, lwrkin, iwrk,
        fsub_ptr, dfsub_ptr, gsub_ptr, dgsub_ptr,
        ckappa1,gamma1,ckappa,pointer_from_objref(rpar),ipar,iflbvp)
end

The Fortran's twpbvpc looks like this (the beginning obviously):

  subroutine twpbvpc(ncomp, nlbc, aleft, aright,
 *       nfxpnt, fixpnt, ntol, ltol, tol,
 *       linear, givmsh, giveu, nmsh,
 *       nxxdim, xx, nudim, u, nmax,
 *       lwrkfl, wrk, lwrkin, iwrk,
 *       fsub, dfsub, gsub, dgsub,
 *       ckappa1,gamma1,ckappa,rpar,ipar,iflbvp)

  implicit double precision (a-h,o-z)
  dimension rpar(*),ipar(*)
  dimension fixpnt(*), ltol(*), tol(*)
  dimension xx(*), u(nudim,*)
  dimension wrk(lwrkfl), iwrk(lwrkin)

  logical linear, givmsh, giveu
  external fsub, dfsub, gsub, dgsub

  logical pdebug, use_c, comp_c
  common/algprs/ nminit, pdebug, iprint, idum, uval0, use_c, comp_c
  ...

Fortran code is compiled with build.jl:

cd(joinpath(Pkg.dir("TWPBVP"), "deps"))
pic = @windows ? "" : "-fPIC"
run(`gfortran -m$WORD_SIZE -fdefault-real-8 -fdefault-integer-8 -ffixed-form $pic -shared -O3 -o libtwpbvpc.so twpbvpc.f`)

So, I'm passing rpar as Any (should be equivalent to Ptr{Void}): Fortran expects an array of floating-point numbers though, but it shouldn't matter.

Now I'm getting segmentation fault when I try to run a simple program (on Pkg.test("TWPBVP")):

signal (11): Segmentation fault
while loading /home/alexey/.julia/v0.5/TWPBVP/test/runtests.jl, in expression starting on line 58
unknown function (ip: 0xffffffffffffffff)
Allocations: 1400208 (Pool: 1399373; Big: 835); GC: 0

As the code is getting very long, here is the link to the full code on github: https://github.com/mobius-eng/TWPBVP.jl

Answer 1

Do I need to change the signature of fsub_julia to accept the pointers and convert them to arrays manually (this is how ODEInterface.jl works at the lower level)?

Yes, the ODEInterface.jl model looks like a pretty good one to follow.

The first thing you need to find out is the size of your fortran INTEGER type (either Int32 or Int64). For the code below I'll borrow from ODEInterface.jl and use FInt (it could either be a type parameter, or a typealias)

The resulting fallback should look something like:

# SUBROUTINE FSUB(NCOMP,X,Z,F,RPAR,IPAR)
# IMPLICIT NONE
# INTEGER NCOMP, IPAR
# DOUBLE PRECISION F, Z, RPAR, X
# DIMENSION Z(*),F(*)
# DIMENSION RPAR(*), IPAR(*)
function unsafe_fsub(ncomp::Ref{FInt}, x::Ref{Float64}, z::Ptr{Float64}, 
        f::Ptr{Float64}, rpar::Ptr{Float64}, ipar::Ptr{FInt})::Void
    xx = x[]
    zz = unsafe_wrap(Array, z, ncomp[])
    ff = unsafe_wrap(Array, f, ncomp[])
    fsub!(xx, zz, ff) # function which updates array ff
    return nothing
end

const fsub_ptr = cfunction(unsafe_fsub, Void,
    (Ref{FInt},Ref{Float64},Ptr{Float64},Ptr{Float64},Ptr{Float64},Ptr{FInt}))