fortrangslfortran-iso-c-bindingfgsl

Passing external function and self-defined data type through type(c_ptr) dummy variable


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);

Solution

  • 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?