Passing external parameters to functions

[giannilmbd]

Hi,
Thanks for making these codes available.
I'm trying to find a way to pass additional parameters (beisdes the argument) to the objective function, eg instead of f(x) have f(x, parameters), where "parameters" could be an ad-hoc structure with approriate number of elements (eg parameters%alpha which is set in the main program).
Is there a quick way to do so in you root_module?
Thanks
Gianni

[giannilmbd]

Just to be clearer, among the parameters there might be data which is read from files in the main calling program

[ivan-pi]

Hi @giannilmbd,

You can do this by extending Jacob's modules and "attaching" a new interface:

module rootx_module

   use root_module, only: root_scalar_original => root_scalar, func2, &
     wp => root_module_rk
   implicit none
   private

   public :: paramfunc, root_scalar, wp

   abstract interface
      function paramfunc(x,params) result(f)
         import wp
         implicit none
         real(wp), intent(in) :: x  ! independent variable
         class(*), intent(inout), optional :: params
         real(wp) :: f
      end function
   end interface

   interface root_scalar
      module procedure root_scalar_original, &
                       root_scalar_by_name_param
   end interface

contains

   subroutine root_scalar_by_name_param(method,fun,ax,bx,xzero,fzero,iflag,&
                                        ftol,rtol,atol,maxiter,fax,fbx,&
                                        bisect_on_failure, &
                                        params)
   implicit none

   character(len=*),intent(in)   :: method   !! the method to use
   procedure(paramfunc)          :: fun      !! user function to find the root of
   real(wp),intent(in)           :: ax       !! left endpoint of initial interval
   real(wp),intent(in)           :: bx       !! right endpoint of initial interval
   real(wp),intent(out)          :: xzero    !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
   real(wp),intent(out)          :: fzero    !! value of `f` at the root (`f(xzero)`)
   integer,intent(out)           :: iflag    !! status flag (`-1`=error, `0`=root found, `-999`=invalid method)
   real(wp),intent(in),optional  :: ftol     !! absolute tolerance for `f=0`
   real(wp),intent(in),optional  :: rtol     !! relative tol for x
   real(wp),intent(in),optional  :: atol     !! absolute tol for x
   integer,intent(in),optional   :: maxiter  !! maximum number of iterations
   real(wp),intent(in),optional  :: fax      !! if `f(ax)` is already known, it can be input here
   real(wp),intent(in),optional  :: fbx      !! if `f(bx)` is already known, it can be input here
   logical,intent(in),optional   :: bisect_on_failure  !! if true, then if the specified method fails,
                                                       !! it will be retried using the bisection method.
                                                       !! (default is False). Note that this can use up
                                                       !! to `maxiter` additional function evaluations.
   class(*), intent(inout) :: params   !! user parameters

   call root_scalar_original(r,fun,ax,bx,xzero,fzero,iflag,&
                             ftol,rtol,atol,maxiter,fax,fbx,&
                             bisect_on_failure)

   contains

      ! pass user parameters using host association
      function fun_adapter(x) result(f)
         real(wp), intent(in) :: x
         real(wp) :: f
         f = fun(x,params)
      end function

   end subroutine root_scalar_by_name_param

end module

Assuming I didn't make any errors, you could now use the module as follows:

program main
use rootx_module, only: root_scalar, wp
implicit none

type :: my_params
  real(wp) :: a, b
end type 
type(my_params) :: p

  real(wp) :: x, f
  integer :: iflag, unit
  
  open(newunit=unit,file="param.in")
  read(unit,*) p%a, p%b
  close(unit)

  call root_scalar('bisection',func,-9.0_wp,31.0_wp,x,f,iflag, &
    params=p)
  ! remember to use the keyword form `params=...` for clarity!

  write(*,*) 'f(',x,') = ', f
  write(*,*) 'iflag    = ', iflag

contains

  function func(x,params) result(f)

    real(wp), intent(in) :: x
    class(*), intent(inout), optional :: params
    real(wp) :: f

    select type(params)
    type is (my_params)
    f = my_params%a * x * exp(my_params%b*x)
    end select type
  end function func

end program

[ivan-pi]

I wonder if it's possible to write a cpp/fpp or fypp macro which would cover most of the wrapping for a specific procedure callback abstract interface.

Aside from the callback adapter shown above, the other options that are always present are host-association and module variables (can be tough to make this thread-safe).

@jacobwilliams, would an example or tutorial fit into the scope of this repository?