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 GaussLegendre 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, 1e12, 1e12, 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 asRef{GSL_Function}
;  the workspace
w
allocated by the C functiongsl_integration_workspace_alloc
can be passed asPtr{Void}
,  the type assertion in
gsl_function_wrap
is requiredjfunc = 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.