Amending from PR feedback

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

The method of providing f0 was changed to be a little bit more user
friendly, and rhobeg and rhoend adopted this approach.

Additionally, the "unset" value of maxfun and npt was changed to 0
based on discussion in the PR.

.github/actions/spelling/allow.txt

@@ -2143,3 +2143,4 @@ nwrite
 stardev
 noutput
 nprobinfo
+interm

.github/workflows/cmake.yml

@@ -160,7 +160,7 @@ jobs:
       matrix:
         toolchain:
           - {compiler: aflang, cflags: '-Wall', fflags: '-Wpedantic -Weverything -Wall -Wextra -Minform=warn -Mstandard -Mrecursive'}
-          - {compiler: nvfortran, cflags: '-Wall', fflags: '-C -Wall -Wextra -Minform=warn -Mstandard -Mrecursive -Mbounds -Mchkstk -Mchkptr'}
+          - {compiler: nvfortran, cflags: '-Wall', fflags: '-C -Wall -Wextra -Minform=warn -Mstandard -Mrecursive -Mbounds -Mchkstk'}
           - {compiler: flang, cflags: '-Wall', fflags: '-Wpedantic -Weverything -Wall -Wextra'}
 
     steps:

c/README.txt

@@ -2,4 +2,18 @@ This is the C interface for calling the Fortran version of PRIMA.
 
 Contributed by Julien Schueller ( schueller@phimeca.com ) in September 2023.
 
+## Usage
 See https://github.com/libprima/prima#c for the usage.
+
+## Developer notes
+Since C does not provide a direct avenue for optional function arguments, we use the following values to indicate that a value is not present:
+
+| type | value |
+| ---- | ----- |
+| integer | 0 |
+| double  | NAN |
+| array   | NULL |
+
+The rationale for using 0 to indicate that an integer argument is not present lies in the fact that most of the time the code expects any integer arguments to be positive integers. See https://github.com/libprima/prima/pull/135 for a more detailed discussion.
+
+For doubles, we set them to NAN but the check for NAN happens in the Fortran code using the `is_nan` function from the `infnan` module. See `rhobeg` and `rhoend` for an example.

c/bobyqa_c.f90

@@ -16,9 +16,11 @@ module bobyqa_c_mod
 subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, 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, C_F_POINTER
+use, non_intrinsic :: bobyqa_mod, only : bobyqa
 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 :: infnan_mod, only : is_nan
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -47,12 +49,14 @@ 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(RP), allocatable :: rhobeg_loc
+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_interm(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_interm(:)
+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 +67,21 @@ 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_interm, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_interm, 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_interm, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_interm, kind(xu_loc))
+end if
+if (.not. is_nan(rhobeg)) then
+    rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
+end if
+if (.not. is_nan(rhoend)) then
+    rhoend_loc = real(rhoend, 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

@@ -16,8 +16,10 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
     & xl, xu, f0, nlconstr0, nf, rhobeg, rhoend, ftarget, maxfun, 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, C_F_POINTER
 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 :: consts_mod, only : RP, IK
+use, non_intrinsic :: infnan_mod, only : is_nan
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -38,7 +40,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
 real(C_DOUBLE), intent(in) :: beq(m_eq)
 type(C_PTR), intent(in), value :: xl
 type(C_PTR), intent(in), value :: xu
-type(C_PTR), intent(in), value :: f0
+real(C_DOUBLE), intent(in), value :: f0
 type(C_PTR), intent(in), value :: nlconstr0
 integer(C_INT), intent(out) :: nf
 real(C_DOUBLE), intent(in), value :: rhobeg
@@ -62,14 +64,17 @@ 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(RP), allocatable :: rhobeg_loc
+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_interm(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_interm(:)
+real(RP), allocatable :: xu_loc(:)
+real(RP), allocatable :: f0_loc
+real(C_DOUBLE), pointer :: nlconstr0_loc_interm(:)
+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 +91,29 @@ 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_interm, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_interm, 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_interm, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_interm, kind(xu_loc))
 end if
-if (C_ASSOCIATED(f0)) then
-    call C_F_POINTER(f0, f0_loc)
-else
-    f0_loc => null()
+if (.not. is_nan(f0)) then
+    f0_loc = real(f0, 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_interm, shape=[m_nlcon])
+    call safealloc(nlconstr0_loc, int(m_nlcon, IK))
+    nlconstr0_loc = real(nlconstr0_loc_interm, kind(nlconstr0_loc))
+end if
+if (.not. is_nan(rhobeg)) then
+    rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
+end if
+if (.not. is_nan(rhoend)) then
+    rhoend_loc = real(rhoend, 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)
 {
@@ -40,6 +41,11 @@ int main(int argc, char * argv[])
   problem.calcfc = &fun;
   problem.x0 = x0;
   problem.m_nlcon = M_NLCON;
+#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

@@ -104,13 +104,13 @@ typedef void (*prima_callback_t)(const int n, const double x[], const double f,
 
 typedef struct {
 
-  // a reasonable initial change to the variables (default=1.0)
+  // a reasonable initial change to the variables
   double rhobeg;
 
-  // required accuracy for the variables (default=1e-6)
+  // required accuracy for the variables
   double rhoend;
 
-  // maximum number of function evaluations (default=-1 interpreted as 500*n)
+  // maximum number of function evaluations (default=500*n)
   int maxfun;
 
   // verbosity level, see the prima_message_t enum (default=PRIMA_MSG_NONE)
@@ -119,7 +119,7 @@ typedef struct {
   // target function value; optimization stops when f <= ftarget for a feasible point (default=-inf)
   double ftarget;
 
-  // number of points in the interpolation set n+2<=npt<=(n+1)(n+2)/2 (default=-1 interpreted as 2*n+1)
+  // number of points in the interpolation set n+2<=npt<=(n+1)(n+2)/2 (default=2*n+1)
   // ignored for uobyqa & cobyla
   int npt;
 
@@ -171,7 +171,7 @@ 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;
 
 } prima_problem_t;

c/lincoa_c.f90

@@ -18,7 +18,9 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
 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 :: infnan_mod, only : is_nan
 use, non_intrinsic :: lincoa_mod, only : lincoa
+use, non_intrinsic :: memory_mod, only : safealloc
 implicit none
 
 ! Compulsory arguments
@@ -59,12 +61,14 @@ 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(RP), allocatable :: rhobeg_loc
+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_interm(:)
+real(RP), allocatable :: xl_loc(:)
+real(C_DOUBLE), pointer :: xu_loc_interm(:)
+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 +85,21 @@ 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_interm, shape=[n])
+    call safealloc(xl_loc, int(n, IK))
+    xl_loc = real(xl_loc_interm, 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_interm, shape=[n])
+    call safealloc(xu_loc, int(n, IK))
+    xu_loc = real(xu_loc_interm, kind(xu_loc))
+end if
+if (.not. is_nan(rhobeg)) then
+    rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
+end if
+if (.not. is_nan(rhoend)) then
+    rhoend_loc = real(rhoend, 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

@@ -16,6 +16,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED, C_F_PROCPOINTER
 use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK
 use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: infnan_mod, only : is_nan
 use, non_intrinsic :: newuoa_mod, only : newuoa
 implicit none
 
@@ -43,8 +44,8 @@ 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(RP), allocatable :: rhobeg_loc
+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 +57,12 @@ 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 (.not. is_nan(rhobeg)) then
+    rhobeg_loc = real(rhobeg, kind(rhobeg_loc))
+end if
+if (.not. is_nan(rhoend)) then
+    rhoend_loc = real(rhoend, 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))

c/prima.c

@@ -15,8 +15,8 @@ int prima_init_options(prima_options_t *options)
   {
     memset(options, 0, sizeof(prima_options_t));
     options->maxfun = 0;  // interpreted as maxfun taking the default value MAXFUN_DIM_DFT*n
-    options->rhobeg = 1.0;
-    options->rhoend = 1e-6;
+    options->rhobeg = NAN;  // interpreted by Fortran as not present
+    options->rhoend = NAN;  // interpreted by Fortran as not present
     options->iprint = PRIMA_MSG_NONE;
     options->ftarget = -INFINITY;
     options->npt = 0;  // interpreted as npt taking the default value 2*n+1
@@ -32,6 +32,7 @@ int prima_init_problem(prima_problem_t *problem, int n)
   {
     memset(problem, 0, sizeof(prima_problem_t));
     problem->n = n;
+    problem->f0 = NAN;
     return 0;
   }
   else
@@ -43,7 +44,7 @@ int cobyla_c(const int m_nlcon, const prima_objcon_t calcfc, const void *data, c
              const int m_ineq, const double Aineq[], const double bineq[],
              const int m_eq, const double Aeq[], const double beq[],
              const double xl[], const double xu[],
-             double *f0, const double nlconstr0[],
+             const double f0, const double nlconstr0[],
              int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, const prima_callback_t callback, int *info);
 int bobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f, const double xl[], const double xu[],
              int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *info);
@@ -79,9 +80,9 @@ int prima_check_problem(prima_problem_t *problem, prima_options_t *options, cons
     return PRIMA_NULL_X0;
   if ((use_constr && !problem->calcfc) || (!use_constr && !problem->calfun))
     return PRIMA_NULL_FUNCTION;
-  if (options->maxfun < 0)
+  if (options->maxfun == 0)
     options->maxfun = MAXFUN_DIM_DFT*problem->n;
-  if (options->npt < 0)
+  if (options->npt == 0)
     options->npt = 2*problem->n+1;
   return 0;
 }