Help with using @cfunction and arrays

Hi, I try to do Julia C API embedding, and I struggle with using @cfunction to work with arguments of the array type. I get a segmentation fault that is not very clear to me.

MWE:

  1. Julia code jl_qeq_solver.jl:
module QeqSolver
export solve

function solve!(a, b, c, roots)
    println("roots[1] = "),
    println(roots[1])
    if a == 0
        roots[1] = -c / b
        roots[2] = -c / b
    else
        D = b^2 - 4 * a * c
        println("D = ", D)
        if b > 0
            roots[1] = (-b - sqrt(D)) / (2 * a)
            roots[2] = c / (a * roots[1])
        else
            roots[1] = (-b + sqrt(D)) / (2 * a)
            roots[2] = c / (a * roots[1])
        end
    end

    return 0
end
end
  1. C code call_julia.c:
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <julia.h>

typedef int64_t (solve_fn)(double, double, double, double *);

int main(int argc, char *argv[])
{
    jl_init();
    double a = 1.0;
    double b = 5.0;
    double c = 4.0;
    double roots[2];

    jl_eval_string("include(\"jl_qeq_solver.jl\")");
    jl_eval_string("using .QeqSolver");

    jl_value_t *v = jl_eval_string("@cfunction(QeqSolver.solve!, Int, (Float64, Float64, Float64, Ref{Float64}))");
    solve_fn *solve = jl_unbox_voidpointer(v);

    int status = solve(a, b, c, roots);
    assert(status == 0);

    printf("roots[0] = %f, roots[1] = %f\n", roots[0], roots[1]);

    jl_atexit_hook(0);
    return 0;
}
  1. Makefile:
JL_SHARE = $(shell julia -e 'print(joinpath(Sys.BINDIR, Base.DATAROOTDIR, "julia"))')

CFLAGS   += $(shell $(JL_SHARE)/julia-config.jl --cflags)
CXXFLAGS += $(shell $(JL_SHARE)/julia-config.jl --cflags)
LDFLAGS  += $(shell $(JL_SHARE)/julia-config.jl --ldflags)
LDLIBS   += $(shell $(JL_SHARE)/julia-config.jl --ldlibs)

.PHONY : all
all: call_julia

The error that I get after (successful) compilation and running the executable ./call_julia:

roots[1] =
0.0
D = 9.0
fatal: error thrown and no exception handler available.
MethodError(f=Base.setindex!, args=(0, -4, 1), world=0x0000000000007ae2)
jl_method_error_bare at /cache/build/builder-amdci5-1/julialang/julia-release-1-dot-10/src/gf.c:2208
jl_method_error at /cache/build/builder-amdci5-1/julialang/julia-release-1-dot-10/src/gf.c:2226
jl_lookup_generic_ at /cache/build/builder-amdci5-1/julialang/julia-release-1-dot-10/src/gf.c:3057 [inlined]
ijl_apply_generic at /cache/build/builder-amdci5-1/julialang/julia-release-1-dot-10/src/gf.c:3072
solve! at /home/dima/dev/learn/2024-04-09-julia-embed-mwe/jl_qeq_solver.jl:14
unknown function (ip: 0x75a524e82422)
main at ./call_julia (unknown line)
unknown function (ip: 0x75a526229d8f)
__libc_start_main at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
_start at ./call_julia (unknown line)

I would highly appreciate if somebody could point me to what is wrong with my approach.

Ref{Float64} is not the correct type for roots. It looks like this should be Vector{Float64}.

On the julia side, you might do as follows.

julia> roots = zeros(2)
2-element Vector{Float64}:
 0.0
 0.0

julia> solve!(1,2,1, roots)
roots[1] =
0.0
D = 0
0

julia> roots
2-element Vector{Float64}:
 -1.0
 -1.0

If you used a Ref{Float64}, you will get an error.

julia> roots = Ref(0.0)
Base.RefValue{Float64}(0.0)

julia> solve!(1,2,1, roots)
roots[1] =
ERROR: MethodError: no method matching getindex(::Base.RefValue{Float64}, ::Int64)

Closest candidates are:
  getindex(::Ref, ::CartesianIndex{0})
   @ Base multidimensional.jl:1932
  getindex(::Base.RefValue)
   @ Base refvalue.jl:59

Stacktrace:
 [1] solve!(a::Int64, b::Int64, c::Int64, roots::Base.RefValue{Float64})
   @ Main ./REPL[1]:2
 [2] top-level scope
   @ REPL[6]:1

The next issue is how to pass roots to Julia from C. You need a version of solve! that will take a pointer and create a Vector.

function solve!(a, b, c, roots::Ptr{Float64})
    roots_vec = unsafe_wrap(Array, roots, 2)
    return solve!(a, b, c, roots_vec)
end

Edit: Fixed second method of solve! to use the argument roots

Dear @mkitti thank you very much! It works :smiley: