I want to use the GSL with the Fortran interface FGSL in order to use the provided multiroots solver. In order to run the simulation I need to define the residual of the root finding problem. In the tests of FSGL there is this Rosenbrock example:
function rosenbrock_f(x, params, f) bind(c)
type(c_ptr), value :: x, params, f
integer(c_int) :: rosenbrock_f
!
type(fgsl_vector) :: fx, ff
real(fgsl_double), pointer :: par(:), xv(:), yv(:)
integer(fgsl_int) :: status
!
call fgsl_obj_c_ptr(fx, x)
call fgsl_obj_c_ptr(ff, f)
call c_f_pointer(params, par, (/ 2 /))
status = fgsl_vector_align(xv, fx)
status = fgsl_vector_align(yv, ff)
yv(1) = par(1) * (1.0_fgsl_double - xv(1))
yv(2) = par(2) * (xv(2) - xv(1)*xv(1))
rosenbrock_f = fgsl_success
end function rosenbrock_f
Here x is the ingoing variable, params are additional variables and f is the outgoing residual. In my code I already have a function for the residual, so I would like to pass it as an external function though params. Unfortunately, I also need to pass another self-defined data type, also through params. I need it this way, because the solver from GSL needs the function with the correct dummy variables and I do not know how to achieve this.
So my question is:
How can I pass an external function and a self-defined data-type through params in order to use it within this function?
Altering the above example led me to
function rosenbrock_f(x, params, f) bind(c)
type(c_ptr), value :: x, params, f
integer(c_int) :: rosenbrock_f
!
type(fgsl_vector) :: fx, ff
real(fgsl_double), pointer :: xv(:), yv(:) ! par removed
integer(fgsl_int) :: status
! NEW DECLARATIONS
type(owntype) :: this
external :: GetRes
!
call fgsl_obj_c_ptr(fx, x)
call fgsl_obj_c_ptr(ff, f)
! In this region I need to get GetRes (external function)
! and this (self-defined data type) from the params variable
! instead of par
call c_f_pointer(params, GetRes))
!call c_f_pointer(params, this)) ! not possible
status = fgsl_vector_align(xv, fx)
status = fgsl_vector_align(yv, ff)
call GetRes(this,xv,yv) ! call the external function
rosenbrock_f = fgsl_success
end function rosenbrock_f
However, this would only let me access GetRes. I have something in mind like a params pointer with shape 2. Similar like above, but allowing for different data types.
Edit 1 (some more code): Here is some more code to it. The external Function would look for the example above like this.
subroutine GetRes(this,y,res)
implicit none
type(owntype) :: this
double precision, dimension(this%neq) :: y
double precision, dimension(this%neq) :: res
intent (in) :: y
intent (out) :: res
res(1) = (1.0 - y(1))
res(2) = (y(2) - y(1)*y(1))
end subroutine GetRes
In order to see how the function rosenbrock_f is used I think the best is to reference the multiroots tests.
The solver is initialized by in line 80:
mroot_f = fgsl_multiroot_function_init(rosenbrock_f,nrt,ptr)
set in line 84
status = fgsl_multiroot_fsolver_set(mroot_fslv, mroot_f, xvec)
and then iteratively called in line 97:
status = fgsl_multiroot_fsolver_iterate(mroot_fslv);
Really more of a question than an answer, but...
Can you just define a type that contains all the C pointers that you need?
type, bind(C) :: T
type(C_PTR) this
type(C_FUNPTR) GetRes
end type T
Then you can create an instance of T
and point its this
and GetRes
pointers at what you wanted, then you have
type(C_PTR) par
!...
par = C_LOC(T_instance)
So now par
will be all loaded up and in your ultimate function you have the declarations
type(C_PTR), value :: params
type(T), pointer :: T_instance
type(owntype), pointer :: this
procedure(whatever), pointer :: GetRes
And then the sequence
call C_F_POINTER(params,T_instance)
call C_F_POINTER(T_instance%this,this)
call C_F_PROCPOINTER(T_instance%GetRes, GetRes)
Or do I have too simplistic an interpretation of your question?