fortran-lang / minpack

Modernized Minpack: for solving nonlinear equations and nonlinear least squares problems
https://fortran-lang.github.io/minpack/
Other
94 stars 20 forks source link

Pass data as pointer #4

Closed Nicholaswogan closed 2 years ago

Nicholaswogan commented 2 years ago

It would be nice if data could be passed to the objective function. This could be done with a void C pointer, type(c_ptr), or any other method. Right now data must be passed globally, which makes programs not thread-safe.

jacobwilliams commented 2 years ago

You don't have to use global variables. You can wrap the call in a class method. So your data lives in the class. Something like this:

module my_module

  use minpack_module
  use iso_fortran_env, only: wp => real64

  implicit none

  type :: my_solver_t
    integer :: some_data ! your data
    contains
    procedure :: solve
  end type my_solver_t

  type(my_solver_t) :: my_solver_1
  type(my_solver_t) :: my_solver_2

  ...

  contains

  subroutine solve(me)

    class(my_solver_t),intent(inout) :: me

    call hybrd(fcn, ...)

  contains

    subroutine fcn(n,x,fvec,iflag)

    integer,intent(in) :: n  
    real(wp),intent(in) :: x(n)  
    real(wp),intent(out) :: fvec(n)  
    integer,intent(inout) :: iflag 

    ! can use me%some_data here !
    ...

    end subroutine fcn

  end subroutine solve

end module my_module

Eventually I might just make this part of the library. This is the modern way to do it.

Nicholaswogan commented 2 years ago

Thanks!!

ivan-pi commented 2 years ago

The pattern shown by @jacobwilliams is quite interesting. On the other hand it's a kind of Frankenstein between data container, callback function, and the actual solver algorithm. Say I want to swap between MINPACK and Ceres (or any other least squares solver). Where is my design breakpoint?

Would I introduce separate procedures, i.e.

subroutine solve_minpack()
! ...
subroutine solve_ceres()
! ...

Or would I introduce a method argument? Here's some pseudo-code

subroutine solve(me,method)
type(my_solver_t), intent(inout) :: me
character(len=*), intent(in), optional :: method
character(len=:), allocatable :: method_

method_ = "minpack"
if (present(method)) method_ = lowercase(method)

select case (method_)
case('minpack')
   call hybrd(...)
case('ceres')
   call ceres_optimize(...) ! ...
case default
  error stop 'my_solver_t%solve: unknown method, available options are "minpack" (default) or "ceres"'
end select

end subroutine