Julia: Calling GSL functions

Question

I'm trying to integrate a simple function in Julia using the GSL numerical integration function QNG. The form of the function in C is given as

int gsl_integration_qng (const gsl_function * f, double a, double b,
 double epsabs, double epsrel, double * result, double * abserr, size_t * neval)

I am trying to perform the following integration

f(x) = x^2
a, b, epsabs, epsrel = 0.0, 1.0, 1.0e-2, 1.0e-2
result, abserr = 0.0, 0.0
neval =  0x123456789abcdef
0x0123456789abcdef
t = ccall( (:gsl_integration_qng, "libgsl"), Int32, 
(Ptr{Void}, Float64, Float64, Float64, Float64, Ptr{Float64}, Ptr{Float64},Csize_t), 
  &f, a, b, epsabs, epsrel, &result, &abserr, neval)

I appreciate any assitance.

Answer 1

After a bit of thought and looking at the existing GSL functions I got this working.

type gsl_function
    function_::Ptr{Void}
    params::Ptr{Void}
end

function function_callback(x::Cdouble, f_::Ptr{Void})
    f = unsafe_pointer_to_objref(f_)::Function
    convert(Cdouble, f(x))::Cdouble
end

const function_callback_ptr = cfunction(function_callback, Cdouble, (Cdouble, Ptr{Void}))

function integration_qng(f,a,b,epsrel,epsabs)
  result = convert(Ptr{Cdouble}, Array(Cdouble, 1))
  abserr = convert(Ptr{Cdouble}, Array(Cdouble, 1))
  neval = convert(Ptr{Csize_t}, Array(Csize_t, 1))
  t = ccall( (:gsl_integration_qng, :libgsl), Cint, (Ptr{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Csize_t}), 
  &gsl_function(function_callback_ptr, pointer_from_objref(f)), a, b, epsabs, epsrel, result, abserr, neval )
  return unsafe_load(result)
end

f(x) = x^2
a = 0
b = 1
epsabs = 0
epsrel = 1e-2

println(integration_qng(f,a,b,epsrel,epsabs))

Also see here http://julialang.org/blog/2013/05/callback/ for a discussion about this kind of thing from someone more qualified than I.

Answer 2

I wrote a blog post awhile ago on how to do this on a different GSL integration function.

The gist of it is:

Define a parametrically typed function

function integrand{T,T2,T3}(x::T,dim::T2,params::T3)
  A = 1.0 / (pi * pi * pi)
  return A / (1.0 - cos(unsafe_load(x,1))*cos(unsafe_load(x,2))*cos(unsafe_load(x,3)))::Cdouble
end

Then get the pointer

integrand_c = cfunction(integrand,Cdouble,(Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}))

and then ccall the GSL function with the function pointer as one of the arguments. In my example:

ccall((:monte_carlo_integrate,"/home/crackauc/Public/libMonte.so"),Int32,(Ptr{Void},Ptr{Cdouble},Ptr{Cdouble},Int32,Int32,Ptr{Cdouble},Ptr{Cdouble},Csize_t),integrand_c,x,y,mode,dim,xl,xu,calls)