我正在尝试为twpbvpc(ODE BVP与重新网格化)Fortran-77解算器编写一个包装器 . 求解器需要带签名的输入函数
subroutine fsub(ncomp, x, u, f, rpar, ipar)
哪里
-
ncomp
是整数(向量的长度), -
x
(in)是一个浮点数, -
u
(in)是一个长度为ncomp
的向量, -
f
(out)是结果的位置,长度为ncomp
的向量 -
rpar
和ipar
是浮点数和整数外部参数的数组;朱莉娅的关闭将是更优选的方式,但显然有困难(见the blog post here) . 但有一会儿他们可以被忽视 .
在Julia,写 fsub
我通常会使用签名
function fsub_julia(x :: Float64, y :: Vector{Float64}, dy :: Vector{Float64})
dy[1] = ...
dy[2] = ...
...
end
ncomp
似乎没有必要,因为可以通过 length
或 size
来获取长度(但是,Julia可以检测从Fortran传递的数组的大小吗?对于测试代码,我明确知道 ncomp
,所以现在它不是问题) .
所以,为了遵守twpbvpc格式,我写了一个包装器:
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
现在,要将此函数传递给Fortran例程,需要使用 cfunction
转换它来声明类型 . 问题是如何?
如果我把它作为:
cf_fsub = cfunction(fsub_par, Void, (Ref{Int64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int64}))
从Fortran调用时,我收到错误:
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})
所以,不知何故签名与数组参数不匹配......
如果我用 Ref{Array{Float64,1}}
替换 Ref{Float64}
数组参数(虽然看起来有点奇怪......):
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}}))
在Fortran代码中调用 fsub_par
( cf_fsub
)时,我得到了一个分段错误(这个位于大约,因为错误没有给出确切的位置) .
用数字 Ptr{Float64}
替换 Ref{Float54}
也不起作用 .
我在Fortran代码中找到的一件有趣的事情是如何调用 fsub
:
call fsub (ncomp, xx(1), u(1,1), fval(1,1),rpar,ipar)
其中 u
和 fval
被声明为:
dimension xx(nmsh), u(nudim,nmsh), fval(ncomp,nmsh)
所以,我猜,它使用了Fortran通过引用传递所有参数的事实,并且对 u(1,1)
的引用应该是指向矩阵第一列的指针(据我记得Fortran,就像Julia一样,存储列中的矩阵 - 第一个订单) .
它的出路是什么?我是否需要更改 fsub_julia
的签名以接受指针并手动将它们转换为数组(这是ODEInterface.jl在较低级别的工作方式)?
更新
按照它在ODEInterface.jl中完成的方式,并结合在C中的 void*
-thunk参数中传递Julia的函数的想法,我提出了这个:
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}))
当我打电话给我做的解算器时(它很长):
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
Fortran的 twpbvpc
看起来像这样(显然是开头):
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代码使用 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`)
所以,我传递 rpar
为 Any
(应该相当于 Ptr{Void}
):Fortran期望一个浮点数的数组,但它应该没关系 .
现在,当我尝试运行一个简单的程序(在 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
随着代码越来越长,这里是github上完整代码的链接:https://github.com/mobius-eng/TWPBVP.jl
1 回答
是的,ODEInterface.jl模型看起来非常好 .
您需要找出的第一件事是您的fortran
INTEGER
类型的大小(Int32
或Int64
) . 对于下面的代码,我将从ODEInterface.jl借用并使用FInt
(它可以是类型参数,也可以是typealias
)产生的后备应该类似于: