diff --git a/src/gen/direct_wrappers/gsl_chebyshev_h.jl b/src/gen/direct_wrappers/gsl_chebyshev_h.jl index ea58fab..407e431 100644 --- a/src/gen/direct_wrappers/gsl_chebyshev_h.jl +++ b/src/gen/direct_wrappers/gsl_chebyshev_h.jl @@ -57,7 +57,7 @@ GSL documentation: > and requires $n$ function evaluations. """ -function cheb_init(cs, func, a, b) +function cheb_init(cs, func::gsl_function, a, b) ccall((:gsl_cheb_init, libgsl), Cint, (Ref{gsl_cheb_series}, Ref{gsl_function}, Cdouble, Cdouble), cs, func, a, b) end diff --git a/src/gen/direct_wrappers/gsl_deriv_h.jl b/src/gen/direct_wrappers/gsl_deriv_h.jl index a6d6e25..545c76f 100644 --- a/src/gen/direct_wrappers/gsl_deriv_h.jl +++ b/src/gen/direct_wrappers/gsl_deriv_h.jl @@ -32,7 +32,7 @@ GSL documentation: > actually used. """ -function deriv_central(f, x, h, result, abserr) +function deriv_central(f::F, x, h, result, abserr) where F ccall((:gsl_deriv_central, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, h, result, abserr) end @@ -58,7 +58,7 @@ GSL documentation: > negative step-size. """ -function deriv_backward(f, x, h, result, abserr) +function deriv_backward(f::F, x, h, result, abserr) where F ccall((:gsl_deriv_backward, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, h, result, abserr) end @@ -89,7 +89,7 @@ GSL documentation: > $x+h/2$, $x+h$. """ -function deriv_forward(f, x, h, result, abserr) +function deriv_forward(f::F, x, h, result, abserr) where F ccall((:gsl_deriv_forward, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, h, result, abserr) end diff --git a/src/gen/direct_wrappers/gsl_diff_h.jl b/src/gen/direct_wrappers/gsl_diff_h.jl index 97f99cf..6b885a4 100644 --- a/src/gen/direct_wrappers/gsl_diff_h.jl +++ b/src/gen/direct_wrappers/gsl_diff_h.jl @@ -12,7 +12,7 @@ C signature: `int gsl_diff_central (const gsl_function *f, double x, double *result, double *abserr)` """ -function diff_central(f, x, result, abserr) +function diff_central(f::F, x, result, abserr) where F ccall((:gsl_diff_central, libgsl), Cint, (Ref{gsl_function}, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, result, abserr) end @@ -22,7 +22,7 @@ end C signature: `int gsl_diff_backward (const gsl_function *f, double x, double *result, double *abserr)` """ -function diff_backward(f, x, result, abserr) +function diff_backward(f::F, x, result, abserr) where F ccall((:gsl_diff_backward, libgsl), Cint, (Ref{gsl_function}, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, result, abserr) end @@ -32,7 +32,7 @@ end C signature: `int gsl_diff_forward (const gsl_function *f, double x, double *result, double *abserr)` """ -function diff_forward(f, x, result, abserr) +function diff_forward(f::F, x, result, abserr) where F ccall((:gsl_diff_forward, libgsl), Cint, (Ref{gsl_function}, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, x, result, abserr) end diff --git a/src/gen/direct_wrappers/gsl_integration_h.jl b/src/gen/direct_wrappers/gsl_integration_h.jl index bbca63a..6b808f1 100644 --- a/src/gen/direct_wrappers/gsl_integration_h.jl +++ b/src/gen/direct_wrappers/gsl_integration_h.jl @@ -168,7 +168,7 @@ GSL documentation: > > texinfo > -> w(x) = sin(omega x) +> w(x) = sin(omega x) > w(x) = cos(omega x) > > The parameter `L` must be the length of the interval over which the @@ -253,7 +253,7 @@ end C signature: `void gsl_integration_qk15 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk15(f, a, b, result, abserr, resabs, resasc) +function integration_qk15(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk15, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -263,7 +263,7 @@ end C signature: `void gsl_integration_qk21 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk21(f, a, b, result, abserr, resabs, resasc) +function integration_qk21(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk21, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -273,7 +273,7 @@ end C signature: `void gsl_integration_qk31 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk31(f, a, b, result, abserr, resabs, resasc) +function integration_qk31(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk31, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -283,7 +283,7 @@ end C signature: `void gsl_integration_qk41 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk41(f, a, b, result, abserr, resabs, resasc) +function integration_qk41(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk41, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -293,7 +293,7 @@ end C signature: `void gsl_integration_qk51 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk51(f, a, b, result, abserr, resabs, resasc) +function integration_qk51(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk51, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -303,7 +303,7 @@ end C signature: `void gsl_integration_qk61 (const gsl_function * f, double a, double b, double *result, double *abserr, double *resabs, double *resasc)` """ -function integration_qk61(f, a, b, result, abserr, resabs, resasc) +function integration_qk61(f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk61, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, result, abserr, resabs, resasc) end @@ -313,7 +313,7 @@ end C signature: `void gsl_integration_qcheb (gsl_function * f, double a, double b, double *cheb12, double *cheb24)` """ -function integration_qcheb(f, a, b, cheb12, cheb24) +function integration_qcheb(f::F, a, b, cheb12, cheb24) where F ccall((:gsl_integration_qcheb, libgsl), Cvoid, (Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}), f, a, b, cheb12, cheb24) end @@ -323,7 +323,7 @@ end C signature: `void gsl_integration_qk (const int n, const double xgk[], const double wg[], const double wgk[], double fv1[], double fv2[], const gsl_function *f, double a, double b, double * result, double * abserr, double * resabs, double * resasc)` """ -function integration_qk(n, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc) +function integration_qk(n, xgk, wg, wgk, fv1, fv2, f::F, a, b, result, abserr, resabs, resasc) where F ccall((:gsl_integration_qk, libgsl), Cvoid, (Cint, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{gsl_function}, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}), n, xgk, wg, wgk, fv1, fv2, f, a, b, result, abserr, resabs, resasc) end @@ -348,7 +348,7 @@ GSL documentation: > of function evaluations. """ -function integration_qng(f, a, b, epsabs, epsrel, result, abserr, neval) +function integration_qng(f::F, a, b, epsabs, epsrel, result, abserr, neval) where F ccall((:gsl_integration_qng, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Csize_t}), f, a, b, epsabs, epsrel, result, abserr, neval) end @@ -394,7 +394,7 @@ GSL documentation: > allocated size of the workspace. """ -function integration_qag(f, a, b, epsabs, epsrel, limit, key, workspace, result, abserr) +function integration_qag(f::F, a, b, epsabs, epsrel, limit, key, workspace, result, abserr) where F ccall((:gsl_integration_qag, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Csize_t, Cint, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, epsabs, epsrel, limit, key, workspace, result, abserr) end @@ -420,7 +420,7 @@ GSL documentation: > In this case a lower-order rule is more efficient. """ -function integration_qagi(f, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qagi(f::F, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qagi, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, epsabs, epsrel, limit, workspace, result, abserr) end @@ -443,7 +443,7 @@ GSL documentation: > and then integrated using the QAGS algorithm. """ -function integration_qagiu(f, a, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qagiu(f::F, a, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qagiu, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, a, epsabs, epsrel, limit, workspace, result, abserr) end @@ -466,7 +466,7 @@ GSL documentation: > and then integrated using the QAGS algorithm. """ -function integration_qagil(f, b, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qagil(f::F, b, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qagil, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, b, epsabs, epsrel, limit, workspace, result, abserr) end @@ -493,7 +493,7 @@ GSL documentation: > which may not exceed the allocated size of the workspace. """ -function integration_qags(f, a, b, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qags(f::F, a, b, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qags, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, epsabs, epsrel, limit, workspace, result, abserr) end @@ -527,7 +527,7 @@ GSL documentation: > region then this routine will be faster than `gsl_integration_qags`. """ -function integration_qagp(f, pts, npts, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qagp(f::F, pts, npts, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qagp, libgsl), Cint, (Ref{gsl_function}, Ref{Cdouble}, Csize_t, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, pts, npts, epsabs, epsrel, limit, workspace, result, abserr) end @@ -547,7 +547,7 @@ GSL documentation: > not texinfo > > $$I = \int_a^b dx\, {f(x) \over x - c} -> = \lim_{\epsilon \to 0} +> = \lim_{\epsilon \to 0} > \left\{ > \int_a^{c-\epsilon} dx\, {f(x) \over x - c} > + @@ -566,7 +566,7 @@ GSL documentation: > ordinary 15-point Gauss-Kronrod integration rule. """ -function integration_qawc(f, a, b, c, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qawc(f::F, a, b, c, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qawc, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, c, epsabs, epsrel, limit, workspace, result, abserr) end @@ -595,7 +595,7 @@ GSL documentation: > Gauss-Kronrod integration rule is used. """ -function integration_qaws(f, a, b, t, epsabs, epsrel, limit, workspace, result, abserr) +function integration_qaws(f::F, a, b, t, epsabs, epsrel, limit, workspace, result, abserr) where F ccall((:gsl_integration_qaws, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Ref{gsl_integration_qaws_table}, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{Cdouble}, Ref{Cdouble}), f, a, b, t, epsabs, epsrel, limit, workspace, result, abserr) end @@ -627,7 +627,7 @@ GSL documentation: > > texinfo > -> I = int\_a^b dx f(x) sin(omega x) +> I = int\_a^b dx f(x) sin(omega x) > I = int\_a^b dx f(x) cos(omega x) > > The results are extrapolated using the epsilon-algorithm to accelerate @@ -645,7 +645,7 @@ GSL documentation: > integration. """ -function integration_qawo(f, a, epsabs, epsrel, limit, workspace, wf, result, abserr) +function integration_qawo(f::F, a, epsabs, epsrel, limit, workspace, wf, result, abserr) where F ccall((:gsl_integration_qawo, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}), f, a, epsabs, epsrel, limit, workspace, wf, result, abserr) end @@ -737,7 +737,7 @@ GSL documentation: > `cycle_workspace` as workspace for the QAWO algorithm. """ -function integration_qawf(f, a, epsabs, limit, workspace, cycle_workspace, wf, result, abserr) +function integration_qawf(f::F, a, epsabs, limit, workspace, cycle_workspace, wf, result, abserr) where F ccall((:gsl_integration_qawf, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Csize_t, Ref{gsl_integration_workspace}, Ref{gsl_integration_workspace}, Ref{gsl_integration_qawo_table}, Ref{Cdouble}, Ref{Cdouble}), f, a, epsabs, limit, workspace, cycle_workspace, wf, result, abserr) end @@ -793,7 +793,7 @@ GSL documentation: > table `t` and returns the result. """ -function integration_glfixed(f, a, b, t) +function integration_glfixed(f::F, a, b, t) where F ccall((:gsl_integration_glfixed, libgsl), Cdouble, (Ref{gsl_function}, Cdouble, Cdouble, Ref{gsl_integration_glfixed_table}), f, a, b, t) end @@ -890,7 +890,7 @@ GSL documentation: > set to `NULL`. """ -function integration_cquad(f, a, b, epsabs, epsrel, ws, result, abserr, nevals) +function integration_cquad(f::F, a, b, epsabs, epsrel, ws, result, abserr, nevals) where F ccall((:gsl_integration_cquad, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Ref{gsl_integration_cquad_workspace}, Ref{Cdouble}, Ref{Cdouble}, Ref{Csize_t}), f, a, b, epsabs, epsrel, ws, result, abserr, nevals) end @@ -957,7 +957,7 @@ GSL documentation: > `neval`. """ -function integration_romberg(f, a, b, epsabs, epsrel, result, neval, w) +function integration_romberg(f::F, a, b, epsabs, epsrel, result, neval, w) where F ccall((:gsl_integration_romberg, libgsl), Cint, (Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Csize_t}, Ref{gsl_integration_romberg_workspace}), f, a, b, epsabs, epsrel, result, neval, w) end @@ -1116,7 +1116,6 @@ GSL documentation: > approximated as """ -function integration_fixed(func, result, w) +function integration_fixed(func::F, result, w) where F ccall((:gsl_integration_fixed, libgsl), Cint, (Ref{gsl_function}, Ref{Cdouble}, Ref{gsl_integration_fixed_workspace}), func, result, w) end - diff --git a/src/gen/direct_wrappers/gsl_min_h.jl b/src/gen/direct_wrappers/gsl_min_h.jl index d2bb693..f7ed0a2 100644 --- a/src/gen/direct_wrappers/gsl_min_h.jl +++ b/src/gen/direct_wrappers/gsl_min_h.jl @@ -67,7 +67,7 @@ GSL documentation: > returns an error code of `GSL_EINVAL`. """ -function min_fminimizer_set(s, f, x_minimum, x_lower, x_upper) +function min_fminimizer_set(s, f::gsl_function, x_minimum, x_lower, x_upper) ccall((:gsl_min_fminimizer_set, libgsl), Cint, (Ref{gsl_min_fminimizer}, Ref{gsl_function}, Cdouble, Cdouble, Cdouble), s, f, x_minimum, x_lower, x_upper) end @@ -86,7 +86,7 @@ GSL documentation: > `f(x_minimum)`, `f(x_lower)` and `f(x_upper)`. """ -function min_fminimizer_set_with_values(s, f, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper) +function min_fminimizer_set_with_values(s, f::gsl_function, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper) ccall((:gsl_min_fminimizer_set_with_values, libgsl), Cint, (Ref{gsl_min_fminimizer}, Ref{gsl_function}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), s, f, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper) end @@ -301,7 +301,7 @@ end C signature: `int gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum, double *x_lower, double * f_lower, double *x_upper, double * f_upper, size_t eval_max)` """ -function min_find_bracket(f, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper, eval_max) +function min_find_bracket(f::F, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper, eval_max) where F ccall((:gsl_min_find_bracket, libgsl), Cint, (Ref{gsl_function}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Csize_t), f, x_minimum, f_minimum, x_lower, f_lower, x_upper, f_upper, eval_max) end diff --git a/src/gen/direct_wrappers/gsl_roots_h.jl b/src/gen/direct_wrappers/gsl_roots_h.jl index a835e33..68a497b 100644 --- a/src/gen/direct_wrappers/gsl_roots_h.jl +++ b/src/gen/direct_wrappers/gsl_roots_h.jl @@ -66,7 +66,7 @@ GSL documentation: > `x_upper`\]. """ -function root_fsolver_set(s, f, x_lower, x_upper) +function root_fsolver_set(s, f::gsl_function, x_lower, x_upper) ccall((:gsl_root_fsolver_set, libgsl), Cint, (Ref{gsl_root_fsolver}, Ref{gsl_function}, Cdouble, Cdouble), s, f, x_lower, x_upper) end diff --git a/src/manual_wrappers.jl b/src/manual_wrappers.jl index c73e789..2a590e2 100644 --- a/src/manual_wrappers.jl +++ b/src/manual_wrappers.jl @@ -1,6 +1,6 @@ # # Manual wrappers for the gsl_* stuff -# +# export wrap_gsl_vector, wrap_gsl_matrix @@ -24,12 +24,36 @@ Return a Julia matrix wrapping the data of a gsl_matrix """ @inline function wrap_gsl_matrix(m::Ptr{gsl_matrix}) M = unsafe_load(m) - @assert M.size2==M.tda "Cannot unsafe_wrap gsl_matrix with tda != size2." + @assert M.size2==M.tda "Cannot unsafe_wrap gsl_matrix with tda != size2." return unsafe_wrap(Array{Float64}, M.data, (M.size1, M.size2)) end ## Root finding +gsl_function_helper(x::Cdouble, (f,))::Cdouble = f(x) + +# The following code relies on `gsl_function` being a mutable type +# (such that we can call `pointer_from_objref` on it) to simplify the object structure +# a little bit and avoid hitting some limitation of the allocation optimizer. +@assert ismutable(gsl_function(C_NULL, C_NULL)) + +function Base.cconvert(::Type{Ref{gsl_function}}, t::T) where {F,T<:Tuple{F}} + # We need to allocate the `gsl_function` here to be kept alive by ccall + # This require us to create the pointer to the function and the callable object + param_ref = Base.cconvert(Ref{T}, t) + fptr = @cfunction(gsl_function_helper, Cdouble, (Cdouble, Ref{T})) + param_ptr = Base.unsafe_convert(Ref{T}, param_ref) + gsl_func = gsl_function(fptr, param_ptr) + return gsl_func, param_ref +end +function Base.unsafe_convert(::Type{Ref{gsl_function}}, + (gsl_func,)::Tuple{gsl_function, F}) where F + return pointer_from_objref(gsl_func) +end + +Base.cconvert(::Type{Ref{gsl_function}}, gslf::gsl_function) = + convert(Ref{gsl_function}, gslf) + # Macros for easier creation of gsl_function and gsl_function_fdf structs export @gsl_function, @gsl_function_fdf @@ -47,6 +71,34 @@ macro gsl_function(f) ) end +gsl_function_f_helper(x::Cdouble, (f,))::Cdouble = f(x) +gsl_function_df_helper(x::Cdouble, (f,df,))::Cdouble = df(x) +gsl_function_fdf_helper(x::Cdouble, (f,df,fdf))::Tuple{Cdouble,Cdouble} = fdf(x) + +@assert ismutable(gsl_function_fdf(C_NULL, C_NULL, C_NULL, C_NULL)) + +function Base.cconvert(::Type{Ref{gsl_function_fdf}}, t::T) where {F,DF,FDF,T<:Tuple{F,DF,FDF}} + # We need to allocate the `gsl_function_fdf` here to be kept alive by ccall + # This require us to create the pointer to the function and the callable object + param_ref = Base.cconvert(Ref{T}, t) + fptr = @cfunction(gsl_function_f_helper, Cdouble, (Cdouble, Ref{T})) + dfptr = @cfunction(gsl_function_df_helper, Cdouble, (Cdouble, Ref{T})) + fdfptr = @cfunction(gsl_function_fdf_helper, Tuple{Cdouble,Cdouble}, (Cdouble, Ref{T})) + param_ptr = Base.unsafe_convert(Ref{T}, param_ref) + gsl_func = gsl_function_fdf(fptr, dfptr, fdfptr, param_ptr) + return gsl_func, param_ref +end +Base.cconvert(T::Type{Ref{gsl_function_fdf}}, (f,df)::Tuple{F,DF}) where {F,DF} = + Base.cconvert(T, (f, df, x -> (f(x), df(x)))) +function Base.unsafe_convert(::Type{Ref{gsl_function_fdf}}, + (gsl_func,)::Tuple{gsl_function_fdf, F}) where F + return pointer_from_objref(gsl_func) +end + + +Base.cconvert(::Type{Ref{gsl_function_fdf}}, gslf::gsl_function_fdf) = + convert(Ref{gsl_function_fdf}, gslf) + """ @gsl_function_fdf(f, df, fdf) @@ -105,6 +157,32 @@ end export @gsl_multiroot_function, @gsl_multiroot_function_fdf +function gsl_multiroot_function_helper(x_vec::Ptr{gsl_vector}, (f,), y_vec::Array{gsl_vector})::Cint + x = GSL.wrap_gsl_vector(x_vec) + y = GSL.wrap_gsl_vector(y_vec) + f(x, y) + return Cint(GSL.GSL_SUCCESS) +end + +@assert ismutable(gsl_multiroot_function(C_NULL, 0, C_NULL)) + +function Base.cconvert(::Type{Ref{gsl_multiroot_function}}, (f,n)::Tuple{F,Integer}) where F + # We need to allocate the `gsl_function_fdf` here to be kept alive by ccall + # This require us to create the pointer to the function and the callable object + param_ref = Base.cconvert(Ref{Tuple{F}}, (f,)) + fptr = @cfunction(gsl_multiroot_function_helper, Cint, (Ptr{gsl_vector}, Ref{Tuple{F}}, Ptr{gsl_vector})) + param_ptr = Base.unsafe_convert(Ref{Tuple{F}}, param_ref) + gsl_func = gsl_multiroot_function(fptr, n, param_ptr) + return gsl_func, param_ref +end +function Base.unsafe_convert(::Type{Ref{gsl_multiroot_function}}, + (gsl_func,)::Tuple{gsl_multiroot_function, F}) where F + return pointer_from_objref(gsl_func) +end + +Base.cconvert(::Type{Ref{gsl_multiroot_function}}, gslf::gsl_multiroot_function) = + convert(Ref{gsl_multiroot_function}, gslf) + """ @gsl_multiroot_function(f, n) @@ -121,7 +199,7 @@ macro gsl_multiroot_function(f, n) x = GSL.wrap_gsl_vector(x_vec) y = GSL.wrap_gsl_vector(y_vec) $f(x, y) - return Cint(GSL.GSL_SUCCESS) + return Cint(GSL.GSL_SUCCESS) end, Cint, (Ptr{gsl_vector}, Ptr{Cvoid}, Ptr{gsl_vector})), # n @@ -132,6 +210,52 @@ macro gsl_multiroot_function(f, n) ) end + +function gsl_multiroot_function_f_helper(x_vec::Ptr{gsl_vector}, (f,), y_vec::Ptr{gsl_vector})::Cint + x = GSL.wrap_gsl_vector(x_vec) + y = GSL.wrap_gsl_vector(y_vec) + f(x, y) + return Cint(GSL.GSL_SUCCESS) +end +function gsl_multiroot_function_df_helper(x_vec::Ptr{gsl_vector}, (f,df), J_mat::Ptr{gsl_matrix})::Cint + x = GSL.wrap_gsl_vector(x_vec) + J = GSL.wrap_gsl_matrix(J_mat) + df(x, J) + return Cint(GSL.GSL_SUCCESS) +end +function gsl_multiroot_function_fdf_helper(x_vec::Ptr{gsl_vector}, (f,df,fdf), y_vec::Ptr{gsl_vector}, J_mat::Ptr{gsl_matrix})::Cint + x = GSL.wrap_gsl_vector(x_vec) + y = GSL.wrap_gsl_vector(y_vec) + J = GSL.wrap_gsl_matrix(J_mat) + fdf(x, y, J) + return Cint(GSL.GSL_SUCCESS) +end + +@assert ismutable(gsl_multiroot_function_fdf(C_NULL, C_NULL, C_NULL, 0, C_NULL)) + +function Base.cconvert(::Type{Ref{gsl_multiroot_function_fdf}}, (f,df,fdf,n)::T) where {F,DF,FDF,T<:Tuple{F,DF,FDF,Integer}} + # We need to allocate the `gsl_function_fdf` here to be kept alive by ccall + # This require us to create the pointer to the function and the callable object + param_ref = Base.cconvert(Ref{Tuple{F,DF,FDF}}, (f,df,fdf)) + fptr = @cfunction(gsl_multiroot_function_f_helper, Cint, (Ptr{gsl_vector}, Ref{Tuple{F,DF,FDF}}, Ptr{gsl_vector})) + dfptr = @cfunction(gsl_multiroot_function_df_helper, Cint, (Ptr{gsl_vector}, Ref{Tuple{F,DF,FDF}}, Ptr{gsl_matrix})) + fdfptr = @cfunction(gsl_multiroot_function_fdf_helper, Cint, (Ptr{gsl_vector}, Ref{Tuple{F,DF,FDF}}, Ptr{gsl_vector}, Ptr{gsl_matrix})) + param_ptr = Base.unsafe_convert(Ref{Tuple{F,DF,FDF}}, param_ref) + gsl_func = gsl_multiroot_function_fdf(fptr, dfptr, fdfptr, n, param_ptr) + return gsl_func, param_ref +end +Base.cconvert(T::Type{Ref{gsl_multiroot_function_fdf}}, (f,df,n)::Tuple{F,DF,Integer}) where {F,DF} = + Base.cconvert(T, (f, df, (x, y, J) -> (f(x, y); df(x, J)), n)) +function Base.unsafe_convert(::Type{Ref{gsl_multiroot_function_fdf}}, + (gsl_func,)::Tuple{gsl_multiroot_function_fdf, F}) where F + return pointer_from_objref(gsl_func) +end + + +Base.cconvert(::Type{Ref{gsl_multiroot_function_fdf}}, gslf::gsl_multiroot_function_fdf) = + convert(Ref{gsl_multiroot_function_fdf}, gslf) + + """ @gsl_multiroot_function_fdf(f, df, n) @@ -225,7 +349,7 @@ macro gsl_multiroot_function_fdf(f, df, fdf, n) 0 ) ) -end +end ## Hypergeometric function wrappers from original GSL.jl #(c) 2013 Jiahao Chen @@ -237,10 +361,10 @@ export hypergeom, hypergeom_e """ hypergeom(a, b, x::Float64) -> Float64 -Computes the appropriate hypergeometric ``{}_p F_q`` function, -where ``p`` and ``p`` are the lengths of the input vectors `a` and `b` respectively. +Computes the appropriate hypergeometric ``{}_p F_q`` function, +where ``p`` and ``p`` are the lengths of the input vectors `a` and `b` respectively. -Singleton `a` and/or `b` may be specified as scalars, +Singleton `a` and/or `b` may be specified as scalars, and length-0 `a` and/or `b` may be input as simply `[]`. Supported values of ``(p, q)`` are ``(0, 0)``, ``(0, 1)``, ``(1, 1)``, ``(2, 0)`` and ``(2, 1)``. @@ -266,10 +390,10 @@ end """ hypergeom_e(a, b, x::Float64) -> gsl_sf_result -Computes the appropriate hypergeometric ``{}_p F_q`` function, -where ``p`` and ``p`` are the lengths of the input vectors `a` and `b` respectively. +Computes the appropriate hypergeometric ``{}_p F_q`` function, +where ``p`` and ``p`` are the lengths of the input vectors `a` and `b` respectively. -Singleton `a` and/or `b` may be specified as scalars, +Singleton `a` and/or `b` may be specified as scalars, and length-0 `a` and/or `b` may be input as simply `[]`. Supported values of ``(p, q)`` are ``(0, 0)``, ``(0, 1)``, ``(1, 1)``, ``(2, 0)`` and ``(2, 1)``. diff --git a/test/numdiff.jl b/test/numdiff.jl index eb1c445..bc52044 100644 --- a/test/numdiff.jl +++ b/test/numdiff.jl @@ -1,7 +1,7 @@ using GSL using Test -func(x) = x^3 +func(x) = x^3 @testset "Numerical differentiation" begin x = 1.0 @@ -14,6 +14,7 @@ func(x) = x^3 df,ddf = Cdouble[0], Cdouble[0] deriv(@gsl_function(func), x, h, df, ddf) @test abs(df_dx - df[]) <= ddf[] <= d2f_dx2*h/2 + 2eps()/h + deriv((func,), x, h, df, ddf) + @test abs(df_dx - df[]) <= ddf[] <= d2f_dx2*h/2 + 2eps()/h end end - diff --git a/test/quadrature.jl b/test/quadrature.jl index b428412..69886f9 100644 --- a/test/quadrature.jl +++ b/test/quadrature.jl @@ -17,5 +17,10 @@ fquad = x -> x^1.5 @test abs(result[] - 1/2.5) < abserr[] + fquad2 = let n=1.5 + x -> x^n + end + integration_cquad((fquad2,), a, b, 1e-10, 1e-10, ws, result, abserr, nevals) + integration_cquad_workspace_free(ws) end