Different output when call Fortran from Julia

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!

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

2 Likes

Is the issue solved?

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))