From 8c045d8cb47a23ec158acd15e8242f11a3c90698 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Tue, 19 Mar 2024 14:17:19 +0800 Subject: [PATCH] Expose ctol to the C interface --- c/cobyla_c.f90 | 11 ++++++++--- c/include/prima/prima.h | 8 ++++++++ c/lincoa_c.f90 | 11 ++++++++--- c/prima.c | 11 +++++++---- 4 files changed, 31 insertions(+), 10 deletions(-) diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 531b7810e7..73a8a34e13 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -13,7 +13,7 @@ module cobyla_c_mod subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ineq, Aineq, bineq, m_eq, Aeq, beq, & - & xl, xu, f0, nlconstr0, nf, rhobeg, rhoend, ftarget, maxfun, iprint, callback_ptr, info) bind(C) + & xl, xu, f0, nlconstr0, nf, rhobeg, rhoend, ftarget, maxfun, iprint, ctol, 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 :: cobyla_mod, only : cobyla @@ -47,6 +47,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(C_DOUBLE), intent(in), value :: ftarget integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: iprint +real(C_DOUBLE), intent(in), value :: ctol type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info @@ -78,6 +79,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP), allocatable :: nlconstr0_loc(:) real(RP), allocatable :: rhobeg_loc real(RP), allocatable :: rhoend_loc +real(RP), allocatable :: ctol_loc real(RP), allocatable :: xl_loc(:) real(RP), allocatable :: xu_loc(:) @@ -129,6 +131,9 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) end if iprint_loc = int(iprint, kind(iprint_loc)) +if (.not. is_nan(ctol)) then + ctol_loc = real(ctol, kind(ctol_loc)) +end if ! Call the Fortran code if (c_associated(callback_ptr)) then @@ -138,12 +143,12 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! We then provide the closure to the algorithm. call cobyla(calcfc, m_nlcon_loc, x_loc, f_loc, cstrv=cstrv_loc, nlconstr=nlconstr_loc, Aineq=Aineq_loc, & & bineq=bineq_loc, Aeq=Aeq_loc, beq=beq_loc, xl=xl_loc, xu=xu_loc, f0=f0_loc, nlconstr0=nlconstr0_loc, & - & nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, maxfun=maxfun_loc, & + & nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, ctol=ctol_loc, maxfun=maxfun_loc, & & iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) else call cobyla(calcfc, m_nlcon_loc, x_loc, f_loc, cstrv=cstrv_loc, nlconstr=nlconstr_loc, Aineq=Aineq_loc, & & bineq=bineq_loc, Aeq=Aeq_loc, beq=beq_loc, xl=xl_loc, xu=xu_loc, f0=f0_loc, nlconstr0=nlconstr0_loc, & - & nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, maxfun=maxfun_loc, & + & nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, ctol=ctol_loc, maxfun=maxfun_loc, & & iprint=iprint_loc, info=info_loc) end if diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h index e317991e6a..73042fe32a 100644 --- a/c/include/prima/prima.h +++ b/c/include/prima/prima.h @@ -220,6 +220,14 @@ typedef struct { // based on the algorithm that will be used int npt; + // ctol: tolerance for the constraint violation (COBYLA and LINCOA only) + // ctol is the tolerance of constraint violation. x is considered feasible if cstrv(x) <= ctol. + // N.B.: 1. ctol is absolute, not relative. + // 2. ctol is used for choosing the returned x. It does not affect the iterations of the algorithm. + // Default: NaN, which will be interpreted in Fortran as not present, in which case a default value + // of machine epsilon will be used. + double ctol; + // data: user data, will be passed through the objective function callback // Default: NULL void *data; diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index 088d576919..a2e9e17754 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -14,7 +14,7 @@ module lincoa_c_mod subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_eq, Aeq, beq, xl, xu, & - & nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, callback_ptr, info) bind(C) + & nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, ctol, 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 : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK @@ -45,6 +45,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: npt integer(C_INT), intent(in), value :: iprint +real(C_DOUBLE), intent(in), value :: ctol type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info @@ -72,6 +73,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP), allocatable :: rhobeg_loc real(RP), allocatable :: rhoend_loc +real(RP), allocatable :: ctol_loc real(RP), allocatable :: xl_loc(:) real(RP), allocatable :: xu_loc(:) @@ -118,6 +120,9 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ npt_loc = int(npt, kind(npt_loc)) end if iprint_loc = int(iprint, kind(iprint_loc)) +if (.not. is_nan(ctol)) then + ctol_loc = real(ctol, kind(ctol_loc)) +end if ! Call the Fortran code if (c_associated(callback_ptr)) then @@ -127,12 +132,12 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! We then provide the closure to the algorithm. call lincoa(calfun, x_loc, f_loc, cstrv=cstrv_loc, Aineq=Aineq_loc, bineq=bineq_loc, Aeq=Aeq_loc, & & beq=beq_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & + & ftarget=ftarget_loc, ctol=ctol_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & & callback_fcn=callback_fcn, info=info_loc) else call lincoa(calfun, x_loc, f_loc, cstrv=cstrv_loc, Aineq=Aineq_loc, bineq=bineq_loc, Aeq=Aeq_loc, & & beq=beq_loc, xl=xl_loc, xu=xu_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, & - & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & + & ftarget=ftarget_loc, ctol=ctol_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, & & info=info_loc) end if diff --git a/c/prima.c b/c/prima.c index 50aaf7bbd9..8cea9e7c78 100644 --- a/c/prima.c +++ b/c/prima.c @@ -57,6 +57,7 @@ int prima_init_options(prima_options_t *const options) options->rhoend = NAN; // Will be interpreted by Fortran as not present options->iprint = PRIMA_MSG_NONE; options->ftarget = -INFINITY; + options->ctol = NAN; // Will be interpreted by Fortran as not present return 0; } @@ -213,7 +214,8 @@ int cobyla_c(const int m_nlcon, const prima_objcon_t calcfc, const void *data, c const int m_eq, const double Aeq[], const double beq[], const double xl[], const double xu[], const double f0, const double nlconstr0[], - int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, const prima_callback_t callback, int *const info); + int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, const double ctol, + const prima_callback_t callback, int *const info); int bobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *const f, const double xl[], const double xu[], int *const 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 *const info); @@ -227,7 +229,8 @@ int uobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], doub int lincoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *const f, double *const cstrv, 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[], - int *const 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 *const info); + int *const nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const double ctol, + const prima_callback_t callback, int *const info); // The function that does the minimization using a PRIMA solver @@ -250,13 +253,13 @@ int prima_minimize(const prima_algorithm_t algorithm, const prima_problem_t prob case PRIMA_COBYLA: cobyla_c(problem.m_nlcon, problem.calcfc, options.data, problem.n, result->x, &(result->f), &(result->cstrv), result->nlconstr, problem.m_ineq, problem.Aineq, problem.bineq, problem.m_eq, problem.Aeq, problem.beq, - problem.xl, problem.xu, problem.f0, problem.nlconstr0, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.iprint, options.callback, &info); + problem.xl, problem.xu, problem.f0, problem.nlconstr0, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.iprint, options.ctol, options.callback, &info); break; case PRIMA_LINCOA: lincoa_c(problem.calfun, options.data, problem.n, result->x, &(result->f), &(result->cstrv), problem.m_ineq, problem.Aineq, problem.bineq, problem.m_eq, problem.Aeq, problem.beq, - problem.xl, problem.xu, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.callback, &info); + problem.xl, problem.xu, &(result->nf), options.rhobeg, options.rhoend, options.ftarget, options.maxfun, options.npt, options.iprint, options.ctol, options.callback, &info); break; case PRIMA_NEWUOA: