Different output when call Fortran from Julia

General Usage
, , ,
1

Hello!

I am trying to use a package in Fortran, and I tested it in Julia - it works perfectly well. Now I am trying to write a function in Julia that I can use it as a normal Julia function. Here is the function code:

function mvndst_j(n::Int64, lower::Vector{Float64}, upper::Vector{Float64}, infin::Vector{Int64},
                correl::Vector{Float64}, maxpts::Int64, abseps::Float64, releps::Float64,
                ift::Int64)
    n_ref = Ref{Cint}(n)
    lower_ptr = pointer(lower)
    upper_ptr = pointer(upper)
    infin_ptr = pointer(infin)
    correl_ptr = pointer(correl)
    maxpts_ref = Ref{Cint}(maxpts)
    abseps_ref = Ref{Cdouble}(abseps)
    releps_ref = Ref{Cdouble}(releps)
    err = zeros(Cdouble, 1)
    val = zeros(Cdouble, 1)
    ift_ref = Ref{Cint}(ift)
    err,val=ccall(mvndst_func, Tuple{Cdouble, Cdouble},
               (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
                Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
               n_ref, lower_ptr, upper_ptr, infin_ptr, correl_ptr, maxpts_ref,
               abseps_ref, releps_ref, err, val, ift_ref)
    return val
end

I am using the value below to test it:

n = 2
lower = [0.0, 0.0]
upper = [1.0, 2.0]
infin = [0, 0]
correl = [0.5]
maxpts = 10000
abseps = 1e-5
releps = 0.0
ift = 0

When I run with the inside of the code,

    lower_ptr = pointer(lower)
    upper_ptr = pointer(upper)
    infin_ptr = pointer(infin)
    correl_ptr = pointer(correl)
    maxpts_ref = Ref{Cint}(maxpts)
    abseps_ref = Ref{Cdouble}(abseps)
    releps_ref = Ref{Cdouble}(releps)
    err = zeros(Cdouble, 1)
    val = zeros(Cdouble, 1)
    ift_ref = Ref{Cint}(ift)
    err,val=ccall(mvndst_func, Tuple{Cdouble, Cdouble},
               (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
                Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
               n_ref, lower_ptr, upper_ptr, infin_ptr, correl_ptr, maxpts_ref,
               abseps_ref, releps_ref, err, val, ift_ref)

It gave me correct value of 0.8318608311308805 but when I tried to use

mvndst_j(n, lower, upper, infin, correl, maxpts, abseps, releps, ift)

It gave me wrong answer of 0.37260179342337485

I have no clue why it happened. (I tried printing values inside the function and it matched with the value I entered.)

Thanks a lot!

2

Update from myself - seems to be a pointer problem in Fortran.

2 Likes
3

Is the issue solved?

5

I am not sure if that’s the reason, but pointer seems to assign new memory after I add

err_ptr = ccall(:malloc, Ptr{Cdouble}, (Csize_t,), sizeof(Cdouble))
  val_ptr = ccall(:malloc, Ptr{Cdouble}, (Csize_t,), sizeof(Cdouble))
6

I’m no expert but I do use ccall.
Have you tried using pointer() on err and val in your ccall ? e.g.

err,val=ccall(mvndst_func, Tuple{Cdouble, Cdouble},
               (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
                Ptr{Cdouble}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble},
                Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
               n_ref, lower_ptr, upper_ptr, infin_ptr, correl_ptr, maxpts_ref,
               abseps_ref, releps_ref, pointer(err), pointer(val), ift_ref)

to make sure you’re not passing a copy in case there is type casting.

7

This is dangerous without GC.@preserve, because it’s not rooted (garbage collection may destroy the underlying object before you are done with the pointer). Better to let ccall do the pointer conversions for you.

That is, I would do something like:

infin = Cint[0, 0] # need a Cint array if you want to pass it for Ptr{Cint}
err, val = Ref{Cdouble}(), Ref{Cdouble}() # need explicit Refs if these are "output" arguments
@ccall mvndst_proc(n::Ref{Cint}, lower::Ptr{Cdouble}, upper::Ptr{Cdouble}, infin::Ptr{Cint},
                   correl::Ref{Cdouble}, maxpts::Ref{Cint}, abseps::Ref{Cdouble}, releps::Ref{Cdouble},
                   err::Ref{Cdouble}, val::Ref{Cdouble}, ift::Ref{Cint})::Cvoid
return err[], val[]

Note that

  • Julia will (safely) convert array arguments to pointer for you and will convert scalars to Ref for you — no need to explicitly construct Ptr and Ref objects. Simply pass lower, maxpts, etc. directly to the ccall. (Declare array arguments as Ptr and Fortran scalar arguments as Ref.)
  • The only time you need to construct an explicit Ref object is for “output” arguments, whose values get changed in the function (and you want to access the new arguments).
  • Returning a Tuple{Cdouble, Cdouble} generally will not work — you cannot easily return a Julia tuple from a C or Fortran routine. Instead, for Fortran I would typically recommend using a subroutine with no return value (declared as Cvoid in Julia), and returning the result via output arguments as in the example above.
  • @ccall is a lot less error prone than the lower-level ccall syntax, because the argument types go next to the argument values.
2 Likes