gradient has the wrong dimension?

[Nicholaswogan]

Looks like the callback function nlopt_mfunc might be slightly wrong? I think the gradient variable should be m*n in length.

https://github.com/grimme-lab/nlopt-f/blob/7e26398dee7bbef656f5fdc626911925a785f12d/src/nlopt_callback.f90#L152

Relevant part of the docs:

https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#vector-valued-constraints

[awvwgk]

Good catch. Would you mind sending a patch for this?

[Nicholaswogan]

Ya, I'll send a patch, but it won't be for a week or two.

[ivan-pi]

I believe this is the relevant section of the NLopt documentation:

In addition, if grad is non-NULL, then grad points to an array of length m*n which should, upon return, be set to the gradients of the constraint functions with respect to x. The n dimension of grad is stored contiguously, so that $\partial c_i/\partial x_j$ is stored in grad[i*n + j].

Would either you be in favor of passing the vector-constraint gradient as a rank-2 array on the Fortran side of the API? In other words, the mfunc callback interface should be:

    subroutine nlopt_mfunc_interface(result, x, gradient, func_data)
      import :: c_int, c_double, c_ptr
      implicit none
      real(c_double), intent(inout) :: result(:)
      real(c_double), intent(in) :: x(:)
      real(c_double), intent(inout), optional :: gradient(:,:)  ! rank-2 instead of rank-1
      class(*), intent(in), optional :: func_data
    end subroutine

In the C-conforming callback wrapper, you would then use pointer remapping before passing it to the Fortran callback

  subroutine wrap_mfunc(m, result, n, x, gradient, func_data) &
      bind(c, name="nlopt_wrap_mfunc")
    integer(c_int), value :: m
    real(c_double), intent(inout) :: result(m)
    integer(c_int), value :: n
    real(c_double), intent(in) :: x(n)
    real(c_double), intent(inout), optional, target :: gradient(n*m)
    type(c_ptr), value :: func_data

    type(nlopt_mfunc), pointer :: mfunc
    real(c_double), pointer :: gradient_(:,:) => null()

    call c_f_pointer(func_data, mfunc)

    ! pointer remapping from rank-1 to rank-2 array expected by Fortran callback interface
    gradient_(1:n,1:m) => gradient

    call mfunc%f(result, x, gradient_, mfunc%f_data)
  end subroutine

I recently used the exact same approach in the libdogleg-f interface (https://github.com/ivan-pi/libdogleg-f/blob/32b65ae4f880a052768c05d24c24b79920710681/src/dogleg_callback.f90#L94), (edit) with the minor difference, that gradient was a required variable.

Caveat, I'm not sure what happens in case gradient is not present... Perhaps a check if gradient is present is needed, as follows:

if (present(gradient)) then
   gradient_(1:n,1:m) => gradient
else
   gradient_ => null()
end if

[awvwgk]

I don't think you need bound remapping, just declare the input array as rank 2 directly:

  subroutine wrap_mfunc(m, result, n, x, gradient, func_data) &
      bind(c, name="nlopt_wrap_mfunc")
    integer(c_int), value :: m
    real(c_double), intent(inout) :: result(m)
    integer(c_int), value :: n
    real(c_double), intent(in) :: x(n)
    real(c_double), intent(inout), optional :: gradient(n, m)
    type(c_ptr), value :: func_data

    type(nlopt_mfunc), pointer :: mfunc

    call c_f_pointer(func_data, mfunc)

    call mfunc%f(result, x, gradient, mfunc%f_data)
  end subroutine

It is only the Fortran side which cares about the shape of an explicit shape array, the C side will pass a pointer to a contiguous slice of n*m doubles regardless.

[ivan-pi]

Right, even better 👍🏻 .

Luckily, the dimensions are part of the prototype. In case of libdogleg, I had to deal with assumed-size arrays, and hence needed pointer remapping or I would get a rank mismatch:

   subroutine dl_callback_adaptor(p, x, J, cookie) bind(c)
      real(c_double), intent(in), target :: p(*)
         !> p is the current state vector
      real(c_double), intent(inout) :: x(*)
      real(c_double), intent(inout), target :: J(*)
      type(c_ptr), value :: cookie

      type(dl_dense_t), pointer :: cb 
      real(c_double), pointer :: Jf(:,:) => null()
      call c_f_pointer(cookie, cb)
      associate(nstate => cb%nstate, nmeas => cb%nmeas)
         Jf(1:nstate,1:nmeas) => J(1:nstate*nmeas)
         call cb%f(p(1:nstate), x(1:nmeas), Jf, cb%fparams)
      end associate
   end subroutine