Передача внешней функции и самоопределяемого типа данных через фиктивную переменную типа (c_ptr) - PullRequest
1 голос
/ 29 марта 2020

Я хочу использовать GSL с интерфейсом Fortran FGSL , чтобы использовать предоставленный multiroots solver . Чтобы запустить симуляцию, мне нужно определить остаток от root проблемы с поиском. В тестах FSGL есть этот пример Розенброка :

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

Здесь x - входящая переменная, params - дополнительные переменные, а f - исходящий остаток. В моем коде у меня уже есть функция для остатка, поэтому я хотел бы передать ее как внешнюю функцию через params. К сожалению, мне также нужно передать другой самоопределенный тип данных, также через params. Мне это нужно, потому что решателю из GSL нужна функция с правильными фиктивными переменными, и я не знаю, как этого добиться.

Итак, мой вопрос :
Как я могу передать через функцию внешнюю функцию и самоопределяемый тип данных, чтобы использовать ее в этой функции?

Изменение приведенного выше примера привело меня к

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

Однако, это только позволило бы мне получить доступ к GetRes. Я имею в виду что-то вроде указателя параметров с формой 2. Аналогично вышеописанному, но с учетом разных типов данных.

Редактировать 1 (еще немного кода): Вот еще немного кода для Это. Внешняя функция будет выглядеть для примера выше, как это.

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

Чтобы увидеть, как используется функция rosenbrock_f, я думаю, лучше всего сослаться на многокорневые тесты .

Решатель инициализируется в строка 80 :

mroot_f = fgsl_multiroot_function_init(rosenbrock_f,nrt,ptr)

устанавливается в строку 84

status = fgsl_multiroot_fsolver_set(mroot_fslv, mroot_f, xvec)

и затем итеративно вызывается в строка 97 :

status = fgsl_multiroot_fsolver_iterate(mroot_fslv);

1 Ответ

2 голосов
/ 29 марта 2020

Действительно больше вопрос, чем ответ, но ...

Можете ли вы просто определить тип, который содержит все нужные вам C указатели?

type, bind(C) :: T
   type(C_PTR) this
   type(C_FUNPTR) GetRes
end type T

Тогда Вы можете создать экземпляр T и указать его указатели this и GetRes на то, что вы хотели, тогда у вас есть

type(C_PTR) par
!...
par = C_LOC(T_instance)

Так что теперь par будет загружен и в вашем Конечная функция у вас есть декларации

type(C_PTR), value :: params
type(T), pointer :: T_instance
type(owntype), pointer :: this
procedure(whatever), pointer :: GetRes

А затем последовательность

call C_F_POINTER(params,T_instance)
call C_F_POINTER(T_instance%this,this)
call C_F_PROCPOINTER(T_instance%GetRes, GetRes)

Или у меня слишком упрощенная c интерпретация вашего вопроса?

...