C callback for GSL

Specific Domains Numerics
1

I would like to add fixed points quadrature to GSL.jl and to understand how callbacks for GSL work, I read the blog post as linked in the documentation.

Here is the code for Gauss-Legendre adaptive quadrature routine that includes the changes to reflect the recent PR by @yuyichao:

function gsl_function_wrap(x::Cdouble, gsl_func::Ptr{Void})
    jfunc = unsafe_pointer_to_objref(gsl_func)::Function
    convert(Cdouble, jfunc(x))::Cdouble
end


struct GSL_Function
    cfunc::Ptr{Void}
    jfunc::Any
    function GSL_Function(jfunc::Function)
        cfunc = cfunction(gsl_function_wrap, Cdouble, (Cdouble, Ptr{Void}))
        new(cfunc, jfunc)
    end
end


function integration_qag(julia_func::Function, a::Real, b::Real)
    w = ccall((:gsl_integration_workspace_alloc, "libgsl"), Ptr{Void}, (Csize_t,), 10^7)
    result = Ref{Cdouble}(0)
    abserr = Ref{Cdouble}(0)
    errno = ccall((:gsl_integration_qag, "libgsl"), Cint,
                  (Ref{GSL_Function}, Cdouble, Cdouble, Cdouble, Cdouble, Csize_t,
                   Cint, Ptr{Void}, Ref{Cdouble}, Ref{Cdouble}),
                  GSL_Function(julia_func), a, b, 1e-12, 1e-12, 10^3,
                  1, w, result, abserr)
    ccall((:gsl_integration_workspace_free, "libgsl"), Void, (Ptr{Void},), w)
    return result[], abserr[]
end

Note that ccall((:gsl_integration_qag,...), ...) in the blog post is missing the epsabs argument.

Now for the question: the code above seems to work, for instance try

f(x)=x^2
integration_qag(f, -2.0, 1.0)

but I would like to make sure this is correct. In particular, I understand that

  • GSL_Function(julia_func) should be passed as Ref{GSL_Function};
  • the workspace w allocated by the C function gsl_integration_workspace_alloc can be passed as Ptr{Void},
  • the type assertion in gsl_function_wrap is required jfunc = unsafe_pointer_to_objref(gsl_func)::Function and
  • putting const cfunc = cfunction(gsl_function_wrap, Cdouble, (Cdouble, Ptr{Void})) in the global scope returns a NULL pointer if the code above is in a module.

The blog post should be updated (I find it too useful to be removed from the documentation) since it uses the & construct that is deprecated and passes Julia objects as Ptr when Ref instead should be used. I can submit a PR, if somebody can confirm the code above is correct (I will add the other missing parameters to the Julia wrapper) and if the syntax has settled.

2

I should add that before writing the post, I was too missing the extra parameter abserr and getting segfaults, which disappeared once I used the correct C function signature.