Amending from PR feedback

An intermediate variable was introduced so that the interface will
still work if RP=4.

The method of providing f0, rhobeg, rhoend was modified to make it
both user-friendly and still able to achieve our goal of passing
NULL to the Fortran code in the event the user does not provide a
value. Documentation was added.

c/bobyqa_c.f90

@@ -19,6 +19,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, &
 use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK
 use, non_intrinsic :: consts_mod, only : RP, IK
 use, non_intrinsic :: bobyqa_mod, only : bobyqa
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -31,8 +32,8 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, &
 type(C_PTR), intent(in), value :: xl
 type(C_PTR), intent(in), value :: xu
 integer(C_INT), intent(out) :: nf
-real(C_DOUBLE), intent(in), value :: rhobeg
-real(C_DOUBLE), intent(in), value :: rhoend
+type(C_PTR), intent(in), value :: rhobeg
+type(C_PTR), intent(in), value :: rhoend
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: npt
@@ -47,12 +48,16 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, &
 integer(IK) :: npt_loc
 integer(IK) :: nf_loc
 real(RP) :: f_loc
-real(RP) :: rhobeg_loc
-real(RP) :: rhoend_loc
+real(C_DOUBLE), pointer :: rhobeg_loc_intermediate
+real(RP), allocatable :: rhobeg_loc
+real(C_DOUBLE), pointer :: rhoend_loc_intermediate
+real(RP), allocatable :: rhoend_loc
 real(RP) :: ftarget_loc
 real(RP) :: x_loc(n)
-real(RP), pointer :: xl_loc(:)
-real(RP), pointer :: xu_loc(:)
+real(C_DOUBLE), pointer :: xl_loc_intermediate(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_intermediate(:)
+real(RP), allocatable :: xu_loc(:)
 ! The initialization to null is necessary to avoid a bug with the newer Intel compiler ifx.
 ! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013
 ! The bug was observed in all versions of ifx up to 2024.0.1. Once this bug is fixed we should remove the
@@ -63,17 +68,23 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, &
 ! Read the inputs and convert them to the Fortran side types
 x_loc = real(x, kind(x_loc))
 if (C_ASSOCIATED(xl)) then
-    call C_F_POINTER(xl, xl_loc, shape=[n])
-else
-    xl_loc => null()
+    call C_F_POINTER(xl, xl_loc_intermediate, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_intermediate, kind(xl_loc))
 end if
 if (C_ASSOCIATED(xu)) then
-    call C_F_POINTER(xu, xu_loc, shape=[n])
-else
-    xu_loc => null()
+    call C_F_POINTER(xu, xu_loc_intermediate, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_intermediate, kind(xu_loc))
+end if
+if (C_ASSOCIATED(rhobeg)) then
+    call C_F_POINTER(rhobeg, rhobeg_loc_intermediate)
+    rhobeg_loc = real(rhobeg_loc_intermediate, kind(rhobeg_loc))
+end if
+if (C_ASSOCIATED(rhoend)) then
+    call C_F_POINTER(rhoend, rhoend_loc_intermediate)
+    rhoend_loc = real(rhoend_loc_intermediate, kind(rhoend_loc))
 end if
-rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
-rhoend_loc = real(rhoend, kind(rhoend_loc))
 ftarget_loc = real(ftarget, kind(ftarget_loc))
 maxfun_loc = int(maxfun, kind(maxfun_loc))
 npt_loc = int(npt, kind(npt_loc))

c/cobyla_c.f90

@@ -18,6 +18,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 use, non_intrinsic :: cintrf_mod, only : COBJCON, CCALLBACK
 use, non_intrinsic :: consts_mod, only : RP, IK
 use, non_intrinsic :: cobyla_mod, only : cobyla
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -41,8 +42,8 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 type(C_PTR), intent(in), value :: f0
 type(C_PTR), intent(in), value :: nlconstr0
 integer(C_INT), intent(out) :: nf
-real(C_DOUBLE), intent(in), value :: rhobeg
-real(C_DOUBLE), intent(in), value :: rhoend
+type(C_PTR), intent(in), value :: rhobeg
+type(C_PTR), intent(in), value :: rhoend
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: iprint
@@ -62,14 +63,20 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 real(RP) :: cstrv_loc
 real(RP) :: nlconstr_loc(m_nlcon)
 real(RP) :: f_loc
-real(RP) :: rhobeg_loc
-real(RP) :: rhoend_loc
+real(C_DOUBLE), pointer :: rhobeg_loc_intermediate
+real(RP), allocatable :: rhobeg_loc
+real(C_DOUBLE), pointer :: rhoend_loc_intermediate
+real(RP), allocatable :: rhoend_loc
 real(RP) :: ftarget_loc
 real(RP) :: x_loc(n)
-real(RP), pointer :: xl_loc(:)
-real(RP), pointer :: xu_loc(:)
-real(RP), pointer :: f0_loc
-real(RP), pointer :: nlconstr0_loc(:)
+real(C_DOUBLE), pointer :: xl_loc_intermediate(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_intermediate(:)
+real(RP), allocatable :: xu_loc(:)
+real(C_DOUBLE), pointer :: f0_loc_intermediate
+real(RP), allocatable :: f0_loc
+real(C_DOUBLE), pointer :: nlconstr0_loc_intermediate(:)
+real(RP), allocatable :: nlconstr0_loc(:)
 ! The initialization to null is necessary to avoid a bug with the newer Intel compiler ifx.
 ! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013
 ! The bug was observed in all versions of ifx up to 2024.0.1. Once this bug is fixed we should remove the
@@ -86,27 +93,32 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 Aeq_loc = real(transpose(Aeq), kind(Aeq_loc))
 beq_loc = real(beq, kind(beq_loc))
 if (C_ASSOCIATED(xl)) then
-    call C_F_POINTER(xl, xl_loc, shape=[n])
-else
-    xl_loc => null()
+    call C_F_POINTER(xl, xl_loc_intermediate, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_intermediate, kind(xl_loc))
 end if
 if (C_ASSOCIATED(xu)) then
-    call C_F_POINTER(xu, xu_loc, shape=[n])
-else
-    xu_loc => null()
+    call C_F_POINTER(xu, xu_loc_intermediate, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_intermediate, kind(xu_loc))
 end if
 if (C_ASSOCIATED(f0)) then
-    call C_F_POINTER(f0, f0_loc)
-else
-    f0_loc => null()
+    call C_F_POINTER(f0, f0_loc_intermediate)
+    f0_loc = real(f0_loc_intermediate, kind(f0_loc))
 end if
 if (C_ASSOCIATED(nlconstr0)) then
-    call C_F_POINTER(nlconstr0, nlconstr0_loc, shape=[m_nlcon])
-else
-    nlconstr0_loc => null()
+    call C_F_POINTER(nlconstr0, nlconstr0_loc_intermediate, shape=[m_nlcon])
+    call safealloc(nlconstr0_loc, int(m_nlcon, IK))
+    nlconstr0_loc = real(nlconstr0_loc_intermediate, kind(nlconstr0_loc))
+end if
+if (C_ASSOCIATED(rhobeg)) then
+    call C_F_POINTER(rhobeg, rhobeg_loc_intermediate)
+    rhobeg_loc = real(rhobeg_loc_intermediate, kind(rhobeg_loc))
+end if
+if (C_ASSOCIATED(rhoend)) then
+    call C_F_POINTER(rhoend, rhoend_loc_intermediate)
+    rhoend_loc = real(rhoend_loc_intermediate, kind(rhoend_loc))
 end if
-rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
-rhoend_loc = real(rhoend, kind(rhoend_loc))
 ftarget_loc = real(ftarget, kind(ftarget_loc))
 maxfun_loc = int(maxfun, kind(maxfun_loc))
 iprint_loc = int(iprint, kind(iprint_loc))

c/examples/cobyla/cobyla_example.c

@@ -6,6 +6,7 @@
 #include <stdint.h>
 
 #define M_NLCON 1
+#define PROVIDE_INITIAL_F_AND_NLCONSTR 0
 
 static void fun(const double x[], double *f, double constr[], const void *data)
 {
@@ -52,6 +53,11 @@ int main(int argc, char * argv[])
   double xu[2] = {6.0, 6.0};
   problem.xl = xl;
   problem.xu = xu;
+#if PROVIDE_INITIAL_F_AND_NLCONSTR
+  double nlconstr0[M_NLCON] = {};
+  fun(x0, &(problem.f0), nlconstr0, NULL);
+  problem.nlconstr0 = nlconstr0;
+#endif
   // set up the options
   prima_options_t options;
   prima_init_options(&options);

c/include/prima/prima.h

@@ -129,6 +129,9 @@ typedef struct {
   // callback function to report algorithm progress (default=NULL)
   prima_callback_t callback;
 
+  double *_rhobeg; // internal - do not use
+  double *_rhoend; // internal - do not use
+
 } prima_options_t;
 
 /* Initialize problem */
@@ -171,8 +174,9 @@ typedef struct {
   int m_nlcon;
 
   // should be set to the objective function value and constraints values of the starting X, cobyla-only
-  double *f0;
+  double f0;
   double *nlconstr0;
+  double *_f0; // internal, do not use
 
 } prima_problem_t;
 

c/lincoa_c.f90

@@ -19,6 +19,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
 use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK
 use, non_intrinsic :: consts_mod, only : RP, IK
 use, non_intrinsic :: lincoa_mod, only : lincoa
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -38,8 +39,8 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
 type(C_PTR), intent(in), value :: xl
 type(C_PTR), intent(in), value :: xu
 integer(C_INT), intent(out) :: nf
-real(C_DOUBLE), intent(in), value :: rhobeg
-real(C_DOUBLE), intent(in), value :: rhoend
+type(C_PTR), intent(in), value :: rhobeg
+type(C_PTR), intent(in), value :: rhoend
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: npt
@@ -59,12 +60,16 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
 real(RP) :: beq_loc(m_eq)
 real(RP) :: cstrv_loc
 real(RP) :: f_loc
-real(RP) :: rhobeg_loc
-real(RP) :: rhoend_loc
+real(C_DOUBLE), pointer :: rhobeg_loc_intermediate
+real(RP), allocatable :: rhobeg_loc
+real(C_DOUBLE), pointer :: rhoend_loc_intermediate
+real(RP), allocatable :: rhoend_loc
 real(RP) :: ftarget_loc
 real(RP) :: x_loc(n)
-real(RP), pointer :: xl_loc(:)
-real(RP), pointer :: xu_loc(:)
+real(C_DOUBLE), pointer :: xl_loc_intermediate(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_intermediate(:)
+real(RP), allocatable :: xu_loc(:)
 ! The initialization to null is necessary to avoid a bug with the newer Intel compiler ifx.
 ! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013
 ! The bug was observed in all versions of ifx up to 2024.0.1. Once this bug is fixed we should remove the
@@ -81,17 +86,23 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
 Aeq_loc = real(transpose(Aeq), kind(Aeq_loc))
 beq_loc = real(beq, kind(beq_loc))
 if (C_ASSOCIATED(xl)) then
-    call C_F_POINTER(xl, xl_loc, shape=[n])
-else
-    xl_loc => null()
+    call C_F_POINTER(xl, xl_loc_intermediate, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_intermediate, kind(xl_loc))
 end if
 if (C_ASSOCIATED(xu)) then
-    call C_F_POINTER(xu, xu_loc, shape=[n])
-else
-    xu_loc => null()
+    call C_F_POINTER(xu, xu_loc_intermediate, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_intermediate, kind(xu_loc))
+end if
+if (C_ASSOCIATED(rhobeg)) then
+    call C_F_POINTER(rhobeg, rhobeg_loc_intermediate)
+    rhobeg_loc = real(rhobeg_loc_intermediate, kind(rhobeg_loc))
+end if
+if (C_ASSOCIATED(rhoend)) then
+    call C_F_POINTER(rhoend, rhoend_loc_intermediate)
+    rhoend_loc = real(rhoend_loc_intermediate, kind(rhoend_loc))
 end if
-rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
-rhoend_loc = real(rhoend, kind(rhoend_loc))
 ftarget_loc = real(ftarget, kind(ftarget_loc))
 maxfun_loc = int(maxfun, kind(maxfun_loc))
 npt_loc = int(npt, kind(npt_loc))

c/newuoa_c.f90

@@ -13,7 +13,7 @@ module newuoa_c_mod
 
 
 subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, callback_ptr, info) bind(C)
-use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED, C_F_PROCPOINTER
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED, C_F_PROCPOINTER, C_F_POINTER
 use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK
 use, non_intrinsic :: consts_mod, only : RP, IK
 use, non_intrinsic :: newuoa_mod, only : newuoa
@@ -27,8 +27,8 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 real(C_DOUBLE), intent(inout) :: x(n)
 real(C_DOUBLE), intent(out) :: f
 integer(C_INT), intent(out) :: nf
-real(C_DOUBLE), intent(in), value :: rhobeg
-real(C_DOUBLE), intent(in), value :: rhoend
+type(C_PTR), intent(in), value :: rhobeg
+type(C_PTR), intent(in), value :: rhoend
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: npt
@@ -43,8 +43,10 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 integer(IK) :: npt_loc
 integer(IK) :: nf_loc
 real(RP) :: f_loc
-real(RP) :: rhobeg_loc
-real(RP) :: rhoend_loc
+real(C_DOUBLE), pointer :: rhobeg_loc_intermediate
+real(RP), allocatable :: rhobeg_loc
+real(C_DOUBLE), pointer :: rhoend_loc_intermediate
+real(RP), allocatable :: rhoend_loc
 real(RP) :: ftarget_loc
 real(RP) :: x_loc(n)
 ! The initialization to null is necessary to avoid a bug with the newer Intel compiler ifx.
@@ -56,8 +58,14 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 
 ! Read the inputs and convert them to the Fortran side types
 x_loc = real(x, kind(x_loc))
-rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
-rhoend_loc = real(rhoend, kind(rhoend_loc))
+if (C_ASSOCIATED(rhobeg)) then
+    call C_F_POINTER(rhobeg, rhobeg_loc_intermediate)
+    rhobeg_loc = real(rhobeg_loc_intermediate, kind(rhobeg_loc))
+end if
+if (C_ASSOCIATED(rhoend)) then
+    call C_F_POINTER(rhoend, rhoend_loc_intermediate)
+    rhoend_loc = real(rhoend_loc_intermediate, kind(rhoend_loc))
+end if
 ftarget_loc = real(ftarget, kind(ftarget_loc))
 maxfun_loc = int(maxfun, kind(maxfun_loc))
 npt_loc = int(npt, kind(npt_loc))