Calling fortran function

New to Julia
,
1

Trying to call a fortran function. The fortran library is qrupdate. The fortran files were compiled with flags

-fimplicit-none -O3 -funroll-loops -fPIC

using gfortran.

My julia command looks like

ccall(("sqrinr_","qrupdate-1.1.2/libqrupdate"), Int32, (Int64, Int64, Array{Float64, 2}, Int64, Array{Float64, 2}, Int64, Int64, Array{Float64, 1}, Array{Float64, 1}), n, q, Q, n+1, R, n+1, n, xp, w)

and I get a segfault

signal (11): Segmentation fault
while loading no file, in expression starting on line 0
sqrinr_ at qrupdate-1.1.2/libqrupdate.so (unknown line)
anonymous at ./<missing> (unknown line)
unknown function (ip: 0x7f084dec901f)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_toplevel_eval_flex at /home/centos/buildbot/slave/package_tarball64/build/src/toplevel.c:569 [inlined]
jl_toplevel_eval at /home/centos/buildbot/slave/package_tarball64/build/src/toplevel.c:580
jl_toplevel_eval_in_warn at /home/centos/buildbot/slave/package_tarball64/build/src/builtins.c:590
eval at ./boot.jl:234
unknown function (ip: 0x7f0a58ebf4df)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
eval_user_input at ./REPL.jl:64
unknown function (ip: 0x7f084de9a0a6)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
macro expansion at ./REPL.jl:95 [inlined]
#3 at ./event.jl:68
unknown function (ip: 0x7f084de971cf)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
jl_apply at /home/centos/buildbot/slave/package_tarball64/build/src/julia.h:1392 [inlined]
start_task at /home/centos/buildbot/slave/package_tarball64/build/src/task.c:253
unknown function (ip: 0xffffffffffffffff)
Allocations: 6546676 (Pool: 6544886; Big: 1790); GC: 10
Segmentation fault

I have no idea how to approach the problem, can someone point me in right direction?

2 Likes
2

I don’t have much experience here, but I think the issue is that fortran expects all arguments to be passed by reference (see here: http://docs.julialang.org/en/stable/manual/calling-c-and-fortran-code/#calling-c-and-fortran-code ). There’s a note in the ccall docs that says as much:

When calling a Fortran function, all inputs must be passed by reference, so all type correspondences above should contain an additional Ptr{…} or Ref{…} wrapper around their type specification.

3

Thanks I will try that out.

4

I changed the ccall to

julia> ccall(("sqrinr_","qrupdate-1.1.2/libqrupdate"), Void, (Ref{Int64}, Ref{Int64}, Ref{Array{Float64, 2}}, Ref{Int64}, Ref{Array{Float64, 2}}, Ref{Int64}, Ref{Int64}
, Ref{Array{Float64, 1}}, Ref{Array{Float64, 1}}), n, q, Q, n+1, R, n+1, n, xp, w)

and it seemed to run without error, but when I look at Q I get a segfault. The fortran function is supposed to modify Q (and R). Any more hints and tips?

julia> Q
4×1 Array{Float64,2}:

signal (11): Segmentation fault
while loading no file, in expression starting on line 0
jl_array_isassigned at /home/centos/buildbot/slave/package_tarball64/build/src/array.c:487
isassigned at ./array.jl:36
alignment at ./show.jl:1277
unknown function (ip: 0x7f49ad3eeaee)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_matrix at ./show.jl:1407
print_matrix at ./show.jl:1379                                                                                                                                  [32/210]
unknown function (ip: 0x7f49ad3ea366)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
jl_apply at /home/centos/buildbot/slave/package_tarball64/build/src/julia.h:1392 [inlined]
jl_f__apply at /home/centos/buildbot/slave/package_tarball64/build/src/builtins.c:547
#showarray#330 at ./show.jl:1618
unknown function (ip: 0x7f49ad3e517d)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./REPL.jl:132
unknown function (ip: 0x7f49ad3e4cb6)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./REPL.jl:135
unknown function (ip: 0x7f49ad3e4986)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
display at ./multimedia.jl:143
unknown function (ip: 0x7f49ad3e47d2)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
print_response at ./REPL.jl:154
unknown function (ip: 0x7f49ad3e3ae8)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942                                                                            [7/210]
print_response at ./REPL.jl:139
unknown function (ip: 0x7f49ad3e3568)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
#22 at ./REPL.jl:652
unknown function (ip: 0x7f49ad3e1c81)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
run_interface at ./LineEdit.jl:1579
unknown function (ip: 0x7f4bb84270bf)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
run_frontend at ./REPL.jl:903
run_repl at ./REPL.jl:188
unknown function (ip: 0x7f49ad3df712)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
_start at ./client.jl:360
unknown function (ip: 0x7f4bb84422e8)
jl_call_method_internal at /home/centos/buildbot/slave/package_tarball64/build/src/julia_internal.h:189 [inlined]
jl_apply_generic at /home/centos/buildbot/slave/package_tarball64/build/src/gf.c:1942
jl_apply at /home/centos/buildbot/slave/package_tarball64/build/ui/../src/julia.h:1392 [inlined]
true_main at /home/centos/buildbot/slave/package_tarball64/build/ui/repl.c:112
main at /home/centos/buildbot/slave/package_tarball64/build/ui/repl.c:232
__libc_start_main at /lib64/libc.so.6 (unknown line)
unknown function (ip: 0x4013fc)
Allocations: 7528308 (Pool: 7526481; Big: 1827); GC: 12
Segmentation fault
5

It would be much easier to give advice if you would show the exact Fortran function signature you’re trying to call. Apart from that: an Array should already be a pointed to memory block, so no need for the Ref{...} - at least as long as it’s not allocatable. You may want to check that you really need Int64.

6

The fortran function is

      subroutine sqrinr(m,n,Q,ldq,R,ldr,j,x,w)                                                                                                                          
c purpose:      updates a QR factorization after inserting a new                                                                                                        
c               row.                                                                                                                                                    
c               i.e., given an m-by-m unitary matrix Q, an m-by-n                                                                                                       
c               upper trapezoidal matrix R and index j in the range                                                                                                     
c               1:m+1, this subroutine updates Q -> Q1  and R -> R1                                                                                                     
c               so that Q1 is again unitary, R1 upper trapezoidal,                                                                                                      
c               and Q1*R1 = [A(1:j-1,:); x; A(j:m,:)], where A = Q*R.                                                                                                   
c               (real version)                                                                                                                                          
c arguments:                                                                                                                                                            
c m (in)        number of rows of the matrix Q.                                                                                                                         
c n (in)        number of columns of the matrix R.                                                                                                                      
c Q (io)        on entry, the unitary matrix Q.                                                                                                                         
c               on exit, the updated matrix Q1.                                                                                                                         
c ldq (in)      leading dimension of Q. ldq >= m+1.                                                                                                                     
c R (io)        on entry, the original matrix R.                                                                                                                        
c               on exit, the updated matrix R1.                                                                                                                         
c ldr (in)      leading dimension of R. ldr >= m+1.                                                                                                                     
c j (in)        the position of the new row in R1                                                                                                                       
c x (io)        on entry, the row being added                                                                                                                           
c               on exit, x is destroyed.                                                                                                                                
c w (out)       a workspace vector of size min(m,n).

I removed the Ref{…} on the Arrays.

7

The signature including the types, which is the important information here, is

subroutine sqrinr(m,n,Q,ldq,R,ldr,j,x,w)
    integer m,n,ldq,ldr
    real Q(ldq,*),R(ldr,*),x(*),w(*)
end subroutine

Since this is Fortran77, and presuming you didn’t give other compiler flags to switch to 64bit default integer and float size, then the ccall should be along the lines of

Q = zeros( Float32, n+1, q )
R = zeros( Float32, n+1, q )
x = zeros( Float32, n)
w = zeros( Float32, min(n,q))
ccall( ("sqrinr_","qrupdate-1.1.2/libqrupdate")
     , Void    # nothing returned by Fortran subroutines
     , ( Ref{Int32}           # m
       , Ref{Int32}           # n
       , Ref{Float32}         # Q, not Array{Float32, 2}
       , Ref{Int32}           # ldq
       , Ref{Float32}         # R, not Array{Float32, 2}
       , Ref{Int32}           # ldr
       , Ref{Int32}           # j
       , Ref{Float32}         # x, not Array{Float32, 1}
       , Ref{Float32}         # w, not Array{Float32, 1}
       )
     , n, q, Q, n+1, R, n+1, n, xp, w)

Edit: Of course, Simon is correct and the array should also be passed as a Ref{T}. Changed the above lines accordingly.

1 Like
8

You should pass arrays as Ref{T} where T is the element type. Array should basically never appear in a ccall signature, as it is a Julia object (and no C libraries that I know of expect to be passed Julia arrays).

i.e. try:

ccall(("sqrinr_","qrupdate-1.1.2/libqrupdate"), Void, 
    (Ref{Int64}, Ref{Int64}, Ref{Float64}, Ref{Int64}, Ref{Float64}, Ref{Int64}, Ref{Int64}
, Ref{Float64}, Ref{Float64}), 
    n, q, Q, n+1, R, n+1, n, xp, w)
9

Also here is a simple example:

subroutine foo1(n,x,y)
        integer*4, intent(in) :: n
        integer*4, intent(in),  dimension(n) :: x
        integer*4, intent(out), dimension(n) :: y
        y = 2*x
end subroutine foo1

subroutine foo2(n,x,y)
        integer*8, intent(in) :: n
        integer*8, intent(in),  dimension(n) :: x
        integer*8, intent(out), dimension(n) :: y
        y = 2*x
end subroutine foo2

Compile into a shared library, and then in Julia, e.g.:

ftest = Libdl.dlopen("ftest.dll")
t = ones(Int32,2)
ccall(Libdl.dlsym(ftest,:foo1_), Void, (Ref{Int32}, Ref{Int32}, Ptr{Int32}), Int32(2), Int32.([3, -1]), t)
println(t)
t = ones(Int64,2)
ccall(Libdl.dlsym(ftest,:foo2_), Void, (Ref{Int64}, Ref{Int64}, Ptr{Int64}), 2, [3, -1], t)
println(t)
Libdl.dlclose(ftest)
10

Thanks guys. It seems to be running, but I’m not getting the results I want. This might be too much to ask but can anyone take a look at my code and find any obvious issues?

The code generates X a n by q matrix, beta=[1,2,3,4,5], and generates y ~ N( X*beta, 1)

Q,R is the QR decomposition of X. I want to use the fortran code to get the updated QR decompositon of X after adding another row (xp) to X. After I run the fortran code, the matrix Q and R are not correct.

using Distributions

n=20
p=1
q=5

X = zeros(Float64, n, q)
for i in 1:n
  X[i,:] = randn(q)
end
beta = [1,2,3,4,5]

y = rand(MvNormal(X * beta, 1.0))

Q, R = qr(X, thin=false)

## resize Q so that its n+1 by n+1, append zeros to R so its n+1 by q
Qp = vcat(hcat(Q, zeros(n)), zeros(n+1)')
Rp = vcat(R, reshape(zeros( (n-q+1)*q), n-q+1,q))

xp = randn(q)
yp = (xp'*beta)[1]
w = zeros(min(n,q))

ccall(("sqrinr_","libqrupdate"), Void,
      (Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Int32}, Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}),
      n, q, Qp, n+1, Rp, n+1, n+1, xp, w)
11

Here is another example for you to try:

subroutine foo3(x,y)
        integer*8, intent(in) :: x(:)
        integer*8, intent(out) :: y(:)
        y = 2*x
end subroutine foo3

This example shows that Fortran doesn’t know the size of the Julia array. The subroutine sqrinr uses arrays which have unspecified dimensions.

1 Like