Fortran arrays in Julia callbacks -
i trying write wrapper twpbvpc (ode bvp remeshing) fortran-77 solver. solver needs input function signature
subroutine fsub(ncomp, x, u, f, rpar, ipar) where
ncompinteger (length of vector),x(in) float,u(in) vector of lengthncomp,f(out) location result, vector of lengthncomprpar,ipararrays of float , integer external parameters; julia's closure more preferred way, apparently there difficulties (see the blog post here). moment can ignored.
in julia, write fsub use signature
function fsub_julia(x :: float64, y :: vector{float64}, dy :: vector{float64}) dy[1] = ... dy[2] = ... ... end ncomp doesn't seem necessary can length via length or size (however, can julia detect size of passed arrays fortran? test code, know ncomp explicitly, it's not problem).
so, comply twpbvpc format wrote 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, pass function fortran routine, needs converted using cfunction declare types. question how?
if put as:
cf_fsub = cfunction(fsub_par, void, (ref{int64}, ref{float64}, ref{float64}, ref{float64}, ref{float64}, ref{int64})) when called fortran, 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 signature doesn't match array arguments...
if replace ref{float64} array arguments ref{array{float64,1}} (it 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 segmentation fault @ point when fsub_par (cf_fsub) called in fortran code (this located approximately, error doesn't give exact location).
replacing ref{float54} ptr{float64} arrays doesn't anything.
one interesting thing found in fortran code how fsub called:
call fsub (ncomp, xx(1), u(1,1), fval(1,1),rpar,ipar) where u , fval declared as:
dimension xx(nmsh), u(nudim,nmsh), fval(ncomp,nmsh) so, guess, uses fact fortran passes arguments reference , reference u(1,1) supposed pointer first column of matrix (as far remember fortran, julia, stores matrices in column-first order).
what way out of it? need change signature of fsub_julia accept pointers , convert them arrays manually (this how odeinterface.jl works @ lower level)?
update
following how done in odeinterface.jl , combining idea of passing julia's functions in void*-thunk parameters in c came 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 call solver (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 0-length any[0] , not float64[0] # local rpar :: vector{float64} = [0.0] local ipar :: vector{int64} = [0] # no need pass these parameters # u matrix solution only! ncomp, nucol = size(u) # maximum of xx nxxdim = length(xx) # max mesh points must same 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 (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 compiled 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 any (should equivalent ptr{void}): fortran expects array of floating-point numbers though, shouldn't matter.
now i'm getting segmentation fault when try run 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 code getting long, here link full code on github: https://github.com/mobius-eng/twpbvp.jl
do need change signature of fsub_julia accept pointers , convert them arrays manually (this how odeinterface.jl works @ lower level)?
yes, odeinterface.jl model looks pretty 1 follow.
the first thing need find out size of fortran integer type (either int32 or int64). code below i'll borrow odeinterface.jl , use fint (it either type parameter, or typealias)
the resulting fallback should 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 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}))
Comments
Post a Comment