Modern Fortran interface for NLopt

Help
15

This weekend I’ve been playing around with a Fortran wrapper of CppNumericalSolvers, a C++ header-only library of numerical optimization methods with an interface striving to be close to the MATLAB’s fminsearch.

I haven’t published the code yet, but wanted to give an example of what the Fortran interface currently looks like. The following callbacks define the problem:

module rosenbrock_problem

   use cppoptlib, only: wp => cppoptlib_wp

   implicit none
   public

contains

   !> Objective F(x)
   function rosenbrock(x,params) result(y)
      real(wp), intent(in) :: x(:)
      class(*), intent(in), optional :: params    ! not used
      real(wp) :: y, t1, t2
      t1 = 1 - x(1)
      t2 = x(2) - x(1)**2
      y = t1**2 + 100*t2*t2
   end function

   !> Objective gradient, dF(x)/dx
   subroutine rosenbrock_grad(x,grad,params)
      real(wp), intent(in) :: x(:)
      real(wp), intent(out) :: grad(:)
      class(*), intent(in), optional :: params    ! not used
      grad(1) = -2*(1-x(1))+200*(x(2)-x(1)**2)*(-2*x(1))
      grad(2) = 200*(x(2)-x(1)**2)
   end subroutine

end module

The driver turned out quite pleasant:

program minimize

    use cppoptlib, only: wp => cppoptlib_wp, problem_type, solver
    use rosenbrock_problem
    implicit none

    type(problem_type) :: problem
    real(wp) :: x(2), res

    call problem%set(rosenbrock,rosenbrock_grad) ! or
    ! problem = problem_type(rosenbrock,rosenbrock_grad)
    
    x = [-1._wp, 2._wp] ! initial guess

    print *, "initial guess:", x

    call problem%solve(res,x,solver=solver%bfgs)

    print *, "objective: ", res
    print *, "solution: ", x

end program

In contrast to NLopt which stores a pointer to an opaque C struct that stores the “state” of the numerical solver, in the present library the problem_type only stores pointers to the callback functions. Any C++ structures only live for the duration of the call to type-bound solve method.

It is also possible to use internal procedures instead of module procedures as the callback arguments:

program simple_problem

   use cppoptlib, only: wp => cppoptlib_wp, problem_type
   implicit none

   real(wp) :: x(2), res
   type(problem_type) :: problem

   x = [-1, 2]

   problem = problem_type(objective,gradient)

   call problem%solve(res,x)

   write(*,'(A,F0.6)') "f in argmin ", res
   write(*,'(A)') "solver status "
   write(*,'(A)') "solver criteria "

contains

   real(wp) function objective(x,params)
      real(wp), intent(in) :: x(:)
      class(*), intent(in), optional :: params  ! not used
      objective = 5*x(1)**2 + 100*x(1)**2 + 5
   end function

   subroutine gradient(x,grad,params)
      real(wp), intent(in) :: x(:)
      real(wp), intent(out) :: grad(:)
      class(*), intent(in), optional :: params  ! not used
      grad(1) = 2*5*x(1)
      grad(2) = 2*100*x(2)
   end subroutine

end program

I am currently investigating how to pass custom stopping criteria and printing callbacks. So far I quite like the interface. Having a “stateless” derived type, which only serves as a collection of callbacks seems to be a useful architectural pattern.

Say I wanted to solve the same parameterized problem multiple times in parallel (an embarrassingly parallel problem), presumably, I could achieve this easily with:

program parallel_minimization
   use cppoptlib, only: wp => cppoptlib_wp, problem_type
   implicit none

   type :: problem_params
      real(wp) :: a, b, c
   end type

   integer, parameter :: nproblems = 100
   type(problem_type) :: problem
   type(problem_params), allocatable :: params(:)
   real(wp), allocatable :: obj(:), x(:,:)
   integer :: i

   problem = problem_type(objective,gradient)

   allocate(params(nproblems),obj(nproblems),x(2,nproblems))
   ! [insert] ... initialize params from file or namelist

! we use a dynamic schedule to parallelize the loop because each solve might 
! take a different number of function evaluations. we can fine-tune this 
! with the chunk-size. to reduce over-head from scheduling, a guided schedule 
! could also be used.

   !$omp parallel do schedule(dynamic) default(private) shared(problem)
   do i = 1, nproblems
      call problem%solve(res(i),x(:,i),params(i))
   end do
   !$omp end parallel do

contains
   ! [insert] ... callbacks
end program

Does anyone know if sharing procedure pointers within a parallel section is a “legal” use? Are there any problems which could pop up? (Assume that the callbacks and type-bound solve routine reference no global data, i.e. they are thread-safe.)

Edit: a second embarassingly parallel usage case would be to start the solver from multiple initial guesses (maybe according to a Latin hypercube sampling). Picking the best solution (presumably the global minimum) is as simple as:

xbest = x(:,minloc(res,dim=1))
4 Likes