From b11ac4338d8070bb22bcb488720039e75389a407 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Thu, 23 Nov 2023 18:10:48 +0800 Subject: [PATCH 01/13] Manually applied jschueller's closed callback PR --- .github/actions/spelling/allow.txt | 4 +- c/bobyqa_c.f90 | 83 +++++++++++++++++++++++--- c/cintrf.f90 | 75 +++++------------------ c/cobyla_c.f90 | 96 ++++++++++++++++++++++++++---- c/examples/bobyqa/bobyqa_example.c | 10 ++++ c/examples/cobyla/cobyla_example.c | 12 ++++ c/examples/lincoa/lincoa_example.c | 10 ++++ c/examples/newuoa/newuoa_example.c | 10 ++++ c/examples/uobyqa/uobyqa_example.c | 10 ++++ c/include/prima/prima.h | 23 ++++++- c/lincoa_c.f90 | 94 +++++++++++++++++++++++++---- c/newuoa_c.f90 | 81 ++++++++++++++++++++++--- c/prima.c | 22 +++---- c/uobyqa_c.f90 | 83 +++++++++++++++++++++++--- fortran/bobyqa/bobyqa.f90 | 20 +++++-- fortran/bobyqa/bobyqb.f90 | 33 +++++++++- fortran/cobyla/cobyla.f90 | 17 ++++-- fortran/cobyla/cobylb.f90 | 28 +++++++-- fortran/common/infos.f90 | 3 +- fortran/common/message.f90 | 4 +- fortran/common/pintrf.f90 | 15 ++++- fortran/lincoa/lincoa.f90 | 23 +++++-- fortran/lincoa/lincob.f90 | 31 +++++++++- fortran/newuoa/newuoa.f90 | 17 ++++-- fortran/newuoa/newuob.f90 | 31 +++++++++- fortran/uobyqa/uobyqa.f90 | 17 ++++-- fortran/uobyqa/uobyqb.f90 | 31 +++++++++- 27 files changed, 726 insertions(+), 157 deletions(-) diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt index c6b95fe345..dfb561edf4 100644 --- a/.github/actions/spelling/allow.txt +++ b/.github/actions/spelling/allow.txt @@ -45,6 +45,7 @@ BQPGASIM BROWNAL Bindel Broyden +CCALLBACK CGRAD CHANDHEQ CHCKTST @@ -2048,8 +2049,6 @@ COBJCON cobjfun cobjfuncon constrc -evalcobj -evalcobjcon execstack FUNPTR PROCPOINTER @@ -2136,3 +2135,4 @@ orthtol nouninit libgfortran chocolatey +fcn diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 7cb90a4466..5dd07cf520 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -13,9 +13,10 @@ module bobyqa_c_mod contains -subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR -use, non_intrinsic :: cintrf_mod, only : COBJ +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 +use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: bobyqa_mod, only : bobyqa implicit none @@ -36,6 +37,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: npt integer(C_INT), intent(in), value :: iprint +type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info ! Local variables @@ -51,6 +53,8 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) +procedure(CCALLBACK), pointer :: cb_ptr +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -62,10 +66,19 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code -call bobyqa(calfun, x_loc, f_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, info=info_loc) +if (C_ASSOCIATED(callback_ptr)) then + ! If a C callback function is provided, we capture it for use in the closure below + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! And then we pass the closure to the Fortran code + call bobyqa(calfun, x_loc, f_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, callback_fcn=callback_fcn, info=info_loc) +else + call bobyqa(calfun, x_loc, f_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, info=info_loc) +end if ! Write the outputs x = real(x_loc, kind(x)) @@ -82,13 +95,69 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj +use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) real(RP), intent(out) :: f_sub -call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub) + +! Local variables +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc + +! Read the inputs and convert them to the types specified in COBJ +x_sub_loc = real(x_sub, kind(x_sub_loc)) + +! Call the C objective function +call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) + +! Write the output +f_sub = real(f_sub_loc, kind(f_sub)) + end subroutine calfun +! We name some variables _sub to avoid masking the parent variables +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +use, non_intrinsic :: consts_mod, only : RP, IK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL +implicit none +real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: f_sub +integer(IK), intent(in) :: nf_sub +integer(IK), intent(in) :: tr +real(RP), intent(in) :: cstrv +real(RP), intent(in) :: nlconstr(:) +logical, intent(out) :: terminate + +! Local variables +integer(C_INT) :: n_sub_loc +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +integer(C_INT) :: nf_sub_loc +integer(C_INT) :: tr_loc +real(C_DOUBLE) :: cstrv_loc +integer(C_INT) :: m_nlconstr +real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +logical(C_BOOL) :: terminate_loc + +! Read the inputs and convert them to the types specified in CCALLBACK +n_sub_loc = size(x_sub) +x_sub_loc = real(x_sub, kind(x_sub_loc)) +f_sub_loc = real(f_sub, kind(f_sub_loc)) +nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) +tr_loc = int(tr, kind(tr_loc)) +cstrv_loc = real(cstrv, kind(cstrv_loc)) +m_nlconstr = size(nlconstr) +nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Call the C objective function +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) + +! Write the output +terminate = logical(terminate_loc, kind(terminate)) + +end subroutine callback_fcn + + end subroutine bobyqa_c diff --git a/c/cintrf.f90 b/c/cintrf.f90 index f837dcc9f7..e66ef17178 100644 --- a/c/cintrf.f90 +++ b/c/cintrf.f90 @@ -7,7 +7,7 @@ module cintrf_mod implicit none private -public :: COBJ, COBJCON, evalcobj, evalcobjcon +public :: COBJ, COBJCON, CCALLBACK abstract interface @@ -31,66 +31,23 @@ subroutine COBJCON(x, f, constr, data_ptr) bind(c) type(C_PTR), intent(in), value :: data_ptr end subroutine COBJCON -end interface - - -contains - - -subroutine evalcobj(cobj_ptr, data_ptr, x, f) -use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_FUNPTR, C_F_PROCPOINTER, C_PTR -implicit none -type(C_FUNPTR), intent(in) :: cobj_ptr -type(C_PTR), intent(in), value :: data_ptr -real(RP), intent(in) :: x(:) -real(RP), intent(out) :: f - -! Local variables -procedure(COBJ), pointer :: obj_ptr -real(C_DOUBLE) :: x_loc(size(x)) -real(C_DOUBLE) :: f_loc - -! Read the inputs and convert them to the types specified in COBJ -x_loc = real(x, kind(x_loc)) -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) - -! Call the C objective function -call obj_ptr(x_loc, f_loc, data_ptr) - -! Write the output -f = real(f_loc, kind(f)) - -end subroutine evalcobj - - -subroutine evalcobjcon(cobjcon_ptr, data_ptr, x, f, constr) -use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_FUNPTR, C_F_PROCPOINTER, C_PTR -implicit none -type(C_FUNPTR), intent(in) :: cobjcon_ptr -type(C_PTR), intent(in), value :: data_ptr -real(RP), intent(in) :: x(:) -real(RP), intent(out) :: f -real(RP), intent(out) :: constr(:) - -! Local variables -procedure(COBJCON), pointer :: objcon_ptr -real(C_DOUBLE) :: x_loc(size(x)) -real(C_DOUBLE) :: f_loc -real(C_DOUBLE) :: constr_loc(size(constr)) - -! Read the inputs and convert them to the types specified in COBJCON -x_loc = real(x, kind(x_loc)) -call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) + subroutine CCALLBACK(n, x, f, nf, tr, cstrv, m_nlcon, nlconstr, terminate) bind(c) + use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_BOOL, C_INT + implicit none + integer(C_INT), intent(in), value :: n + ! We cannot use assumed-shape arrays for C interoperability + real(C_DOUBLE), intent(in) :: x(*) + real(C_DOUBLE), intent(in), value :: f + integer(C_INT), intent(in), value :: nf + integer(C_INT), intent(in), value :: tr + real(C_DOUBLE), intent(in), value :: cstrv + integer(C_INT), intent(in), value :: m_nlcon + real(C_DOUBLE), intent(in) :: nlconstr(*) + logical(C_BOOL), intent(out) :: terminate + end subroutine CCALLBACK -! Call the C objective function -call objcon_ptr(x_loc, f_loc, constr_loc, data_ptr) -! Write the output -f = real(f_loc, kind(f)) -constr = real(constr_loc, kind(constr)) +end interface -end subroutine evalcobjcon end module cintrf_mod diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 35714a1884..8fed66b2e1 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -12,10 +12,10 @@ module cobyla_c_mod contains -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, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR -use, non_intrinsic :: cintrf_mod, only : COBJCON +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) +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 : COBJCON, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: cobyla_mod, only : cobyla implicit none @@ -46,6 +46,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 +type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info ! Local variables @@ -69,6 +70,8 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP) :: xu_loc(n) real(RP) :: f0_loc real(RP) :: nlconstr0_loc(m_nlcon) +procedure(CCALLBACK), pointer :: cb_ptr +procedure(COBJCON), pointer :: objcon_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -88,13 +91,25 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) m_nlcon_loc = int(m_nlcon, kind(m_nlcon_loc)) +call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) ! Call the Fortran code -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, & - & iprint=iprint_loc, info=info_loc) +if (C_ASSOCIATED(callback_ptr)) then + ! If a C callback function is provided, we capture it for use in the closure below + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! And then we pass the closure to the Fortran code + 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, & + & 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, & + & iprint=iprint_loc, info=info_loc) +end if ! Write the outputs x = real(x_loc, kind(x)) @@ -113,14 +128,73 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ !--------------------------------------------------------------------------------------------------! subroutine calcfc(x_sub, f_sub, constr_sub) use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobjcon +use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) real(RP), intent(out) :: f_sub real(RP), intent(out) :: constr_sub(:) -call evalcobjcon(cobjcon_ptr, data_ptr, x_sub, f_sub, constr_sub) + +! Local variables +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +real(C_DOUBLE) :: constr_sub_loc(size(constr_sub)) + +! Read the inputs and convert them to the types specified in COBJCON +x_sub_loc = real(x_sub, kind(x_sub_loc)) + +! Call the C objective function +call objcon_ptr(x_sub_loc, f_sub_loc, constr_sub_loc, data_ptr) + +! Write the output +f_sub = real(f_sub_loc, kind(f_sub)) +constr_sub = real(constr_sub_loc, kind(constr_sub)) + end subroutine calcfc + +! We name some variables _sub to avoid masking the parent variables +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +use, non_intrinsic :: consts_mod, only : RP, IK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL +implicit none +real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: f_sub +integer(IK), intent(in) :: nf_sub +integer(IK), intent(in) :: tr +real(RP), intent(in) :: cstrv +real(RP), intent(in) :: nlconstr(:) +logical, intent(out) :: terminate + +! Local variables +integer(C_INT) :: n_sub_loc +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +integer(C_INT) :: nf_sub_loc +integer(C_INT) :: tr_loc +real(C_DOUBLE) :: cstrv_loc +integer(C_INT) :: m_nlconstr +real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +logical(C_BOOL) :: terminate_loc + +! Read the inputs and convert them to the types specified in CCALLBACK +n_sub_loc = size(x_sub) +x_sub_loc = real(x_sub, kind(x_sub_loc)) +f_sub_loc = real(f_sub, kind(f_sub_loc)) +nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) +tr_loc = int(tr, kind(tr_loc)) +cstrv_loc = real(cstrv, kind(cstrv_loc)) +m_nlconstr = size(nlconstr) +nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Call the C objective function +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) + +! Write the output +terminate = logical(terminate_loc, kind(terminate)) + +end subroutine callback_fcn + + end subroutine cobyla_c diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index d871944961..1b8cb617d0 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data) (void)data; } +static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) +{ + (void)n; + printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + *terminate = x[1] > 1.8; // test early termination + (void)m_nlcon; + (void)nlconstr; +} + int main(int argc, char * argv[]) { (void)argc; @@ -27,6 +36,7 @@ int main(int argc, char * argv[]) options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; + options.callback = &callback; prima_result_t result; const int rc = prima_minimize(PRIMA_BOBYQA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c index 8053eabc2b..ee0b0f6baa 100644 --- a/c/examples/cobyla/cobyla_example.c +++ b/c/examples/cobyla/cobyla_example.c @@ -16,6 +16,17 @@ static void fun(const double x[], double *f, double constr[], const void *data) (void)data; } + +static void callback(const int n, const double x[], const double f, const int nf, const int tr, const double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) +{ + (void)n; + printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + *terminate = 0; + (void)m_nlcon; + (void)nlconstr; +} + + int main(int argc, char * argv[]) { (void)argc; @@ -31,6 +42,7 @@ int main(int argc, char * argv[]) options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; + options.callback = &callback; problem.m_nlcon = M_NLCON; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c index e5aa89d05e..a1b4626ada 100644 --- a/c/examples/lincoa/lincoa_example.c +++ b/c/examples/lincoa/lincoa_example.c @@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data) (void)data; } +static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate) +{ + (void)n; + printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + (void)m_nlcon; + (void)nlconstr; + *terminate = 0; +} + int main(int argc, char * argv[]) { (void)argc; @@ -27,6 +36,7 @@ int main(int argc, char * argv[]) options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; + options.callback = &callback; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; double Aineq[3*2] = {1.0, 0.0, diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c index 00c4aa709b..c05d370b45 100644 --- a/c/examples/newuoa/newuoa_example.c +++ b/c/examples/newuoa/newuoa_example.c @@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data) (void)data; } +static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate) +{ + (void)n; + printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + (void)m_nlcon; + (void)nlconstr; + *terminate = 0; +} + int main(int argc, char * argv[]) { (void)argc; @@ -27,6 +36,7 @@ int main(int argc, char * argv[]) options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; + options.callback = &callback; prima_result_t result; const int rc = prima_minimize(PRIMA_NEWUOA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c index be53e87c21..90f4e02c22 100644 --- a/c/examples/uobyqa/uobyqa_example.c +++ b/c/examples/uobyqa/uobyqa_example.c @@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data) (void)data; } +static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) +{ + (void)n; + printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + *terminate = 0; + (void)m_nlcon; + (void)nlconstr; +} + int main(int argc, char * argv[]) { (void)argc; @@ -27,6 +36,7 @@ int main(int argc, char * argv[]) options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; + options.callback = &callback; prima_result_t result; const int rc = prima_minimize(PRIMA_UOBYQA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h index bb38c17961..babc335b2f 100644 --- a/c/include/prima/prima.h +++ b/c/include/prima/prima.h @@ -3,6 +3,8 @@ #ifndef PRIMA_H #define PRIMA_H +#include + #ifdef __cplusplus extern "C" { #endif @@ -49,6 +51,7 @@ typedef enum PRIMA_NO_SPACE_BETWEEN_BOUNDS = 6, PRIMA_DAMAGING_ROUNDING = 7, PRIMA_ZERO_LINEAR_CONSTRAINT = 8, + PRIMA_USER_STOP = 9, PRIMA_INVALID_INPUT = 100, PRIMA_ASSERTION_FAILS = 101, PRIMA_VALIDATION_FAILS = 102, @@ -67,7 +70,7 @@ PRIMAC_API const char *prima_get_rc_string(const prima_rc_t rc); /* - * A function as required by solvers + * The objective function as required by solvers * * x : on input, then vector of variables (should not be modified) * f : on output, the value of the function @@ -80,6 +83,21 @@ const char *prima_get_rc_string(const prima_rc_t rc); typedef void (*prima_obj_t)(const double x[], double *f, const void *data); typedef void (*prima_objcon_t)(const double x[], double *f, double constr[], const void *data); +/* An optional callback function to report algorithm progress + * + * n : number of variables + * x : the current best point + * f : the function value of the best point + * nf : number of objective function calls + * tr : iteration number + * cstrv : the constraint value verified by the current best point + * m_nlcon : number of non-linear constraints (cobyla only) + * nlconstr : non-linear constraints values verified by the current best point (cobyla only) + * terminate : a boolean to ask from early optimization exit +*/ +typedef void (*prima_callback_t)(const int n, const double x[], const double f, int nf, int tr, + const double cstrv, int m_nlcon, const double nlconstr[], bool *terminate); + typedef struct { @@ -105,6 +123,9 @@ typedef struct { // user-data, will be passed through the objective function callback void *data; + // callback function to report algorithm progress (default=NULL) + prima_callback_t callback; + } prima_options_t; /* Initialize problem */ diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index b41746394d..c4107c7dad 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -14,9 +14,9 @@ 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, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR -use, non_intrinsic :: cintrf_mod, only : COBJ + & 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, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: lincoa_mod, only : lincoa implicit none @@ -44,6 +44,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 +type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info ! Local variables @@ -64,6 +65,8 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) +procedure(CCALLBACK), pointer :: cb_ptr +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -81,13 +84,28 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code -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, info=info_loc) +if (C_ASSOCIATED(callback_ptr)) then + ! If a C callback function is provided, we capture it for use in the closure below + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! And then we pass the closure to the Fortran code + 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, 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, info=info_loc) +end if + ! Write the outputs x = real(x_loc, kind(x)) @@ -105,13 +123,69 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj +use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) real(RP), intent(out) :: f_sub -call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub) + +! Local variables +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc + +! Read the inputs and convert them to the types specified in COBJ +x_sub_loc = real(x_sub, kind(x_sub_loc)) + +! Call the C objective function +call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) + +! Write the output +f_sub = real(f_sub_loc, kind(f_sub)) + end subroutine calfun +! We name some variables _sub to avoid masking the parent variables +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +use, non_intrinsic :: consts_mod, only : RP, IK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL +implicit none +real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: f_sub +integer(IK), intent(in) :: nf_sub +integer(IK), intent(in) :: tr +real(RP), intent(in) :: cstrv +real(RP), intent(in) :: nlconstr(:) +logical, intent(out) :: terminate + +! Local variables +integer(C_INT) :: n_sub_loc +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +integer(C_INT) :: nf_sub_loc +integer(C_INT) :: tr_loc +real(C_DOUBLE) :: cstrv_loc +integer(C_INT) :: m_nlconstr +real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +logical(C_BOOL) :: terminate_loc + +! Read the inputs and convert them to the types specified in CCALLBACK +n_sub_loc = size(x_sub) +x_sub_loc = real(x_sub, kind(x_sub_loc)) +f_sub_loc = real(f_sub, kind(f_sub_loc)) +nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) +tr_loc = int(tr, kind(tr_loc)) +cstrv_loc = real(cstrv, kind(cstrv_loc)) +m_nlconstr = size(nlconstr) +nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Call the C objective function +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) + +! Write the output +terminate = logical(terminate_loc, kind(terminate)) + +end subroutine callback_fcn + + end subroutine lincoa_c diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index 5177d1970c..0d1d265cc2 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -12,9 +12,9 @@ module newuoa_c_mod contains -subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR -use, non_intrinsic :: cintrf_mod, only : COBJ +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, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: newuoa_mod, only : newuoa implicit none @@ -33,6 +33,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: npt integer(C_INT), intent(in), value :: iprint +type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info ! Local variables @@ -46,6 +47,8 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) +procedure(CCALLBACK), pointer :: cb_ptr +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -55,10 +58,19 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code -call newuoa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & - & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, info=info_loc) +if (C_ASSOCIATED(callback_ptr)) then + ! If a C callback function is provided, we capture it for use in the closure below + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! And then we pass the closure to the Fortran code + call newuoa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & + & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) +else + call newuoa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & + & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, info=info_loc) +end if ! Write the outputs x = real(x_loc, kind(x)) @@ -75,13 +87,68 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj +use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) real(RP), intent(out) :: f_sub -call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub) + +! Local variables +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc + +! Read the inputs and convert them to the types specified in COBJ +x_sub_loc = real(x_sub, kind(x_sub_loc)) + +! Call the C objective function +call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) + +! Write the output +f_sub = real(f_sub_loc, kind(f_sub)) + end subroutine calfun +! We name some variables _sub to avoid masking the parent variables +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +use, non_intrinsic :: consts_mod, only : RP, IK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL +implicit none +real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: f_sub +integer(IK), intent(in) :: nf_sub +integer(IK), intent(in) :: tr +real(RP), intent(in) :: cstrv +real(RP), intent(in) :: nlconstr(:) +logical, intent(out) :: terminate + +! Local variables +integer(C_INT) :: n_sub_loc +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +integer(C_INT) :: nf_sub_loc +integer(C_INT) :: tr_loc +real(C_DOUBLE) :: cstrv_loc +integer(C_INT) :: m_nlconstr +real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +logical(C_BOOL) :: terminate_loc + +! Read the inputs and convert them to the types specified in CCALLBACK +n_sub_loc = size(x_sub) +x_sub_loc = real(x_sub, kind(x_sub_loc)) +f_sub_loc = real(f_sub, kind(f_sub_loc)) +nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) +tr_loc = int(tr, kind(tr_loc)) +cstrv_loc = real(cstrv, kind(cstrv_loc)) +m_nlconstr = size(nlconstr) +nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Call the C objective function +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) + +! Write the output +terminate = logical(terminate_loc, kind(terminate)) + +end subroutine callback_fcn + end subroutine newuoa_c diff --git a/c/prima.c b/c/prima.c index 38d2d96127..a7ecfa75b2 100644 --- a/c/prima.c +++ b/c/prima.c @@ -44,19 +44,19 @@ 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[], double f0, const double nlconstr0[], - int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, int *info); + 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, int *info); + 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); int newuoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f, - int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, int *info); + 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); int uobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f, - int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, int *info); + 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 lincoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f, double *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 *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, int *info); + 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); int prima_check_problem(prima_problem_t *problem, prima_options_t *options, int alloc_bounds, int use_constr) { @@ -203,27 +203,27 @@ int prima_minimize(const prima_algorithm_t algorithm, prima_problem_t *problem, switch (algorithm) { case PRIMA_BOBYQA: - bobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, &info); + bobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info); break; 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, &info); + problem->xl, problem->xu, problem->f0, problem->nlconstr0, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, 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, &info); + problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info); break; case PRIMA_NEWUOA: - newuoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, &info); + newuoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info); break; case PRIMA_UOBYQA: - uobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, &info); + uobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, options->callback, &info); break; default: @@ -261,6 +261,8 @@ const char *prima_get_rc_string(const prima_rc_t rc) return "Rounding errors are becoming damaging"; case PRIMA_ZERO_LINEAR_CONSTRAINT: return "One of the linear constraints has a zero gradient"; + case PRIMA_USER_STOP: + return "User requested end of optimization"; case PRIMA_INVALID_INPUT: return "Invalid input"; case PRIMA_ASSERTION_FAILS: diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index 07d920dc13..ff117a8ba8 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -12,9 +12,9 @@ module uobyqa_c_mod contains -subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, iprint, info) bind(C) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR -use, non_intrinsic :: cintrf_mod, only : COBJ +subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, 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 +use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: uobyqa_mod, only : uobyqa implicit none @@ -22,6 +22,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr type(C_PTR), intent(in), value :: data_ptr + integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -32,6 +33,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(C_DOUBLE), intent(in), value :: ftarget integer(C_INT), intent(in), value :: maxfun integer(C_INT), intent(in), value :: iprint +type(C_FUNPTR), intent(in), value :: callback_ptr integer(C_INT), intent(out) :: info ! Local variables @@ -44,6 +46,8 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) +procedure(CCALLBACK), pointer :: cb_ptr +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -52,10 +56,19 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ftarget_loc = real(ftarget, kind(ftarget_loc)) maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code -call uobyqa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & - & maxfun=maxfun_loc, iprint=iprint_loc, info=info_loc) +if (C_ASSOCIATED(callback_ptr)) then + ! If a C callback function is provided, we capture it for use in the closure below + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! And then we pass the closure to the Fortran code + call uobyqa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & + & maxfun=maxfun_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) +else + call uobyqa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & + & maxfun=maxfun_loc, iprint=iprint_loc, info=info_loc) +end if ! Write the outputs x = real(x_loc, kind(x)) @@ -72,13 +85,69 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj +use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) real(RP), intent(out) :: f_sub -call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub) + +! Local variables +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc + +! Read the inputs and convert them to the types specified in COBJ +x_sub_loc = real(x_sub, kind(x_sub_loc)) + +! Call the C objective function +call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) + +! Write the output +f_sub = real(f_sub_loc, kind(f_sub)) + end subroutine calfun +! We name some variables _sub to avoid masking the parent variables +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +use, non_intrinsic :: consts_mod, only : RP, IK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL +implicit none +real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: f_sub +integer(IK), intent(in) :: nf_sub +integer(IK), intent(in) :: tr +real(RP), intent(in) :: cstrv +real(RP), intent(in) :: nlconstr(:) +logical, intent(out) :: terminate + +! Local variables +integer(C_INT) :: n_sub_loc +real(C_DOUBLE) :: x_sub_loc(size(x_sub)) +real(C_DOUBLE) :: f_sub_loc +integer(C_INT) :: nf_sub_loc +integer(C_INT) :: tr_loc +real(C_DOUBLE) :: cstrv_loc +integer(C_INT) :: m_nlconstr +real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +logical(C_BOOL) :: terminate_loc + +! Read the inputs and convert them to the types specified in CCALLBACK +n_sub_loc = size(x_sub) +x_sub_loc = real(x_sub, kind(x_sub_loc)) +f_sub_loc = real(f_sub, kind(f_sub_loc)) +nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) +tr_loc = int(tr, kind(tr_loc)) +cstrv_loc = real(cstrv, kind(cstrv_loc)) +m_nlconstr = size(nlconstr) +nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Call the C objective function +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) + +! Write the output +terminate = logical(terminate_loc, kind(terminate)) + +end subroutine callback_fcn + + end subroutine uobyqa_c diff --git a/fortran/bobyqa/bobyqa.f90 b/fortran/bobyqa/bobyqa.f90 index 24e7f71bc5..6662834203 100644 --- a/fortran/bobyqa/bobyqa.f90 +++ b/fortran/bobyqa/bobyqa.f90 @@ -39,7 +39,7 @@ module bobyqa_mod subroutine bobyqa(calfun, x, & & f, xl, xu, & & nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, & - & eta1, eta2, gamma1, gamma2, xhist, fhist, maxhist, honour_x0, info) + & eta1, eta2, gamma1, gamma2, xhist, fhist, maxhist, honour_x0, callback_fcn, info) !--------------------------------------------------------------------------------------------------! ! Among all the arguments, only CALFUN and X are obligatory. The others are OPTIONAL and you can ! neglect them unless you are familiar with the algorithm. Any unspecified optional input will take @@ -168,6 +168,9 @@ subroutine bobyqa(calfun, x, & ! if this requirement is not met. If HONOUR_X0 == TRUE, revise RHOBEG if needed; otherwise, revise ! X0 if needed. See the PREPROC subroutine for more information. ! +! CALLBACK +! Input, function to report progress and request termination. +! ! INFO ! Output, INTEGER(IK) scalar. ! INFO is the exit flag. It will be set to one of the following values defined in the module @@ -198,7 +201,7 @@ subroutine bobyqa(calfun, x, & use, non_intrinsic :: infos_mod, only : NO_SPACE_BETWEEN_BOUNDS use, non_intrinsic :: linalg_mod, only : trueloc use, non_intrinsic :: memory_mod, only : safealloc -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: preproc_mod, only : preproc use, non_intrinsic :: string_mod, only : num2str @@ -212,6 +215,7 @@ subroutine bobyqa(calfun, x, & real(RP), intent(inout) :: x(:) ! X(N) ! Optional inputs +procedure(CALLBACK), optional :: callback_fcn integer(IK), intent(in), optional :: iprint integer(IK), intent(in), optional :: maxfun integer(IK), intent(in), optional :: maxhist @@ -415,9 +419,15 @@ subroutine bobyqa(calfun, x, & !-------------------- Call BOBYQB, which performs the real calculations. --------------------------! -call bobyqb(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, & - & gamma1_loc, gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, f_loc, & - & fhist_loc, xhist_loc, info_loc) +if (present(callback_fcn)) then + call bobyqb(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, & + & gamma1_loc, gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, f_loc, & + & fhist_loc, xhist_loc, info_loc, callback_fcn) +else + call bobyqb(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, & + & gamma1_loc, gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, f_loc, & + & fhist_loc, xhist_loc, info_loc) +end if !--------------------------------------------------------------------------------------------------! ! Write the outputs. diff --git a/fortran/bobyqa/bobyqb.f90 b/fortran/bobyqa/bobyqb.f90 index 57fb718116..87c11e08ab 100644 --- a/fortran/bobyqa/bobyqb.f90 +++ b/fortran/bobyqa/bobyqb.f90 @@ -44,7 +44,7 @@ module bobyqb_mod subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamma2, rhobeg, rhoend, & - & xl, xu, x, nf, f, fhist, xhist, info) + & xl, xu, x, nf, f, fhist, xhist, info, callback_fcn) !--------------------------------------------------------------------------------------------------! ! This subroutine performs the major calculations of BOBYQA. ! @@ -82,10 +82,10 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED +use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP use, non_intrinsic :: linalg_mod, only : norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: powalg_mod, only : quadinc, calden, calvlag!, errquad use, non_intrinsic :: ratio_mod, only : redrat use, non_intrinsic :: redrho_mod, only : redrho @@ -102,6 +102,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm implicit none ! Inputs +procedure(CALLBACK), optional :: callback_fcn procedure(OBJ) :: calfun ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER integer(IK), intent(in) :: iprint integer(IK), intent(in) :: maxfun @@ -150,6 +151,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm logical :: reduce_rho logical :: shortd logical :: small_trrad +logical :: terminate logical :: trfail logical :: ximproved real(RP) :: bmat(size(x), npt + size(x)) @@ -216,6 +218,18 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm call initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x, ij, kopt, nf, fhist, fval, & & sl, su, xbase, xhist, xpt, subinfo) +! Report the current best value, and check if user asks for early termination. +terminate = .false. +if (present(callback_fcn)) then + ! Since this algorithm is a feasible method, the constraint violation is always 0. + ! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set + ! to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + subinfo = USER_STOP + end if +endif + ! Initialize X and F according to KOPT. x = xinbd(xbase, xpt(:, kopt), xl, xu, sl, su) ! In precise arithmetic, X = XBASE + XOPT. f = fval(kopt) @@ -571,6 +585,19 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm call shiftbase(kopt, xbase, xpt, zmat, bmat, pq, hq) xbase = max(xl, min(xu, xbase)) end if + + ! Report the current best value, and check if user asks for early termination. + if (present(callback_fcn)) then + ! Since this algorithm is a feasible method, the constraint violation is always 0. + ! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set + ! to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + info = USER_STOP + exit + end if + end if + end do ! End of DO TR = 1, MAXTR. The iterative procedure ends. ! Return from the calculation, after trying the Newton-Raphson step if it has not been tried yet. diff --git a/fortran/cobyla/cobyla.f90 b/fortran/cobyla/cobyla.f90 index e9dca05b5d..acbd81e6f7 100644 --- a/fortran/cobyla/cobyla.f90 +++ b/fortran/cobyla/cobyla.f90 @@ -60,7 +60,7 @@ subroutine cobyla(calcfc, m_nlcon, x, & & xl, xu, & & f0, nlconstr0, & & nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, iprint, eta1, eta2, gamma1, gamma2, & - & xhist, fhist, chist, nlchist, maxhist, maxfilt, info) + & xhist, fhist, chist, nlchist, maxhist, maxfilt, callback_fcn, info) !--------------------------------------------------------------------------------------------------! ! Among all the arguments, only CALCFC, M_NLCON, and X are obligatory. The others are OPTIONAL and ! you can neglect them unless you are familiar with the algorithm. Any unspecified optional input @@ -234,6 +234,9 @@ subroutine cobyla(calcfc, m_nlcon, x, & ! the returned solution; default: MAXFILT_DFT (a value lower than MIN_MAXFILT is not recommended); ! see common/consts.F90 for the definitions of MAXFILT_DFT and MIN_MAXFILT. ! +! CALLBACK +! Input, function to report progress and request termination. +! ! INFO ! Output, INTEGER(IK) scalar. ! INFO is the exit flag. It will be set to one of the following values defined in the module @@ -264,7 +267,7 @@ subroutine cobyla(calcfc, m_nlcon, x, & use, non_intrinsic :: infos_mod, only : INVALID_INPUT use, non_intrinsic :: linalg_mod, only : trueloc, matprod use, non_intrinsic :: memory_mod, only : safealloc -use, non_intrinsic :: pintrf_mod, only : OBJCON +use, non_intrinsic :: pintrf_mod, only : OBJCON, CALLBACK use, non_intrinsic :: selectx_mod, only : isbetter use, non_intrinsic :: preproc_mod, only : preproc use, non_intrinsic :: string_mod, only : num2str @@ -280,6 +283,7 @@ subroutine cobyla(calcfc, m_nlcon, x, & integer(IK), intent(in) :: m_nlcon ! Number of constraints defined in CALCFC ! Optional inputs +procedure(CALLBACK), optional :: callback_fcn integer(IK), intent(in), optional :: iprint integer(IK), intent(in), optional :: maxfilt integer(IK), intent(in), optional :: maxfun @@ -601,10 +605,15 @@ subroutine cobyla(calcfc, m_nlcon, x, & !-------------------- Call COBYLB, which performs the real calculations. --------------------------! -!call cobylb(calcfc_internal, iprint_loc, maxfilt_loc, maxfun_loc, ctol_loc, cweight_loc, eta1_loc, eta2_loc, & -call cobylb(calcfc, iprint_loc, maxfilt_loc, maxfun_loc, amat, bvec, ctol_loc, cweight_loc, & +if (present(callback_fcn)) then + call cobylb(calcfc, iprint_loc, maxfilt_loc, maxfun_loc, amat, bvec, ctol_loc, cweight_loc, & + & eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, gamma2_loc, rhobeg_loc, rhoend_loc, constr_loc, & + & f_loc, x, nf_loc, chist_loc, conhist_loc, cstrv_loc, fhist_loc, xhist_loc, info_loc, callback_fcn) +else + call cobylb(calcfc, iprint_loc, maxfilt_loc, maxfun_loc, amat, bvec, ctol_loc, cweight_loc, & & eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, gamma2_loc, rhobeg_loc, rhoend_loc, constr_loc, & & f_loc, x, nf_loc, chist_loc, conhist_loc, cstrv_loc, fhist_loc, xhist_loc, info_loc) +end if !--------------------------------------------------------------------------------------------------! ! Deallocate variables not needed any more. Indeed, automatic allocation will take place at exit. diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index b83cef5d7c..2f21b053af 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -29,12 +29,12 @@ module cobylb_mod subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, eta1, eta2, ftarget, & - & gamma1, gamma2, rhobeg, rhoend, constr, f, x, nf, chist, conhist, cstrv, fhist, xhist, info) + & gamma1, gamma2, rhobeg, rhoend, constr, f, x, nf, chist, conhist, cstrv, fhist, xhist, info, callback_fcn) !--------------------------------------------------------------------------------------------------! ! This subroutine performs the actual calculations of COBYLA. ! ! IPRINT, MAXFILT, MAXFUN, MAXHIST, CTOL, CWEIGHT, ETA1, ETA2, FTARGET, GAMMA1, GAMMA2, RHOBEG, -! RHOEND, X, NF, F, XHIST, FHIST, CHIST, CONHIST, CSTRV and INFO are identical to the corresponding +! RHOEND, X, NF, F, XHIST, FHIST, CHIST, CONHIST, CSTRV, INFO and CALLBACK are identical to the corresponding ! arguments in subroutine COBYLA. !--------------------------------------------------------------------------------------------------! @@ -46,10 +46,10 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING, USER_STOP use, non_intrinsic :: linalg_mod, only : inprod, matprod, norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg -use, non_intrinsic :: pintrf_mod, only : OBJCON +use, non_intrinsic :: pintrf_mod, only : OBJCON, CALLBACK use, non_intrinsic :: ratio_mod, only : redrat use, non_intrinsic :: redrho_mod, only : redrho use, non_intrinsic :: selectx_mod, only : savefilt, selectx, isbetter @@ -63,6 +63,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et implicit none ! Inputs +procedure(CALLBACK), optional :: callback_fcn procedure(OBJCON) :: calcfc ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER integer(IK), intent(in) :: iprint integer(IK), intent(in) :: maxfilt @@ -119,6 +120,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et logical :: improve_geo logical :: reduce_rho logical :: shortd +logical :: terminate logical :: trfail logical :: ximproved real(RP) :: A(size(x), size(constr)) ! A contains the approximate gradient for the constraints @@ -213,6 +215,15 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et call initxfc(calcfc_internal, iprint, maxfun, constr, ctol, f, ftarget, rhobeg, x, nf, chist, conhist, & & conmat, cval, fhist, fval, sim, simi, xhist, evaluated, subinfo) +! Report the current best value, and check if user asks for early termination. +terminate = .false. +if (present(callback_fcn)) then + call callback_fcn(sim(:, n+1), fval(n+1), nf, 0, cval(n+1), conmat(:, n+1), terminate) + if (terminate) then + subinfo = USER_STOP + end if +end if + ! Initialize the filter, including XFILT, FFILT, CONFILT, CFILT, and NFILT. ! N.B.: The filter is used only when selecting which iterate to return. It does not interfere with ! the iterations. COBYLA is NOT a filter method but a trust-region method based on an L-infinity @@ -592,6 +603,15 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et end if end if ! End of IF (REDUCE_RHO). The procedure of reducing RHO ends. + ! Report the current best value, and check if user asks for early termination. + if (present(callback_fcn)) then + call callback_fcn(sim(:, n+1), fval(n+1), nf, tr, cval(n+1), conmat(:, n+1), terminate) + if (terminate) then + info = USER_STOP + exit + end if + end if + end do ! End of DO TR = 1, MAXTR. The iterative procedure ends. ! Return from the calculation, after trying the last trust-region step if it has not been tried yet. diff --git a/fortran/common/infos.f90 b/fortran/common/infos.f90 index 72698fd48f..95344a4664 100644 --- a/fortran/common/infos.f90 +++ b/fortran/common/infos.f90 @@ -25,12 +25,12 @@ module infos_mod public :: DAMAGING_ROUNDING public :: NO_SPACE_BETWEEN_BOUNDS public :: ZERO_LINEAR_CONSTRAINT +public :: USER_STOP public :: INVALID_INPUT public :: ASSERTION_FAILS public :: VALIDATION_FAILS public :: MEMORY_ALLOCATION_FAILS - integer(IK), parameter :: INFO_DFT = 0 integer(IK), parameter :: SMALL_TR_RADIUS = 0 integer(IK), parameter :: FTARGET_ACHIEVED = 1 @@ -43,6 +43,7 @@ module infos_mod integer(IK), parameter :: NO_SPACE_BETWEEN_BOUNDS = 6 integer(IK), parameter :: DAMAGING_ROUNDING = 7 integer(IK), parameter :: ZERO_LINEAR_CONSTRAINT = 8 +integer(IK), parameter :: USER_STOP = 9 ! Stop-codes. ! The following codes are used by ERROR STOP as stop-codes, which should be default integers. diff --git a/fortran/common/message.f90 b/fortran/common/message.f90 index 779259e6a6..c254a3a005 100644 --- a/fortran/common/message.f90 +++ b/fortran/common/message.f90 @@ -33,7 +33,7 @@ subroutine retmsg(solver, info, iprint, nf, f, x, cstrv, constr) use, non_intrinsic :: fprint_mod, only : fprint use, non_intrinsic :: infos_mod, only : FTARGET_ACHIEVED, MAXFUN_REACHED, MAXTR_REACHED, & & SMALL_TR_RADIUS, TRSUBP_FAILED, NAN_INF_X, NAN_INF_F, NAN_INF_MODEL, DAMAGING_ROUNDING, & - & NO_SPACE_BETWEEN_BOUNDS, ZERO_LINEAR_CONSTRAINT + & NO_SPACE_BETWEEN_BOUNDS, ZERO_LINEAR_CONSTRAINT, USER_STOP use, non_intrinsic :: string_mod, only : strip, num2str implicit none @@ -124,6 +124,8 @@ subroutine retmsg(solver, info, iprint, nf, f, x, cstrv, constr) reason = 'there is no space between the lower and upper bounds of variable.' case (ZERO_LINEAR_CONSTRAINT) reason = 'one of the linear constraints has a zero gradient' +case (USER_STOP) + reason = 'user requested end of optimization' case default reason = 'UNKNOWN EXIT FLAG' end select diff --git a/fortran/common/pintrf.f90 b/fortran/common/pintrf.f90 index 8a610e0c60..3674f6dfdc 100644 --- a/fortran/common/pintrf.f90 +++ b/fortran/common/pintrf.f90 @@ -15,7 +15,7 @@ module pintrf_mod implicit none private -public :: OBJ, OBJCON +public :: OBJ, OBJCON, CALLBACK abstract interface @@ -36,7 +36,18 @@ subroutine OBJCON(x, f, constr) real(RP), intent(out) :: constr(:) end subroutine OBJCON -end interface + subroutine CALLBACK(x, f, nf, tr, cstrv, nlconstr, terminate) + use consts_mod, only : RP, IK + implicit none + real(RP), intent(in) :: x(:) + real(RP), intent(in) :: f + integer(IK), intent(in) :: nf + integer(IK), intent(in) :: tr + real(RP), intent(in) :: cstrv + real(RP), intent(in) :: nlconstr(:) + logical, intent(out) :: terminate + end subroutine CALLBACK +end interface end module pintrf_mod diff --git a/fortran/lincoa/lincoa.f90 b/fortran/lincoa/lincoa.f90 index c6aac8a3c7..e5a1574a29 100644 --- a/fortran/lincoa/lincoa.f90 +++ b/fortran/lincoa/lincoa.f90 @@ -53,7 +53,7 @@ subroutine lincoa(calfun, x, & & Aeq, beq, & & xl, xu, & & nf, rhobeg, rhoend, ftarget, ctol, cweight, maxfun, npt, iprint, eta1, eta2, gamma1, gamma2, & - & xhist, fhist, chist, maxhist, maxfilt, info) + & xhist, fhist, chist, maxhist, maxfilt, callback_fcn, info) !--------------------------------------------------------------------------------------------------! ! Among all the arguments, only CALFUN, and X are obligatory. The others are OPTIONAL and you can ! neglect them unless you are familiar with the algorithm. Any unspecified optional input will take @@ -195,6 +195,9 @@ subroutine lincoa(calfun, x, & ! CONSTS_MOD (see consts.F90 under the directory named "common"). ! Use *HIST with caution! (N.B.: the algorithm is NOT designed for large problems). ! +! CALLBACK +! Input, function to report progress and request termination. +! ! INFO ! Output, INTEGER(IK) scalar. ! INFO is the exit flag. It will be set to one of the following values defined in the module @@ -224,7 +227,7 @@ subroutine lincoa(calfun, x, & use, non_intrinsic :: history_mod, only : prehist use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf use, non_intrinsic :: memory_mod, only : safealloc -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: preproc_mod, only : preproc use, non_intrinsic :: selectx_mod, only : isbetter use, non_intrinsic :: string_mod, only : num2str @@ -239,6 +242,7 @@ subroutine lincoa(calfun, x, & real(RP), intent(inout) :: x(:) ! X(N) ! Optional inputs +procedure(CALLBACK), optional :: callback_fcn integer(IK), intent(in), optional :: iprint integer(IK), intent(in), optional :: maxfilt integer(IK), intent(in), optional :: maxfun @@ -509,10 +513,17 @@ subroutine lincoa(calfun, x, & call get_lincon(Aeq_loc, Aineq_loc, beq_loc, bineq_loc, rhoend_loc, xl_loc, xu_loc, x, amat, bvec) !-------------------- Call LINCOB, which performs the real calculations. --------------------------! -call lincob(calfun, iprint_loc, maxfilt_loc, maxfun_loc, npt_loc, Aeq_loc, Aineq_loc, amat, & - & beq_loc, bineq_loc, bvec, ctol_loc, cweight_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & - & gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, chist_loc, cstrv_loc, & - & f_loc, fhist_loc, xhist_loc, info_loc) +if (present(callback_fcn)) then + call lincob(calfun, iprint_loc, maxfilt_loc, maxfun_loc, npt_loc, Aeq_loc, Aineq_loc, amat, & + & beq_loc, bineq_loc, bvec, ctol_loc, cweight_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, chist_loc, cstrv_loc, & + & f_loc, fhist_loc, xhist_loc, info_loc, callback_fcn) +else + call lincob(calfun, iprint_loc, maxfilt_loc, maxfun_loc, npt_loc, Aeq_loc, Aineq_loc, amat, & + & beq_loc, bineq_loc, bvec, ctol_loc, cweight_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, xl_loc, xu_loc, x, nf_loc, chist_loc, cstrv_loc, & + & f_loc, fhist_loc, xhist_loc, info_loc) +end if !--------------------------------------------------------------------------------------------------! ! Deallocate variables not needed any more. Indeed, automatic allocation will take place at exit. diff --git a/fortran/lincoa/lincob.f90 b/fortran/lincoa/lincob.f90 index 1efdf7e33d..446ed4892c 100644 --- a/fortran/lincoa/lincob.f90 +++ b/fortran/lincoa/lincob.f90 @@ -28,7 +28,7 @@ module lincob_mod subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, bineq, bvec, & & ctol, cweight, eta1, eta2, ftarget, gamma1, gamma2, rhobeg, rhoend, xl, xu, x, nf, chist, & - & cstrv, f, fhist, xhist, info) + & cstrv, f, fhist, xhist, info, callback_fcn) !--------------------------------------------------------------------------------------------------! ! This subroutine performs the actual calculations of LINCOA. ! @@ -78,11 +78,11 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP use, non_intrinsic :: linalg_mod, only : matprod, maximum, eye, trueloc, linspace, norm, trueloc use, non_intrinsic :: memory_mod, only : safealloc use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: powalg_mod, only : quadinc, omega_mul, hess_mul, updateh use, non_intrinsic :: ratio_mod, only : redrat use, non_intrinsic :: redrho_mod, only : redrho @@ -98,6 +98,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b implicit none ! Inputs +procedure(CALLBACK), optional :: callback_fcn procedure(OBJ) :: calfun ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER integer(IK), intent(in) :: iprint integer(IK), intent(in) :: maxfilt @@ -169,6 +170,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b logical :: reduce_rho logical :: shortd logical :: small_trrad +logical :: terminate logical :: trfail logical :: ximproved real(RP) :: b(size(bvec)) @@ -253,6 +255,17 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b call initxf(calfun, iprint, maxfun, Aeq, Aineq, amat, beq, bineq, ctol, ftarget, rhobeg, xl, xu, & & x, b, ij, kopt, nf, chist, cval, fhist, fval, xbase, xhist, xpt, evaluated, subinfo) +! Report the current best value, and check if user asks for early termination. +terminate = .false. +if (present(callback_fcn)) then + ! Since we only solve linearly constrained problems, the nonlinear constraint is set + ! to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), [real(RP) ::], terminate) + if (terminate) then + subinfo = USER_STOP + end if +end if + ! Initialize X, F, CONSTR, and CSTRV according to KOPT. ! N.B.: We must set CONSTR and CSTRV. Otherwise, if REDUCE_RHO is TRUE after the very first ! iteration due to SHORTD, then RHOMSG will be called with CONSTR and CSTRV uninitialized. @@ -635,6 +648,18 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b pqalt = omega_mul(idz, zmat, fval - fval(kopt)) galt = matprod(bmat(:, 1:npt), fval - fval(kopt)) + hess_mul(xpt(:, kopt), xpt, pqalt) end if + + ! Report the current best value, and check if user asks for early termination. + if (present(callback_fcn)) then + ! Since we only solve linearly constrained problems, the nonlinear constraint is set + ! to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), [real(RP) ::], terminate) + if (terminate) then + info = USER_STOP + exit + end if + end if + end do ! End of DO TR = 1, MAXTR. The iterative procedure ends. ! Return from the calculation, after trying the Newton-Raphson step if it has not been tried yet. diff --git a/fortran/newuoa/newuoa.f90 b/fortran/newuoa/newuoa.f90 index f2c7b2d16f..f53222310e 100644 --- a/fortran/newuoa/newuoa.f90 +++ b/fortran/newuoa/newuoa.f90 @@ -35,7 +35,7 @@ module newuoa_mod subroutine newuoa(calfun, x, & & f, nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, eta1, eta2, gamma1, gamma2, & - & xhist, fhist, maxhist, info) + & xhist, fhist, maxhist, callback_fcn, info) !--------------------------------------------------------------------------------------------------! ! Among all the arguments, only CALFUN and X are obligatory. The others are OPTIONAL and you can ! neglect them unless you are familiar with the algorithm. Any unspecified optional input will take @@ -142,6 +142,9 @@ subroutine newuoa(calfun, x, & ! CONSTS_MOD (see consts.F90 under the directory named "common"). ! Use *HIST with caution! (N.B.: the algorithm is NOT designed for large problems). ! +! CALLBACK +! Input, function to report progress and request termination. +! ! INFO ! Output, INTEGER(IK) scalar. ! INFO is the exit flag. It will be set to one of the following values defined in the module @@ -169,7 +172,7 @@ subroutine newuoa(calfun, x, & use, non_intrinsic :: history_mod, only : prehist use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf use, non_intrinsic :: memory_mod, only : safealloc -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: preproc_mod, only : preproc use, non_intrinsic :: string_mod, only : num2str @@ -183,6 +186,7 @@ subroutine newuoa(calfun, x, & real(RP), intent(inout) :: x(:) ! Optional inputs +procedure(CALLBACK), optional :: callback_fcn integer(IK), intent(in), optional :: iprint integer(IK), intent(in), optional :: maxfun integer(IK), intent(in), optional :: maxhist @@ -333,8 +337,13 @@ subroutine newuoa(calfun, x, & !-------------------- Call NEWUOB, which performs the real calculations. --------------------------! -call newuob(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & - & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc) +if (present(callback_fcn)) then + call newuob(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc, callback_fcn) +else + call newuob(calfun, iprint_loc, maxfun_loc, npt_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc) +end if !--------------------------------------------------------------------------------------------------! diff --git a/fortran/newuoa/newuob.f90 b/fortran/newuoa/newuob.f90 index 0f6a2e61ec..1e67590fc8 100644 --- a/fortran/newuoa/newuob.f90 +++ b/fortran/newuoa/newuob.f90 @@ -20,7 +20,7 @@ module newuob_mod subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamma2, rhobeg, & - & rhoend, x, nf, f, fhist, xhist, info) + & rhoend, x, nf, f, fhist, xhist, info, callback_fcn) !--------------------------------------------------------------------------------------------------! ! This subroutine performs the actual calculations of NEWUOA. ! @@ -57,10 +57,10 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP use, non_intrinsic :: linalg_mod, only : norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: powalg_mod, only : quadinc, updateh use, non_intrinsic :: ratio_mod, only : redrat use, non_intrinsic :: redrho_mod, only : redrho @@ -75,6 +75,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm implicit none ! Inputs +procedure(CALLBACK), optional :: callback_fcn procedure(OBJ) :: calfun ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER integer(IK), intent(in) :: iprint integer(IK), intent(in) :: maxfun @@ -121,6 +122,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm logical :: reduce_rho logical :: shortd logical :: small_trrad +logical :: terminate logical :: trfail logical :: ximproved real(RP) :: bmat(size(x), npt + size(x)) @@ -175,6 +177,17 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Initialize XBASE, XPT, FVAL, and KOPT, together with the history, NF, and IJ. call initxf(calfun, iprint, maxfun, ftarget, rhobeg, x, ij, kopt, nf, fhist, fval, xbase, xhist, xpt, subinfo) +! Report the current best value, and check if user asks for early termination. +terminate = .false. +if (present(callback_fcn)) then + ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear + ! constraint is set to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + subinfo = USER_STOP + end if +end if + ! Initialize X and F according to KOPT. x = xbase + xpt(:, kopt) f = fval(kopt) @@ -574,6 +587,18 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm if (sum(xpt(:, kopt)**2) >= 1.0E2_RP * delta**2) then ! 1.0E2 works better than 1.0E3 on 20230227. call shiftbase(kopt, xbase, xpt, zmat, bmat, pq, hq, idz) end if + + ! Report the current best value, and check if user asks for early termination. + if (present(callback_fcn)) then + ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear + ! constraint is set to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + info = USER_STOP + exit + end if + end if + end do ! End of DO TR = 1, MAXTR. The iterative procedure ends. ! Return from the calculation, after trying the Newton-Raphson step if it has not been tried yet. diff --git a/fortran/uobyqa/uobyqa.f90 b/fortran/uobyqa/uobyqa.f90 index 55261e3398..297ec04673 100644 --- a/fortran/uobyqa/uobyqa.f90 +++ b/fortran/uobyqa/uobyqa.f90 @@ -32,7 +32,7 @@ module uobyqa_mod subroutine uobyqa(calfun, x, & & f, nf, rhobeg, rhoend, ftarget, maxfun, iprint, eta1, eta2, gamma1, gamma2, & - & xhist, fhist, maxhist, info) + & xhist, fhist, maxhist, callback_fcn, info) !--------------------------------------------------------------------------------------------------! ! Among all the arguments, only CALFUN and X are obligatory. The others are OPTIONAL and you can ! neglect them unless you are familiar with the algorithm. Any unspecified optional input will take @@ -134,6 +134,9 @@ subroutine uobyqa(calfun, x, & ! CONSTS_MOD (see consts.F90 under the directory named "common"). ! Use *HIST with caution!!! (N.B.: the algorithm is NOT designed for large problems). ! +! CALLBACK +! Input, function to report progress and request termination. +! ! INFO ! Output, INTEGER(IK) scalar. ! INFO is the exit flag. It will be set to one of the following values defined in the module @@ -161,7 +164,7 @@ subroutine uobyqa(calfun, x, & use, non_intrinsic :: history_mod, only : prehist use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf use, non_intrinsic :: memory_mod, only : safealloc -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: preproc_mod, only : preproc use, non_intrinsic :: string_mod, only : num2str @@ -175,6 +178,7 @@ subroutine uobyqa(calfun, x, & real(RP), intent(inout) :: x(:) ! X(N) ! Optional inputs +procedure(CALLBACK), optional :: callback_fcn integer(IK), intent(in), optional :: iprint integer(IK), intent(in), optional :: maxfun integer(IK), intent(in), optional :: maxhist @@ -316,8 +320,13 @@ subroutine uobyqa(calfun, x, & !-------------------- Call UOBYQB, which performs the real calculations. --------------------------! -call uobyqb(calfun, iprint_loc, maxfun_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & - & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc) +if (present(callback_fcn)) then + call uobyqb(calfun, iprint_loc, maxfun_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc, callback_fcn) +else + call uobyqb(calfun, iprint_loc, maxfun_loc, eta1_loc, eta2_loc, ftarget_loc, gamma1_loc, & + & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, info_loc) +end if !--------------------------------------------------------------------------------------------------! diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90 index 0f9c36015c..f854664707 100644 --- a/fortran/uobyqa/uobyqb.f90 +++ b/fortran/uobyqa/uobyqb.f90 @@ -20,7 +20,7 @@ module uobyqb_mod subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, rhobeg, rhoend, & - & x, nf, f, fhist, xhist, info) + & x, nf, f, fhist, xhist, info, callback_fcn) !--------------------------------------------------------------------------------------------------! ! This subroutine performs the major calculations of UOBYQA. ! @@ -50,11 +50,11 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED +use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP use, non_intrinsic :: linalg_mod, only : vec2smat, smat_mul_vec, norm use, non_intrinsic :: memory_mod, only : safealloc use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg -use, non_intrinsic :: pintrf_mod, only : OBJ +use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK use, non_intrinsic :: powalg_mod, only : quadinc use, non_intrinsic :: ratio_mod, only : redrat use, non_intrinsic :: redrho_mod, only : redrho @@ -69,6 +69,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r implicit none ! Inputs +procedure(CALLBACK), optional :: callback_fcn procedure(OBJ) :: calfun ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER integer(IK), intent(in) :: iprint integer(IK), intent(in) :: maxfun @@ -112,6 +113,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r logical :: reduce_rho logical :: shortd logical :: small_trrad +logical :: terminate logical :: trfail logical :: ximproved real(RP) :: crvmin @@ -165,6 +167,18 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r ! Initialize XBASE, XPT, FVAL, and KOPT, together with the history and NF. call initxf(calfun, iprint, maxfun, ftarget, rhobeg, x, kopt, nf, fhist, fval, xbase, xhist, xpt, subinfo) + +! Report the current best value, and check if user asks for early termination. +terminate = .false. +if (present(callback_fcn)) then + ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear + ! constraint is set to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + subinfo = USER_STOP + end if +end if + ! Initialize X and F according to KOPT. x = xbase + xpt(:, kopt) f = fval(kopt) @@ -458,6 +472,17 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r if (sum(xpt(:, kopt)**2) >= 1.0E3_RP * delta**2) then call shiftbase(kopt, pl, pq, xbase, xpt) end if + + ! Report the current best value, and check if user asks for early termination. + if (present(callback_fcn)) then + ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear + ! constraint is set to an empty array. + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + if (terminate) then + info = USER_STOP + exit + end if + end if end do ! End of DO TR = 1, MAXTR. The iterative procedure ends. ! Deallocate PL. F2003 automatically deallocate local ALLOCATABLE variables at exit, yet we prefer From 55985e1b59de9daa1abac024600c7cb92ee443d9 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Thu, 14 Dec 2023 10:43:14 +0800 Subject: [PATCH 02/13] Add gdb to tests run via CI in order to obtain more information upon segfault Also disabled parallelism in cmake.yml so that the build log is more readable. The extra time it takes is not a concern. --- .github/workflows/cmake.yml | 26 +++++++++++++++++-------- CMakeLists.txt | 9 +++++++++ c/CMakeLists.txt | 19 ++++++++++++++++++- c/tests/CMakeLists.txt | 38 ++++++++++++++++++++++++++++++++----- fortran/CMakeLists.txt | 19 ++++++++++++++++++- 5 files changed, 96 insertions(+), 15 deletions(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index adfe78991a..f459f2042e 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -83,6 +83,9 @@ jobs: ssh-key: ${{ secrets.SSH_PRIVATE_KEY_ACT }} # This forces checkout to use SSH, not HTTPS submodules: recursive + - name: Miscellaneous setup + run: bash .github/scripts/misc_setup + - name: Install Ninja / Ubuntu if: ${{ matrix.os == 'ubuntu-latest' }} run: sudo apt update && sudo apt install ninja-build @@ -116,9 +119,10 @@ jobs: run: | cmake --version cmake -G Ninja -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_INSTALL_PREFIX=. -LAH -DCMAKE_C_FLAGS="${{ matrix.toolchain.cflags }}" -DCMAKE_Fortran_FLAGS="${{ matrix.toolchain.fflags }}" . - cmake --build . --target install --parallel 4 - cmake --build . --target tests --parallel 4 - ctest --output-on-failure -V -j4 -E stress + cmake --build . --target install + cmake --build . --target tests + ctest --output-on-failure -V -E stress + shell: bash env: FC: ${{ steps.setup-fortran.outputs.fc }} shell: bash @@ -126,7 +130,8 @@ jobs: - name: Stress test if: ${{ github.event_name == 'schedule' || github.event.inputs.stress-test == 'true' }} run: | - ctest --output-on-failure -V -j4 -R stress + ctest --output-on-failure -V -R stress + shell: bash cmake-other: @@ -156,6 +161,9 @@ jobs: ssh-key: ${{ secrets.SSH_PRIVATE_KEY_ACT }} # This forces checkout to use SSH, not HTTPS submodules: recursive + - name: Miscellaneous setup + run: bash .github/scripts/misc_setup + - name: Install AOCC if: ${{ matrix.toolchain.compiler == 'aflang' }} run: bash .github/scripts/install_aocc @@ -171,17 +179,19 @@ jobs: - name: Build run: | cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_INSTALL_PREFIX=. -LAH -DCMAKE_C_FLAGS="${{ matrix.toolchain.cflags }}" -DCMAKE_Fortran_FLAGS="${{ matrix.toolchain.fflags }}" . - cmake --build . --target install --parallel 4 - cmake --build . --target tests --parallel 4 + cmake --build . --target install + cmake --build . --target tests # cobyla test does not pass on AOCC: https://github.com/libprima/prima/issues/41 - ctest --output-on-failure -V -j4 -E "stress|cobyla" + ctest --output-on-failure -V -E "stress|cobyla" + shell: bash env: FC: ${{ matrix.toolchain.compiler }} - name: Stress test if: ${{ github.event_name == 'schedule' || github.event.inputs.stress-test == 'true' }} run: | - ctest --output-on-failure -V -j4 -R stress -E cobyla + ctest --output-on-failure -V -R stress -E cobyla + shell: bash # The following job check whether the tests were successful or cancelled due to timeout. diff --git a/CMakeLists.txt b/CMakeLists.txt index 78fcd98166..8361a369d0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,15 @@ if (PRIMA_HEAP_ARRAYS) endif () endif () +# For running tests with gdb. $_exitcode == -1 means the program ran without exiting +# normally, and in this case we want to show a stack trace +file(WRITE ${CMAKE_BINARY_DIR}/cmdfile.gdb "init-if-undefined $_exitcode = -1 +run +if $_exitcode == -1 + where +end +quit $_exitcode") + option(PRIMA_ENABLE_EXAMPLES "build examples by default" OFF) add_custom_target (examples) enable_testing () diff --git a/c/CMakeLists.txt b/c/CMakeLists.txt index 2eadfb114b..4440796e3c 100644 --- a/c/CMakeLists.txt +++ b/c/CMakeLists.txt @@ -35,7 +35,24 @@ macro (prima_add_c_test name) if (WIN32) set_target_properties(example_${name}_c_exe PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin) endif() - add_test (NAME example_${name}_c COMMAND example_${name}_c_exe) + + # Outside of CI we don't want to force people to run examples with gdb, so we test the executables by themselves. + # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up + # before we set up the examples for CI + if(APPLE OR UNIX) + add_test (NAME example_${name}_c COMMAND example_${name}_c_exe) + elseif(WIN32) + add_test (NAME example_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_c_exe.exe) + endif() + + # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. + # Of course this can be run locally as well if you define CI in your environment. + if(NOT APPLE AND UNIX AND DEFINED ENV{CI}) # Apple security policy will not allow running gdb in CI + add_test (NAME example_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb example_${name}_c_exe) + elseif(WIN32 AND DEFINED ENV{CI}) + add_test (NAME example_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_c_exe.exe) + endif() + add_dependencies(examples example_${name}_c_exe) endmacro () diff --git a/c/tests/CMakeLists.txt b/c/tests/CMakeLists.txt index 6fc0d612df..ed7f678116 100644 --- a/c/tests/CMakeLists.txt +++ b/c/tests/CMakeLists.txt @@ -9,11 +9,39 @@ macro (prima_add_c_test_multi name) endif() target_link_libraries(${name}_c_exe PRIVATE primac) add_dependencies(tests ${name}_c_exe) -add_test(NAME bobyqa_${name}_c COMMAND ${name}_c_exe bobyqa) - add_test(NAME cobyla_${name}_c COMMAND ${name}_c_exe cobyla) - add_test(NAME lincoa_${name}_c COMMAND ${name}_c_exe lincoa) - add_test(NAME newuoa_${name}_c COMMAND ${name}_c_exe newuoa) - add_test(NAME uobyqa_${name}_c COMMAND ${name}_c_exe uobyqa) + # Outside of CI we don't want to force people to run tests with gdb, so we test the executables by themselves. + # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up + # before we set up the tests for CI + if(APPLE OR UNIX) + add_test(NAME bobyqa_${name}_c COMMAND ${name}_c_exe bobyqa) + add_test(NAME cobyla_${name}_c COMMAND ${name}_c_exe cobyla) + add_test(NAME lincoa_${name}_c COMMAND ${name}_c_exe lincoa) + add_test(NAME newuoa_${name}_c COMMAND ${name}_c_exe newuoa) + add_test(NAME uobyqa_${name}_c COMMAND ${name}_c_exe uobyqa) + elseif(WIN32) + add_test(NAME bobyqa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe bobyqa) + add_test(NAME cobyla_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe cobyla) + add_test(NAME lincoa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe lincoa) + add_test(NAME newuoa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe newuoa) + add_test(NAME uobyqa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe uobyqa) + endif() + + # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. + # Of course this can be run locally as well if you define CI in your environment. + if(NOT APPLE AND UNIX AND DEFINED ENV{CI}) # Apple security policy will not allow running gdb in CI + add_test(NAME bobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe bobyqa) + add_test(NAME cobyla_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe cobyla) + add_test(NAME lincoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe lincoa) + add_test(NAME newuoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe newuoa) + add_test(NAME uobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe uobyqa) + elseif(WIN32 AND DEFINED ENV{CI}) + add_test(NAME bobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe bobyqa) + add_test(NAME cobyla_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe cobyla) + add_test(NAME lincoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe lincoa) + add_test(NAME newuoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe newuoa) + add_test(NAME uobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe uobyqa) + endif() + endmacro () prima_add_c_test_multi(stress) diff --git a/fortran/CMakeLists.txt b/fortran/CMakeLists.txt index 1a926fd4b4..3d1eb05dfe 100644 --- a/fortran/CMakeLists.txt +++ b/fortran/CMakeLists.txt @@ -114,7 +114,24 @@ macro (prima_add_f_test name) if (PRIMA_ENABLE_EXAMPLES) set_target_properties (example_${name}_fortran_exe PROPERTIES EXCLUDE_FROM_ALL FALSE) endif () - add_test (NAME example_${name}_fortran COMMAND example_${name}_fortran_exe) + + # Outside of CI we don't want to force people to run examples with gdb, so we test the executables by themselves. + # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up + # before we set up the examples for CI + if(APPLE OR UNIX) + add_test (NAME example_${name}_fortran COMMAND example_${name}_fortran_exe) + elseif(WIN32) + add_test (NAME example_${name}_fortran COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_fortran_exe.exe) + endif() + + # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. + # Of course this can be run locally as well if you define CI in your environment. + if(NOT APPLE AND UNIX AND DEFINED ENV{CI}) # Apple security policy will not allow running gdb in CI + add_test (NAME example_${name}_fortran_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb example_${name}_fortran_exe) + elseif(WIN32 AND DEFINED ENV{CI}) + add_test (NAME example_${name}_fortran_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_fortran_exe.exe) + endif() + add_dependencies(examples example_${name}_fortran_exe) endmacro () From d7bc0f439c3484a21b846f552c1b2788966c2d8b Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Thu, 14 Dec 2023 23:01:16 +0800 Subject: [PATCH 03/13] Turned cstrv, nlconstr, and terminate into optional arguments THe flang-derived compiler were beginning to have problem with the empty array used for nlconstr for the non-COBYLA algorithms. By making them optional arguments, we are able to handle them in a different way which satisfies all the compilers. --- .github/actions/spelling/allow.txt | 2 ++ c/bobyqa_c.f90 | 56 ++++++++++++++++++++++-------- c/cobyla_c.f90 | 54 ++++++++++++++++++++-------- c/lincoa_c.f90 | 53 ++++++++++++++++++++-------- c/newuoa_c.f90 | 55 +++++++++++++++++++++-------- c/uobyqa_c.f90 | 53 ++++++++++++++++++++-------- fortran/bobyqa/bobyqb.f90 | 10 ++---- fortran/common/pintrf.f90 | 6 ++-- fortran/lincoa/lincob.f90 | 8 ++--- fortran/newuoa/newuob.f90 | 8 ++--- fortran/primaf-GNU.def | 2 ++ fortran/primaf-Intel.def | 2 ++ fortran/uobyqa/uobyqb.f90 | 8 ++--- 13 files changed, 217 insertions(+), 100 deletions(-) diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt index dfb561edf4..57e1859dc8 100644 --- a/.github/actions/spelling/allow.txt +++ b/.github/actions/spelling/allow.txt @@ -2136,3 +2136,5 @@ nouninit libgfortran chocolatey fcn +BINDIR +cmdfile diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 5dd07cf520..be5d22b708 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -92,12 +92,13 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub ! Local variables @@ -115,18 +116,25 @@ subroutine calfun(x_sub, f_sub) end subroutine calfun -! We name some variables _sub to avoid masking the parent variables -subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) + +!--------------------------------------------------------------------------------------------------! +! This subroutine defines `callback_fcn` using the C function pointer with an internal subroutine. +! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. +! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across all 5 algorithms. +!--------------------------------------------------------------------------------------------------! +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) use, non_intrinsic :: consts_mod, only : RP, IK +use, non_intrinsic :: memory_mod, only : safealloc use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub integer(IK), intent(in) :: nf_sub integer(IK), intent(in) :: tr -real(RP), intent(in) :: cstrv -real(RP), intent(in) :: nlconstr(:) -logical, intent(out) :: terminate +real(RP), intent(in), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: terminate ! Local variables integer(C_INT) :: n_sub_loc @@ -134,9 +142,9 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) real(C_DOUBLE) :: f_sub_loc integer(C_INT) :: nf_sub_loc integer(C_INT) :: tr_loc -real(C_DOUBLE) :: cstrv_loc +real(C_DOUBLE) :: cstrv_sub_loc integer(C_INT) :: m_nlconstr -real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc ! Read the inputs and convert them to the types specified in CCALLBACK @@ -145,18 +153,36 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) f_sub_loc = real(f_sub, kind(f_sub_loc)) nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) tr_loc = int(tr, kind(tr_loc)) -cstrv_loc = real(cstrv, kind(cstrv_loc)) -m_nlconstr = size(nlconstr) -nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Set the constraint violation to a sensible default value if it is not provided. +if (present(cstrv_sub)) then + cstrv_sub_loc = real(cstrv_sub, kind(cstrv_sub_loc)) +else + cstrv_sub_loc = 0.0_C_DOUBLE +end if + +! Set the nonlinear constraints to a sensible default value if it is not provided. +if (present(nlconstr_sub)) then + m_nlconstr = int(size(nlconstr_sub), C_INT) + call safealloc(nlconstr_sub_loc, int(m_nlconstr, IK)) + nlconstr_sub_loc = real(nlconstr_sub, kind(nlconstr_sub_loc)) +else + m_nlconstr = 0_C_INT + nlconstr_sub_loc = [real(C_DOUBLE) ::] +end if ! Call the C objective function -call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output -terminate = logical(terminate_loc, kind(terminate)) +if ( present(terminate) ) then + terminate = logical(terminate_loc, kind(terminate)) +end if -end subroutine callback_fcn +! Deallocate resources +if (allocated(nlconstr_sub_loc)) deallocate(nlconstr_sub_loc) +end subroutine callback_fcn end subroutine bobyqa_c diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 8fed66b2e1..2ba1fd59ce 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -125,12 +125,13 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! This subroutine defines `calcfc` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calcfc(x_sub, f_sub, constr_sub) use, non_intrinsic :: consts_mod, only : RP use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub real(RP), intent(out) :: constr_sub(:) @@ -152,18 +153,24 @@ subroutine calcfc(x_sub, f_sub, constr_sub) end subroutine calcfc -! We name some variables _sub to avoid masking the parent variables -subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) +!--------------------------------------------------------------------------------------------------! +! This subroutine defines `callback_fcn` using the C function pointer with an internal subroutine. +! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. +! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across all 5 algorithms. +!--------------------------------------------------------------------------------------------------! +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) use, non_intrinsic :: consts_mod, only : RP, IK +use, non_intrinsic :: memory_mod, only : safealloc use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub integer(IK), intent(in) :: nf_sub integer(IK), intent(in) :: tr -real(RP), intent(in) :: cstrv -real(RP), intent(in) :: nlconstr(:) -logical, intent(out) :: terminate +real(RP), intent(in), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: terminate ! Local variables integer(C_INT) :: n_sub_loc @@ -171,9 +178,9 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) real(C_DOUBLE) :: f_sub_loc integer(C_INT) :: nf_sub_loc integer(C_INT) :: tr_loc -real(C_DOUBLE) :: cstrv_loc +real(C_DOUBLE) :: cstrv_sub_loc integer(C_INT) :: m_nlconstr -real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc ! Read the inputs and convert them to the types specified in CCALLBACK @@ -182,15 +189,34 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) f_sub_loc = real(f_sub, kind(f_sub_loc)) nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) tr_loc = int(tr, kind(tr_loc)) -cstrv_loc = real(cstrv, kind(cstrv_loc)) -m_nlconstr = size(nlconstr) -nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Set the constraint violation to a sensible default value if it is not provided. +if (present(cstrv_sub)) then + cstrv_sub_loc = real(cstrv_sub, kind(cstrv_sub_loc)) +else + cstrv_sub_loc = 0.0_C_DOUBLE +end if + +! Set the nonlinear constraints to a sensible default value if it is not provided. +if (present(nlconstr_sub)) then + m_nlconstr = int(size(nlconstr_sub), C_INT) + call safealloc(nlconstr_sub_loc, int(m_nlconstr, IK)) + nlconstr_sub_loc = real(nlconstr_sub, kind(nlconstr_sub_loc)) +else + m_nlconstr = 0_C_INT + nlconstr_sub_loc = [real(C_DOUBLE) ::] +end if ! Call the C objective function -call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output -terminate = logical(terminate_loc, kind(terminate)) +if ( present(terminate) ) then + terminate = logical(terminate_loc, kind(terminate)) +end if + +! Deallocate resources +if (allocated(nlconstr_sub_loc)) deallocate(nlconstr_sub_loc) end subroutine callback_fcn diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index c4107c7dad..0ddb65a4d3 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -143,18 +143,25 @@ subroutine calfun(x_sub, f_sub) end subroutine calfun -! We name some variables _sub to avoid masking the parent variables -subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) + +!--------------------------------------------------------------------------------------------------! +! This subroutine defines `callback_fcn` using the C function pointer with an internal subroutine. +! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. +! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across all 5 algorithms. +!--------------------------------------------------------------------------------------------------! +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) use, non_intrinsic :: consts_mod, only : RP, IK +use, non_intrinsic :: memory_mod, only : safealloc use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub integer(IK), intent(in) :: nf_sub integer(IK), intent(in) :: tr -real(RP), intent(in) :: cstrv -real(RP), intent(in) :: nlconstr(:) -logical, intent(out) :: terminate +real(RP), intent(in), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: terminate ! Local variables integer(C_INT) :: n_sub_loc @@ -162,9 +169,9 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) real(C_DOUBLE) :: f_sub_loc integer(C_INT) :: nf_sub_loc integer(C_INT) :: tr_loc -real(C_DOUBLE) :: cstrv_loc +real(C_DOUBLE) :: cstrv_sub_loc integer(C_INT) :: m_nlconstr -real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc ! Read the inputs and convert them to the types specified in CCALLBACK @@ -173,18 +180,36 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) f_sub_loc = real(f_sub, kind(f_sub_loc)) nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) tr_loc = int(tr, kind(tr_loc)) -cstrv_loc = real(cstrv, kind(cstrv_loc)) -m_nlconstr = size(nlconstr) -nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Set the constraint violation to a sensible default value if it is not provided. +if (present(cstrv_sub)) then + cstrv_sub_loc = real(cstrv_sub, kind(cstrv_sub_loc)) +else + cstrv_sub_loc = 0.0_C_DOUBLE +end if + +! Set the nonlinear constraints to a sensible default value if it is not provided. +if (present(nlconstr_sub)) then + m_nlconstr = int(size(nlconstr_sub), C_INT) + call safealloc(nlconstr_sub_loc, int(m_nlconstr, IK)) + nlconstr_sub_loc = real(nlconstr_sub, kind(nlconstr_sub_loc)) +else + m_nlconstr = 0_C_INT + nlconstr_sub_loc = [real(C_DOUBLE) ::] +end if ! Call the C objective function -call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output -terminate = logical(terminate_loc, kind(terminate)) +if ( present(terminate) ) then + terminate = logical(terminate_loc, kind(terminate)) +end if -end subroutine callback_fcn +! Deallocate resources +if (allocated(nlconstr_sub_loc)) deallocate(nlconstr_sub_loc) +end subroutine callback_fcn end subroutine lincoa_c diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index 0d1d265cc2..b079a03fb0 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -84,12 +84,13 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, non_intrinsic :: consts_mod, only : RP use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub ! Local variables @@ -107,18 +108,25 @@ subroutine calfun(x_sub, f_sub) end subroutine calfun -! We name some variables _sub to avoid masking the parent variables -subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) + +!--------------------------------------------------------------------------------------------------! +! This subroutine defines `callback_fcn` using the C function pointer with an internal subroutine. +! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. +! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across all 5 algorithms. +!--------------------------------------------------------------------------------------------------! +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) use, non_intrinsic :: consts_mod, only : RP, IK +use, non_intrinsic :: memory_mod, only : safealloc use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub integer(IK), intent(in) :: nf_sub integer(IK), intent(in) :: tr -real(RP), intent(in) :: cstrv -real(RP), intent(in) :: nlconstr(:) -logical, intent(out) :: terminate +real(RP), intent(in), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: terminate ! Local variables integer(C_INT) :: n_sub_loc @@ -126,9 +134,9 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) real(C_DOUBLE) :: f_sub_loc integer(C_INT) :: nf_sub_loc integer(C_INT) :: tr_loc -real(C_DOUBLE) :: cstrv_loc +real(C_DOUBLE) :: cstrv_sub_loc integer(C_INT) :: m_nlconstr -real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc ! Read the inputs and convert them to the types specified in CCALLBACK @@ -137,15 +145,34 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) f_sub_loc = real(f_sub, kind(f_sub_loc)) nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) tr_loc = int(tr, kind(tr_loc)) -cstrv_loc = real(cstrv, kind(cstrv_loc)) -m_nlconstr = size(nlconstr) -nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Set the constraint violation to a sensible default value if it is not provided. +if (present(cstrv_sub)) then + cstrv_sub_loc = real(cstrv_sub, kind(cstrv_sub_loc)) +else + cstrv_sub_loc = 0.0_C_DOUBLE +end if + +! Set the nonlinear constraints to a sensible default value if it is not provided. +if (present(nlconstr_sub)) then + m_nlconstr = int(size(nlconstr_sub), C_INT) + call safealloc(nlconstr_sub_loc, int(m_nlconstr, IK)) + nlconstr_sub_loc = real(nlconstr_sub, kind(nlconstr_sub_loc)) +else + m_nlconstr = 0_C_INT + nlconstr_sub_loc = [real(C_DOUBLE) ::] +end if ! Call the C objective function -call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output -terminate = logical(terminate_loc, kind(terminate)) +if ( present(terminate) ) then + terminate = logical(terminate_loc, kind(terminate)) +end if + +! Deallocate resources +if (allocated(nlconstr_sub_loc)) deallocate(nlconstr_sub_loc) end subroutine callback_fcn diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index ff117a8ba8..9c744b0124 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -105,18 +105,25 @@ subroutine calfun(x_sub, f_sub) end subroutine calfun -! We name some variables _sub to avoid masking the parent variables -subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) + +!--------------------------------------------------------------------------------------------------! +! This subroutine defines `callback_fcn` using the C function pointer with an internal subroutine. +! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. +! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across all 5 algorithms. +!--------------------------------------------------------------------------------------------------! +subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) use, non_intrinsic :: consts_mod, only : RP, IK +use, non_intrinsic :: memory_mod, only : safealloc use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub integer(IK), intent(in) :: nf_sub integer(IK), intent(in) :: tr -real(RP), intent(in) :: cstrv -real(RP), intent(in) :: nlconstr(:) -logical, intent(out) :: terminate +real(RP), intent(in), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: terminate ! Local variables integer(C_INT) :: n_sub_loc @@ -124,9 +131,9 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) real(C_DOUBLE) :: f_sub_loc integer(C_INT) :: nf_sub_loc integer(C_INT) :: tr_loc -real(C_DOUBLE) :: cstrv_loc +real(C_DOUBLE) :: cstrv_sub_loc integer(C_INT) :: m_nlconstr -real(C_DOUBLE) :: nlconstr_loc(size(nlconstr)) +real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc ! Read the inputs and convert them to the types specified in CCALLBACK @@ -135,18 +142,36 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv, nlconstr, terminate) f_sub_loc = real(f_sub, kind(f_sub_loc)) nf_sub_loc = int(nf_sub, kind(nf_sub_loc)) tr_loc = int(tr, kind(tr_loc)) -cstrv_loc = real(cstrv, kind(cstrv_loc)) -m_nlconstr = size(nlconstr) -nlconstr_loc = real(nlconstr, kind(nlconstr_loc)) + +! Set the constraint violation to a sensible default value if it is not provided. +if (present(cstrv_sub)) then + cstrv_sub_loc = real(cstrv_sub, kind(cstrv_sub_loc)) +else + cstrv_sub_loc = 0.0_C_DOUBLE +end if + +! Set the nonlinear constraints to a sensible default value if it is not provided. +if (present(nlconstr_sub)) then + m_nlconstr = int(size(nlconstr_sub), C_INT) + call safealloc(nlconstr_sub_loc, int(m_nlconstr, IK)) + nlconstr_sub_loc = real(nlconstr_sub, kind(nlconstr_sub_loc)) +else + m_nlconstr = 0_C_INT + nlconstr_sub_loc = [real(C_DOUBLE) ::] +end if ! Call the C objective function -call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_loc, m_nlconstr, nlconstr_loc, terminate_loc) +call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output -terminate = logical(terminate_loc, kind(terminate)) +if ( present(terminate) ) then + terminate = logical(terminate_loc, kind(terminate)) +end if -end subroutine callback_fcn +! Deallocate resources +if (allocated(nlconstr_sub_loc)) deallocate(nlconstr_sub_loc) +end subroutine callback_fcn end subroutine uobyqa_c diff --git a/fortran/bobyqa/bobyqb.f90 b/fortran/bobyqa/bobyqb.f90 index 87c11e08ab..cc7f12e64a 100644 --- a/fortran/bobyqa/bobyqb.f90 +++ b/fortran/bobyqa/bobyqb.f90 @@ -221,10 +221,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Report the current best value, and check if user asks for early termination. terminate = .false. if (present(callback_fcn)) then - ! Since this algorithm is a feasible method, the constraint violation is always 0. - ! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set - ! to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then subinfo = USER_STOP end if @@ -588,10 +585,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Report the current best value, and check if user asks for early termination. if (present(callback_fcn)) then - ! Since this algorithm is a feasible method, the constraint violation is always 0. - ! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set - ! to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then info = USER_STOP exit diff --git a/fortran/common/pintrf.f90 b/fortran/common/pintrf.f90 index 3674f6dfdc..cbc1f1490c 100644 --- a/fortran/common/pintrf.f90 +++ b/fortran/common/pintrf.f90 @@ -43,9 +43,9 @@ subroutine CALLBACK(x, f, nf, tr, cstrv, nlconstr, terminate) real(RP), intent(in) :: f integer(IK), intent(in) :: nf integer(IK), intent(in) :: tr - real(RP), intent(in) :: cstrv - real(RP), intent(in) :: nlconstr(:) - logical, intent(out) :: terminate + real(RP), intent(in), optional :: cstrv + real(RP), intent(in), optional :: nlconstr(:) + logical, intent(out), optional :: terminate end subroutine CALLBACK end interface diff --git a/fortran/lincoa/lincob.f90 b/fortran/lincoa/lincob.f90 index 446ed4892c..f93a59c2b9 100644 --- a/fortran/lincoa/lincob.f90 +++ b/fortran/lincoa/lincob.f90 @@ -258,9 +258,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b ! Report the current best value, and check if user asks for early termination. terminate = .false. if (present(callback_fcn)) then - ! Since we only solve linearly constrained problems, the nonlinear constraint is set - ! to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), terminate=terminate) if (terminate) then subinfo = USER_STOP end if @@ -651,9 +649,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b ! Report the current best value, and check if user asks for early termination. if (present(callback_fcn)) then - ! Since we only solve linearly constrained problems, the nonlinear constraint is set - ! to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), terminate=terminate) if (terminate) then info = USER_STOP exit diff --git a/fortran/newuoa/newuob.f90 b/fortran/newuoa/newuob.f90 index 1e67590fc8..93fbd79df1 100644 --- a/fortran/newuoa/newuob.f90 +++ b/fortran/newuoa/newuob.f90 @@ -180,9 +180,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Report the current best value, and check if user asks for early termination. terminate = .false. if (present(callback_fcn)) then - ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear - ! constraint is set to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then subinfo = USER_STOP end if @@ -590,9 +588,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Report the current best value, and check if user asks for early termination. if (present(callback_fcn)) then - ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear - ! constraint is set to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then info = USER_STOP exit diff --git a/fortran/primaf-GNU.def b/fortran/primaf-GNU.def index 9b3eb465fc..b59f09635d 100644 --- a/fortran/primaf-GNU.def +++ b/fortran/primaf-GNU.def @@ -1,4 +1,6 @@ EXPORTS + __memory_mod_MOD_alloc_rvector_sp + __memory_mod_MOD_alloc_rvector_dp __bobyqa_mod_MOD_bobyqa __cobyla_mod_MOD_cobyla __lincoa_mod_MOD_lincoa diff --git a/fortran/primaf-Intel.def b/fortran/primaf-Intel.def index 936abb5af4..023ae973a6 100644 --- a/fortran/primaf-Intel.def +++ b/fortran/primaf-Intel.def @@ -1,4 +1,6 @@ EXPORTS + MEMORY_MOD_mp_ALLOC_RVECTOR_SP + MEMORY_MOD_mp_ALLOC_RVECTOR_DP BOBYQA_MOD_mp_BOBYQA COBYLA_MOD_mp_COBYLA LINCOA_MOD_mp_LINCOA diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90 index f854664707..c0ef0b3fc2 100644 --- a/fortran/uobyqa/uobyqb.f90 +++ b/fortran/uobyqa/uobyqb.f90 @@ -171,9 +171,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r ! Report the current best value, and check if user asks for early termination. terminate = .false. if (present(callback_fcn)) then - ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear - ! constraint is set to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then subinfo = USER_STOP end if @@ -475,9 +473,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r ! Report the current best value, and check if user asks for early termination. if (present(callback_fcn)) then - ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear - ! constraint is set to an empty array. - call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate) + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then info = USER_STOP exit From b2efd7bda186baf78902330a30b914cc71bf410b Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Fri, 15 Dec 2023 17:59:47 +0800 Subject: [PATCH 04/13] Apparently the intel compiler struggles... It struggles with capturing a procedure pointer inside a closure, although it appears capable of capturing a C function pointer. As a result we move the C_F_PROCPOINTER call back inside the closure. --- c/bobyqa_c.f90 | 27 +++++++++++++++++---------- c/cobyla_c.f90 | 27 +++++++++++++++++---------- c/lincoa_c.f90 | 29 ++++++++++++++++++----------- c/newuoa_c.f90 | 27 +++++++++++++++++---------- c/uobyqa_c.f90 | 30 +++++++++++++++++++----------- 5 files changed, 88 insertions(+), 52 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index be5d22b708..54eb1cbe03 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -15,8 +15,7 @@ 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 -use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: bobyqa_mod, only : bobyqa implicit none @@ -53,8 +52,6 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -procedure(CCALLBACK), pointer :: cb_ptr -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -66,13 +63,11 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture it for use in the closure below - call C_F_PROCPOINTER(callback_ptr, cb_ptr) - ! And then we pass the closure to the Fortran code + ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, + ! and then we pass the closure to the Fortran code. call bobyqa(calfun, x_loc, f_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, callback_fcn=callback_fcn, info=info_loc) else @@ -95,8 +90,9 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : COBJ use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub @@ -104,10 +100,15 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) +! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! here as opposed to within the parent function, otherwise it segfaults. +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) + ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -124,9 +125,10 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub @@ -146,6 +148,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -171,6 +174,10 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if +! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! within the parent function, once. +call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! Call the C objective function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 2ba1fd59ce..9dbbe6751d 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -14,8 +14,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) -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 : COBJCON, CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: cobyla_mod, only : cobyla implicit none @@ -70,8 +69,6 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP) :: xu_loc(n) real(RP) :: f0_loc real(RP) :: nlconstr0_loc(m_nlcon) -procedure(CCALLBACK), pointer :: cb_ptr -procedure(COBJCON), pointer :: objcon_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -91,13 +88,11 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) m_nlcon_loc = int(m_nlcon, kind(m_nlcon_loc)) -call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture it for use in the closure below - call C_F_PROCPOINTER(callback_ptr, cb_ptr) - ! And then we pass the closure to the Fortran code + ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, + ! and then we pass the closure to the Fortran code. 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, & @@ -128,8 +123,9 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calcfc(x_sub, f_sub, constr_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : COBJCON use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub @@ -139,10 +135,15 @@ subroutine calcfc(x_sub, f_sub, constr_sub) real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc real(C_DOUBLE) :: constr_sub_loc(size(constr_sub)) +procedure(COBJCON), pointer :: objcon_ptr ! Read the inputs and convert them to the types specified in COBJCON x_sub_loc = real(x_sub, kind(x_sub_loc)) +! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! here as opposed to within the parent function, otherwise it segfaults. +call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) + ! Call the C objective function call objcon_ptr(x_sub_loc, f_sub_loc, constr_sub_loc, data_ptr) @@ -160,9 +161,10 @@ end subroutine calcfc ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub @@ -182,6 +184,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -207,6 +210,10 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if +! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! within the parent function, once. +call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! Call the C objective function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index 0ddb65a4d3..79d7669d80 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -15,8 +15,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) -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, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: lincoa_mod, only : lincoa implicit none @@ -65,8 +64,6 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -procedure(CCALLBACK), pointer :: cb_ptr -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -84,13 +81,11 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture it for use in the closure below - call C_F_PROCPOINTER(callback_ptr, cb_ptr) - ! And then we pass the closure to the Fortran code + ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, + ! and then we pass the closure to the Fortran code. 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, & @@ -122,19 +117,25 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! A possible security downside is that the compiler must allow for an executable stack. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : COBJ use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) +! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! here as opposed to within the parent function, otherwise it segfaults. +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) + ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -151,9 +152,10 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub @@ -173,6 +175,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -198,6 +201,10 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if +! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! within the parent function, once. +call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! Call the C objective function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index b079a03fb0..f1c80dfb8b 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -13,8 +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, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: newuoa_mod, only : newuoa implicit none @@ -47,8 +46,6 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) -procedure(CCALLBACK), pointer :: cb_ptr -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -58,13 +55,11 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture it for use in the closure below - call C_F_PROCPOINTER(callback_ptr, cb_ptr) - ! And then we pass the closure to the Fortran code + ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, + ! and then we pass the closure to the Fortran code. call newuoa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) else @@ -87,8 +82,9 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : COBJ use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub @@ -96,10 +92,15 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) +! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! here as opposed to within the parent function, otherwise it segfaults. +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) + ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -116,9 +117,10 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub @@ -138,6 +140,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -163,6 +166,10 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if +! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! within the parent function, once. +call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! Call the C objective function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index 9c744b0124..ead0e6c804 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -13,8 +13,7 @@ module uobyqa_c_mod subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, 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 -use, non_intrinsic :: cintrf_mod, only : COBJ, CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: uobyqa_mod, only : uobyqa implicit none @@ -46,8 +45,6 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) -procedure(CCALLBACK), pointer :: cb_ptr -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -56,13 +53,11 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ftarget_loc = real(ftarget, kind(ftarget_loc)) maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture it for use in the closure below - call C_F_PROCPOINTER(callback_ptr, cb_ptr) - ! And then we pass the closure to the Fortran code + ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, + ! and then we pass the closure to the Fortran code. call uobyqa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & & maxfun=maxfun_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) else @@ -82,21 +77,28 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : COBJ use, non_intrinsic :: consts_mod, only : RP -use, intrinsic :: iso_c_binding, only : C_DOUBLE implicit none -real(RP), intent(in) :: x_sub(:) +real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(out) :: f_sub ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc +procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) +! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! here as opposed to within the parent function, otherwise it segfaults. +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) + ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -113,9 +115,10 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER +use, non_intrinsic :: cintrf_mod, only : CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables real(RP), intent(in) :: f_sub @@ -135,6 +138,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -160,6 +164,10 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if +! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! within the parent function, once. +call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! Call the C objective function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) From b6930e5f46a13948d807fb526775324ec9240d53 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Fri, 15 Dec 2023 19:26:45 +0800 Subject: [PATCH 05/13] And now the other intel compiler... Now that we've moved the procedure pointer into the closure so that the new intel compiler, ifx, works, the old intel compiler, ifort, has started to have issues with that data_ptr that gets passed around. The issue was observed while running data_c_exe. It passes a data pointer and in the callback it checks to make sure that it gets back the same data pointer. It started complaining. Removing the value attribute seems to have fixed it, and this seems to work for other compilers as well. --- c/bobyqa_c.f90 | 2 +- c/cintrf.f90 | 4 ++-- c/cobyla_c.f90 | 2 +- c/lincoa_c.f90 | 2 +- c/newuoa_c.f90 | 2 +- c/uobyqa_c.f90 | 3 +-- 6 files changed, 7 insertions(+), 8 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 54eb1cbe03..fd0fc3980f 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -22,7 +22,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in), value :: data_ptr +type(C_PTR), intent(in) :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) diff --git a/c/cintrf.f90 b/c/cintrf.f90 index e66ef17178..7d972525d1 100644 --- a/c/cintrf.f90 +++ b/c/cintrf.f90 @@ -18,7 +18,7 @@ subroutine COBJ(x, f, data_ptr) bind(c) ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(in) :: x(*) real(C_DOUBLE), intent(out) :: f - type(C_PTR), intent(in), value :: data_ptr + type(C_PTR), intent(in) :: data_ptr end subroutine COBJ subroutine COBJCON(x, f, constr, data_ptr) bind(c) @@ -28,7 +28,7 @@ subroutine COBJCON(x, f, constr, data_ptr) bind(c) real(C_DOUBLE), intent(in) :: x(*) real(C_DOUBLE), intent(out) :: f real(C_DOUBLE), intent(out) :: constr(*) - type(C_PTR), intent(in), value :: data_ptr + type(C_PTR), intent(in) :: data_ptr end subroutine COBJCON subroutine CCALLBACK(n, x, f, nf, tr, cstrv, m_nlcon, nlconstr, terminate) bind(c) diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 9dbbe6751d..b4f9f08b36 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -22,7 +22,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! Compulsory arguments integer(C_INT), intent(in), value :: m_nlcon type(C_FUNPTR), intent(in), value :: cobjcon_ptr -type(C_PTR), intent(in), value :: data_ptr +type(C_PTR), intent(in) :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index 79d7669d80..b68430f4b2 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -22,7 +22,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in), value :: data_ptr +type(C_PTR), intent(in) :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index f1c80dfb8b..9f8a6f930a 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -20,7 +20,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in), value :: data_ptr +type(C_PTR), intent(in) :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index ead0e6c804..692bd57ce7 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -20,8 +20,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in), value :: data_ptr - +type(C_PTR), intent(in) :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) From aa5e57723d085c29c8da346cbd80a4f7019fbf1b Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Mon, 18 Dec 2023 22:19:00 +0800 Subject: [PATCH 06/13] Reverting the change to the value attribute of type(C_PTR) and other cosmetic fixes After reviewing the Fortran 2008, 2018, and 2023 standards, it looks like value is meant to be an attribute for type(C_PTR) in the way that we are using it, so it is probably best to leave it in and accept the fact that the classic Intel fortran compiler does not work in this case. It's worth noting that the only case where there is an issue is when sending custom data to the objective function, which is not a common use case, and hopefully those who need that functionality are able to use a compiler other than the classic Intel one. Minimal example: https://fortran-lang.discourse.group/t/problem-with-ifort-compiler-and-type-c-ptr/7011 Other various minor changes are made based on feedback from an in-person review session with @zaikunzhang --- c/bobyqa_c.f90 | 9 +++++---- c/cintrf.f90 | 4 ++-- c/cobyla_c.f90 | 23 +++++++++++------------ c/examples/bobyqa/bobyqa_example.c | 2 +- c/include/prima/prima.h | 6 +++--- c/lincoa_c.f90 | 28 +++++++++++++--------------- c/newuoa_c.f90 | 9 +++++---- c/prima.c | 4 ++-- c/uobyqa_c.f90 | 9 +++++---- fortran/bobyqa/bobyqa.f90 | 2 +- fortran/bobyqa/bobyqb.f90 | 6 +++--- fortran/cobyla/cobyla.f90 | 2 +- fortran/cobyla/cobylb.f90 | 6 +++--- fortran/common/infos.f90 | 4 ++-- fortran/common/message.f90 | 6 +++--- fortran/lincoa/lincoa.f90 | 2 +- fortran/lincoa/lincob.f90 | 6 +++--- fortran/newuoa/newuoa.f90 | 2 +- fortran/newuoa/newuob.f90 | 6 +++--- fortran/uobyqa/uobyqa.f90 | 2 +- fortran/uobyqa/uobyqb.f90 | 6 +++--- 21 files changed, 72 insertions(+), 72 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index fd0fc3980f..fc0318ec7e 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -22,7 +22,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in) :: data_ptr +type(C_PTR), intent(in), value :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -87,7 +87,8 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. -! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. +! This subroutine is identical across 4 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER @@ -105,7 +106,7 @@ subroutine calfun(x_sub, f_sub) ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer ! here as opposed to within the parent function, otherwise it segfaults. call C_F_PROCPOINTER(cobj_ptr, obj_ptr) @@ -174,7 +175,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to ! within the parent function, once. call C_F_PROCPOINTER(callback_ptr, cb_ptr) diff --git a/c/cintrf.f90 b/c/cintrf.f90 index 7d972525d1..e66ef17178 100644 --- a/c/cintrf.f90 +++ b/c/cintrf.f90 @@ -18,7 +18,7 @@ subroutine COBJ(x, f, data_ptr) bind(c) ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(in) :: x(*) real(C_DOUBLE), intent(out) :: f - type(C_PTR), intent(in) :: data_ptr + type(C_PTR), intent(in), value :: data_ptr end subroutine COBJ subroutine COBJCON(x, f, constr, data_ptr) bind(c) @@ -28,7 +28,7 @@ subroutine COBJCON(x, f, constr, data_ptr) bind(c) real(C_DOUBLE), intent(in) :: x(*) real(C_DOUBLE), intent(out) :: f real(C_DOUBLE), intent(out) :: constr(*) - type(C_PTR), intent(in) :: data_ptr + type(C_PTR), intent(in), value :: data_ptr end subroutine COBJCON subroutine CCALLBACK(n, x, f, nf, tr, cstrv, m_nlcon, nlconstr, terminate) bind(c) diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index b4f9f08b36..33f43a2431 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -22,7 +22,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! Compulsory arguments integer(C_INT), intent(in), value :: m_nlcon type(C_FUNPTR), intent(in), value :: cobjcon_ptr -type(C_PTR), intent(in) :: data_ptr +type(C_PTR), intent(in), value :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -93,16 +93,14 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ if (C_ASSOCIATED(callback_ptr)) then ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, ! and then we pass the closure to the Fortran code. - 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, & + 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, & & 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, & + 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, & & iprint=iprint_loc, info=info_loc) end if @@ -120,7 +118,8 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! This subroutine defines `calcfc` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. -! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. +! This subroutine is identical across 4 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calcfc(x_sub, f_sub, constr_sub) use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER @@ -140,7 +139,7 @@ subroutine calcfc(x_sub, f_sub, constr_sub) ! Read the inputs and convert them to the types specified in COBJCON x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer ! here as opposed to within the parent function, otherwise it segfaults. call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) @@ -210,7 +209,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to ! within the parent function, once. call C_F_PROCPOINTER(callback_ptr, cb_ptr) diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index 1b8cb617d0..cd8adfae32 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -16,7 +16,7 @@ static void callback(int n, const double x[], double f, int nf, int tr, double c { (void)n; printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); - *terminate = x[1] > 1.8; // test early termination + *terminate = 0; (void)m_nlcon; (void)nlconstr; } diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h index babc335b2f..e2d63df0ea 100644 --- a/c/include/prima/prima.h +++ b/c/include/prima/prima.h @@ -51,7 +51,7 @@ typedef enum PRIMA_NO_SPACE_BETWEEN_BOUNDS = 6, PRIMA_DAMAGING_ROUNDING = 7, PRIMA_ZERO_LINEAR_CONSTRAINT = 8, - PRIMA_USER_STOP = 9, + PRIMA_CALLBACK_TERMINATE = 30, PRIMA_INVALID_INPUT = 100, PRIMA_ASSERTION_FAILS = 101, PRIMA_VALIDATION_FAILS = 102, @@ -75,7 +75,7 @@ const char *prima_get_rc_string(const prima_rc_t rc); * x : on input, then vector of variables (should not be modified) * f : on output, the value of the function * a NaN value can be passed to signal an evaluation error - * data : user-data + * data : user data * constr : on output, the value of the constraints (of size m_nlcon) * NaN values can be passed to signal evaluation errors * only for cobyla @@ -120,7 +120,7 @@ typedef struct { // ignored for uobyqa & cobyla int npt; - // user-data, will be passed through the objective function callback + // user data, will be passed through the objective function callback void *data; // callback function to report algorithm progress (default=NULL) diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index b68430f4b2..d508f6e298 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -22,7 +22,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in) :: data_ptr +type(C_PTR), intent(in), value :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -86,19 +86,15 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ if (C_ASSOCIATED(callback_ptr)) then ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, ! and then we pass the closure to the Fortran code. - 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, callback_fcn=callback_fcn, info=info_loc) + 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, & + & 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, info=info_loc) + 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, & + & info=info_loc) end if @@ -115,6 +111,8 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. +! This subroutine is identical across 4 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER @@ -132,7 +130,7 @@ subroutine calfun(x_sub, f_sub) ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer ! here as opposed to within the parent function, otherwise it segfaults. call C_F_PROCPOINTER(cobj_ptr, obj_ptr) @@ -201,7 +199,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to ! within the parent function, once. call C_F_PROCPOINTER(callback_ptr, cb_ptr) diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index 9f8a6f930a..da7e7c9d4a 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -20,7 +20,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in) :: data_ptr +type(C_PTR), intent(in), value :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -79,7 +79,8 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. -! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. +! This subroutine is identical across 4 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER @@ -97,7 +98,7 @@ subroutine calfun(x_sub, f_sub) ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer ! here as opposed to within the parent function, otherwise it segfaults. call C_F_PROCPOINTER(cobj_ptr, obj_ptr) @@ -166,7 +167,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to ! within the parent function, once. call C_F_PROCPOINTER(callback_ptr, cb_ptr) diff --git a/c/prima.c b/c/prima.c index a7ecfa75b2..c250c6d7ae 100644 --- a/c/prima.c +++ b/c/prima.c @@ -261,8 +261,8 @@ const char *prima_get_rc_string(const prima_rc_t rc) return "Rounding errors are becoming damaging"; case PRIMA_ZERO_LINEAR_CONSTRAINT: return "One of the linear constraints has a zero gradient"; - case PRIMA_USER_STOP: - return "User requested end of optimization"; + case PRIMA_CALLBACK_TERMINATE: + return "Callback function requested early termination of optimization"; case PRIMA_INVALID_INPUT: return "Invalid input"; case PRIMA_ASSERTION_FAILS: diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index 692bd57ce7..e46b2245ef 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -20,7 +20,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! Compulsory arguments type(C_FUNPTR), intent(IN), value :: cobj_ptr -type(C_PTR), intent(in) :: data_ptr +type(C_PTR), intent(in), value :: data_ptr integer(C_INT), intent(in), value :: n ! We cannot use assumed-shape arrays for C interoperability real(C_DOUBLE), intent(inout) :: x(n) @@ -76,7 +76,8 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! This subroutine defines `calfun` using the C function pointer with an internal subroutine. ! This allows to avoid passing the C function pointer by a module variable, which is thread-unsafe. ! A possible security downside is that the compiler must allow for an executable stack. -! This subroutine is identical across 4/5 algorithms; COBYLA requires a slightly different signature. +! This subroutine is identical across 4 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER @@ -94,7 +95,7 @@ subroutine calfun(x_sub, f_sub) ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The intel compiler insists that we convert the C function pointer to a Fortran function pointer +! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer ! here as opposed to within the parent function, otherwise it segfaults. call C_F_PROCPOINTER(cobj_ptr, obj_ptr) @@ -163,7 +164,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the intel compiler insists on doing this conversion here, every time, as opposed to +! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to ! within the parent function, once. call C_F_PROCPOINTER(callback_ptr, cb_ptr) diff --git a/fortran/bobyqa/bobyqa.f90 b/fortran/bobyqa/bobyqa.f90 index 6662834203..8092eae92d 100644 --- a/fortran/bobyqa/bobyqa.f90 +++ b/fortran/bobyqa/bobyqa.f90 @@ -169,7 +169,7 @@ subroutine bobyqa(calfun, x, & ! X0 if needed. See the PREPROC subroutine for more information. ! ! CALLBACK -! Input, function to report progress and request termination. +! Input, function to report progress and optionally request termination. ! ! INFO ! Output, INTEGER(IK) scalar. diff --git a/fortran/bobyqa/bobyqb.f90 b/fortran/bobyqa/bobyqb.f90 index cc7f12e64a..2a3a8aabac 100644 --- a/fortran/bobyqa/bobyqb.f90 +++ b/fortran/bobyqa/bobyqb.f90 @@ -82,7 +82,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP +use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, CALLBACK_TERMINATE use, non_intrinsic :: linalg_mod, only : norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK @@ -223,7 +223,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then - subinfo = USER_STOP + subinfo = CALLBACK_TERMINATE end if endif @@ -587,7 +587,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then - info = USER_STOP + info = CALLBACK_TERMINATE exit end if end if diff --git a/fortran/cobyla/cobyla.f90 b/fortran/cobyla/cobyla.f90 index acbd81e6f7..bde122dc48 100644 --- a/fortran/cobyla/cobyla.f90 +++ b/fortran/cobyla/cobyla.f90 @@ -235,7 +235,7 @@ subroutine cobyla(calcfc, m_nlcon, x, & ! see common/consts.F90 for the definitions of MAXFILT_DFT and MIN_MAXFILT. ! ! CALLBACK -! Input, function to report progress and request termination. +! Input, function to report progress and optionally request termination. ! ! INFO ! Output, INTEGER(IK) scalar. diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index 2f21b053af..2aca8b1365 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -46,7 +46,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING, USER_STOP +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING, CALLBACK_TERMINATE use, non_intrinsic :: linalg_mod, only : inprod, matprod, norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg use, non_intrinsic :: pintrf_mod, only : OBJCON, CALLBACK @@ -220,7 +220,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et if (present(callback_fcn)) then call callback_fcn(sim(:, n+1), fval(n+1), nf, 0, cval(n+1), conmat(:, n+1), terminate) if (terminate) then - subinfo = USER_STOP + subinfo = CALLBACK_TERMINATE end if end if @@ -607,7 +607,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et if (present(callback_fcn)) then call callback_fcn(sim(:, n+1), fval(n+1), nf, tr, cval(n+1), conmat(:, n+1), terminate) if (terminate) then - info = USER_STOP + info = CALLBACK_TERMINATE exit end if end if diff --git a/fortran/common/infos.f90 b/fortran/common/infos.f90 index 95344a4664..e474a5b98b 100644 --- a/fortran/common/infos.f90 +++ b/fortran/common/infos.f90 @@ -25,7 +25,7 @@ module infos_mod public :: DAMAGING_ROUNDING public :: NO_SPACE_BETWEEN_BOUNDS public :: ZERO_LINEAR_CONSTRAINT -public :: USER_STOP +public :: CALLBACK_TERMINATE public :: INVALID_INPUT public :: ASSERTION_FAILS public :: VALIDATION_FAILS @@ -43,7 +43,7 @@ module infos_mod integer(IK), parameter :: NO_SPACE_BETWEEN_BOUNDS = 6 integer(IK), parameter :: DAMAGING_ROUNDING = 7 integer(IK), parameter :: ZERO_LINEAR_CONSTRAINT = 8 -integer(IK), parameter :: USER_STOP = 9 +integer(IK), parameter :: CALLBACK_TERMINATE = 30 ! Stop-codes. ! The following codes are used by ERROR STOP as stop-codes, which should be default integers. diff --git a/fortran/common/message.f90 b/fortran/common/message.f90 index c254a3a005..5f7b99d226 100644 --- a/fortran/common/message.f90 +++ b/fortran/common/message.f90 @@ -33,7 +33,7 @@ subroutine retmsg(solver, info, iprint, nf, f, x, cstrv, constr) use, non_intrinsic :: fprint_mod, only : fprint use, non_intrinsic :: infos_mod, only : FTARGET_ACHIEVED, MAXFUN_REACHED, MAXTR_REACHED, & & SMALL_TR_RADIUS, TRSUBP_FAILED, NAN_INF_X, NAN_INF_F, NAN_INF_MODEL, DAMAGING_ROUNDING, & - & NO_SPACE_BETWEEN_BOUNDS, ZERO_LINEAR_CONSTRAINT, USER_STOP + & NO_SPACE_BETWEEN_BOUNDS, ZERO_LINEAR_CONSTRAINT, CALLBACK_TERMINATE use, non_intrinsic :: string_mod, only : strip, num2str implicit none @@ -124,8 +124,8 @@ subroutine retmsg(solver, info, iprint, nf, f, x, cstrv, constr) reason = 'there is no space between the lower and upper bounds of variable.' case (ZERO_LINEAR_CONSTRAINT) reason = 'one of the linear constraints has a zero gradient' -case (USER_STOP) - reason = 'user requested end of optimization' +case (CALLBACK_TERMINATE) + reason = 'callback function requested early termination of optimization' case default reason = 'UNKNOWN EXIT FLAG' end select diff --git a/fortran/lincoa/lincoa.f90 b/fortran/lincoa/lincoa.f90 index e5a1574a29..8cd7a54f8b 100644 --- a/fortran/lincoa/lincoa.f90 +++ b/fortran/lincoa/lincoa.f90 @@ -196,7 +196,7 @@ subroutine lincoa(calfun, x, & ! Use *HIST with caution! (N.B.: the algorithm is NOT designed for large problems). ! ! CALLBACK -! Input, function to report progress and request termination. +! Input, function to report progress and optionally request termination. ! ! INFO ! Output, INTEGER(IK) scalar. diff --git a/fortran/lincoa/lincob.f90 b/fortran/lincoa/lincob.f90 index f93a59c2b9..35ce812e44 100644 --- a/fortran/lincoa/lincob.f90 +++ b/fortran/lincoa/lincob.f90 @@ -78,7 +78,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, CALLBACK_TERMINATE use, non_intrinsic :: linalg_mod, only : matprod, maximum, eye, trueloc, linspace, norm, trueloc use, non_intrinsic :: memory_mod, only : safealloc use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg @@ -260,7 +260,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), terminate=terminate) if (terminate) then - subinfo = USER_STOP + subinfo = CALLBACK_TERMINATE end if end if @@ -651,7 +651,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), terminate=terminate) if (terminate) then - info = USER_STOP + info = CALLBACK_TERMINATE exit end if end if diff --git a/fortran/newuoa/newuoa.f90 b/fortran/newuoa/newuoa.f90 index f53222310e..05318d4194 100644 --- a/fortran/newuoa/newuoa.f90 +++ b/fortran/newuoa/newuoa.f90 @@ -143,7 +143,7 @@ subroutine newuoa(calfun, x, & ! Use *HIST with caution! (N.B.: the algorithm is NOT designed for large problems). ! ! CALLBACK -! Input, function to report progress and request termination. +! Input, function to report progress and optionally request termination. ! ! INFO ! Output, INTEGER(IK) scalar. diff --git a/fortran/newuoa/newuob.f90 b/fortran/newuoa/newuob.f90 index 93fbd79df1..8b948f287f 100644 --- a/fortran/newuoa/newuob.f90 +++ b/fortran/newuoa/newuob.f90 @@ -57,7 +57,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP +use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, CALLBACK_TERMINATE use, non_intrinsic :: linalg_mod, only : norm use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK @@ -182,7 +182,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then - subinfo = USER_STOP + subinfo = CALLBACK_TERMINATE end if end if @@ -590,7 +590,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then - info = USER_STOP + info = CALLBACK_TERMINATE exit end if end if diff --git a/fortran/uobyqa/uobyqa.f90 b/fortran/uobyqa/uobyqa.f90 index 297ec04673..546c8d84d6 100644 --- a/fortran/uobyqa/uobyqa.f90 +++ b/fortran/uobyqa/uobyqa.f90 @@ -135,7 +135,7 @@ subroutine uobyqa(calfun, x, & ! Use *HIST with caution!!! (N.B.: the algorithm is NOT designed for large problems). ! ! CALLBACK -! Input, function to report progress and request termination. +! Input, function to report progress and optionally request termination. ! ! INFO ! Output, INTEGER(IK) scalar. diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90 index c0ef0b3fc2..244782ab5e 100644 --- a/fortran/uobyqa/uobyqb.f90 +++ b/fortran/uobyqa/uobyqb.f90 @@ -50,7 +50,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r use, non_intrinsic :: evaluate_mod, only : evaluate use, non_intrinsic :: history_mod, only : savehist, rangehist use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf -use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP +use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, CALLBACK_TERMINATE use, non_intrinsic :: linalg_mod, only : vec2smat, smat_mul_vec, norm use, non_intrinsic :: memory_mod, only : safealloc use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg @@ -173,7 +173,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) if (terminate) then - subinfo = USER_STOP + subinfo = CALLBACK_TERMINATE end if end if @@ -475,7 +475,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r if (present(callback_fcn)) then call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) if (terminate) then - info = USER_STOP + info = CALLBACK_TERMINATE exit end if end if From 86b1fe9da102ceda6cd1e863719c53150381d463 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Tue, 19 Dec 2023 22:53:50 +0800 Subject: [PATCH 07/13] One more attempt to get things working for all compilers --- c/CMakeLists.txt | 14 ++++++++++++++ c/bobyqa_c.f90 | 30 ++++++++++++------------------ c/cobyla_c.f90 | 30 ++++++++++++------------------ c/lincoa_c.f90 | 30 ++++++++++++------------------ c/newuoa_c.f90 | 30 ++++++++++++------------------ c/uobyqa_c.f90 | 30 ++++++++++++------------------ 6 files changed, 74 insertions(+), 90 deletions(-) diff --git a/c/CMakeLists.txt b/c/CMakeLists.txt index 4440796e3c..144118cc7d 100644 --- a/c/CMakeLists.txt +++ b/c/CMakeLists.txt @@ -6,6 +6,20 @@ if (WIN32) set_target_properties(primac PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin) endif() +# The newer intel compiler has issues when trying to capture a Fortran procedure pointer in a closure. +# We can fix this by appending norecursion to our compile options when we compile the C interface. +# More details may be found here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +if (CMAKE_Fortran_COMPILER_ID MATCHES "IntelLLVM") + if (WIN32) + # The $<$ part is so that this only applies when compiling Fortran + # files and *not* when compiling prima.c (because 'assume' is not a valid compiler flag when + # compiling C). + target_compile_options(primac PRIVATE $<$:/assume:norecursion>) + else () + target_compile_options(primac PRIVATE $<$:-assume norecursion>) + endif () +endif () + target_include_directories (primac PUBLIC $ $ diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index fc0318ec7e..23a769c187 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -15,7 +15,8 @@ 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 +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 :: bobyqa_mod, only : bobyqa implicit none @@ -52,6 +53,8 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) +procedure(COBJ), pointer :: obj_ptr +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -63,11 +66,14 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, - ! and then we pass the closure to the Fortran code. + ! If a C callback function is provided, we convert it to a Fortran procedure pointer and capture + ! that pointer in the closure below. + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! We then provide the closure to the algorithm. call bobyqa(calfun, x_loc, f_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, callback_fcn=callback_fcn, info=info_loc) else @@ -91,8 +97,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & ! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : COBJ +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables @@ -101,15 +106,10 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer -! here as opposed to within the parent function, otherwise it segfaults. -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) - ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -126,8 +126,7 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc implicit none @@ -149,7 +148,6 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc -procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -175,11 +173,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to -! within the parent function, once. -call C_F_PROCPOINTER(callback_ptr, cb_ptr) - -! Call the C objective function +! Call the C callback function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 33f43a2431..7f70e8bcbe 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -14,7 +14,8 @@ 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) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED +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 : COBJCON, CCALLBACK use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: cobyla_mod, only : cobyla implicit none @@ -69,6 +70,8 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP) :: xu_loc(n) real(RP) :: f0_loc real(RP) :: nlconstr0_loc(m_nlcon) +procedure(COBJCON), pointer :: objcon_ptr +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -88,11 +91,14 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) m_nlcon_loc = int(m_nlcon, kind(m_nlcon_loc)) +call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, - ! and then we pass the closure to the Fortran code. + ! If a C callback function is provided, we convert it to a Fortran procedure pointer and capture + ! that pointer in the closure below. + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! 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, & @@ -122,8 +128,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ ! signature. !--------------------------------------------------------------------------------------------------! subroutine calcfc(x_sub, f_sub, constr_sub) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : COBJCON +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables @@ -134,15 +139,10 @@ subroutine calcfc(x_sub, f_sub, constr_sub) real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc real(C_DOUBLE) :: constr_sub_loc(size(constr_sub)) -procedure(COBJCON), pointer :: objcon_ptr ! Read the inputs and convert them to the types specified in COBJCON x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer -! here as opposed to within the parent function, otherwise it segfaults. -call C_F_PROCPOINTER(cobjcon_ptr, objcon_ptr) - ! Call the C objective function call objcon_ptr(x_sub_loc, f_sub_loc, constr_sub_loc, data_ptr) @@ -160,8 +160,7 @@ end subroutine calcfc ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc implicit none @@ -183,7 +182,6 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc -procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -209,11 +207,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to -! within the parent function, once. -call C_F_PROCPOINTER(callback_ptr, cb_ptr) - -! Call the C objective function +! Call the C callback function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index d508f6e298..4e44f6f50d 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -15,7 +15,8 @@ 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) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_FUNPTR, C_PTR, C_ASSOCIATED +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 :: lincoa_mod, only : lincoa implicit none @@ -64,6 +65,8 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) +procedure(COBJ), pointer :: obj_ptr +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major @@ -81,11 +84,14 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, - ! and then we pass the closure to the Fortran code. + ! If a C callback function is provided, we convert it to a Fortran procedure pointer and capture + ! that pointer in the closure below. + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! 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, & @@ -115,8 +121,7 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ ! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : COBJ +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables @@ -125,15 +130,10 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer -! here as opposed to within the parent function, otherwise it segfaults. -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) - ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -150,8 +150,7 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc implicit none @@ -173,7 +172,6 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc -procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -199,11 +197,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to -! within the parent function, once. -call C_F_PROCPOINTER(callback_ptr, cb_ptr) - -! Call the C objective function +! Call the C callback function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index da7e7c9d4a..bcba37edfc 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -13,7 +13,8 @@ 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 +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 :: newuoa_mod, only : newuoa implicit none @@ -46,6 +47,8 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) +procedure(COBJ), pointer :: obj_ptr +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -55,11 +58,14 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma maxfun_loc = int(maxfun, kind(maxfun_loc)) npt_loc = int(npt, kind(npt_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, - ! and then we pass the closure to the Fortran code. + ! If a C callback function is provided, we convert it to a Fortran procedure pointer and capture + ! that pointer in the closure below. + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! We then provide the closure to the algorithm. call newuoa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) else @@ -83,8 +89,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : COBJ +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables @@ -93,15 +98,10 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer -! here as opposed to within the parent function, otherwise it segfaults. -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) - ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -118,8 +118,7 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc implicit none @@ -141,7 +140,6 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc -procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -167,11 +165,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to -! within the parent function, once. -call C_F_PROCPOINTER(callback_ptr, cb_ptr) - -! Call the C objective function +! Call the C callback function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index e46b2245ef..ecb7ddae91 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -13,7 +13,8 @@ module uobyqa_c_mod subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, 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 +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 :: uobyqa_mod, only : uobyqa implicit none @@ -44,6 +45,8 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) +procedure(COBJ), pointer :: obj_ptr +procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -52,11 +55,14 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ftarget_loc = real(ftarget, kind(ftarget_loc)) maxfun_loc = int(maxfun, kind(maxfun_loc)) iprint_loc = int(iprint, kind(iprint_loc)) +call C_F_PROCPOINTER(cobj_ptr, obj_ptr) ! Call the Fortran code if (C_ASSOCIATED(callback_ptr)) then - ! If a C callback function is provided, we capture the callback_ptr for use in the closure below, - ! and then we pass the closure to the Fortran code. + ! If a C callback function is provided, we convert it to a Fortran procedure pointer and capture + ! that pointer in the closure below. + call C_F_PROCPOINTER(callback_ptr, cb_ptr) + ! We then provide the closure to the algorithm. call uobyqa(calfun, x_loc, f_loc, nf=nf_loc, rhobeg=rhobeg_loc, rhoend=rhoend_loc, ftarget=ftarget_loc, & & maxfun=maxfun_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc) else @@ -80,8 +86,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma ! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : COBJ +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP implicit none real(RP), intent(in) :: x_sub(:) ! We name some variables _sub to avoid masking the parent variables @@ -90,15 +95,10 @@ subroutine calfun(x_sub, f_sub) ! Local variables real(C_DOUBLE) :: x_sub_loc(size(x_sub)) real(C_DOUBLE) :: f_sub_loc -procedure(COBJ), pointer :: obj_ptr ! Read the inputs and convert them to the types specified in COBJ x_sub_loc = real(x_sub, kind(x_sub_loc)) -! The Intel compiler ifx insists that we convert the C function pointer to a Fortran function pointer -! here as opposed to within the parent function, otherwise it segfaults. -call C_F_PROCPOINTER(cobj_ptr, obj_ptr) - ! Call the C objective function call obj_ptr(x_sub_loc, f_sub_loc, data_ptr) @@ -115,8 +115,7 @@ end subroutine calfun ! This subroutine is identical across all 5 algorithms. !--------------------------------------------------------------------------------------------------! subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, terminate) -use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL, C_F_PROCPOINTER -use, non_intrinsic :: cintrf_mod, only : CCALLBACK +use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_BOOL use, non_intrinsic :: consts_mod, only : RP, IK use, non_intrinsic :: memory_mod, only : safealloc implicit none @@ -138,7 +137,6 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi integer(C_INT) :: m_nlconstr real(C_DOUBLE), allocatable :: nlconstr_sub_loc(:) logical(C_BOOL) :: terminate_loc -procedure(CCALLBACK), pointer :: cb_ptr ! Read the inputs and convert them to the types specified in CCALLBACK n_sub_loc = size(x_sub) @@ -164,11 +162,7 @@ subroutine callback_fcn(x_sub, f_sub, nf_sub, tr, cstrv_sub, nlconstr_sub, termi nlconstr_sub_loc = [real(C_DOUBLE) ::] end if -! As above, the Intel compiler ifx insists on doing this conversion here, every time, as opposed to -! within the parent function, once. -call C_F_PROCPOINTER(callback_ptr, cb_ptr) - -! Call the C objective function +! Call the C callback function call cb_ptr(n_sub_loc, x_sub_loc, f_sub_loc, nf_sub_loc, tr_loc, cstrv_sub_loc, m_nlconstr, nlconstr_sub_loc, terminate_loc) ! Write the output From 74fb158d98e93a6b5c9328bccb39588be3b6d027 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 10:18:28 +0800 Subject: [PATCH 08/13] Yet another attempt at making things work for all compilers, this one slightly less clunky --- c/CMakeLists.txt | 13 ------------- c/bobyqa_c.f90 | 6 ++++-- c/cobyla_c.f90 | 6 ++++-- c/lincoa_c.f90 | 6 ++++-- c/newuoa_c.f90 | 6 ++++-- c/uobyqa_c.f90 | 6 ++++-- 6 files changed, 20 insertions(+), 23 deletions(-) diff --git a/c/CMakeLists.txt b/c/CMakeLists.txt index 144118cc7d..960e131173 100644 --- a/c/CMakeLists.txt +++ b/c/CMakeLists.txt @@ -6,19 +6,6 @@ if (WIN32) set_target_properties(primac PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin) endif() -# The newer intel compiler has issues when trying to capture a Fortran procedure pointer in a closure. -# We can fix this by appending norecursion to our compile options when we compile the C interface. -# More details may be found here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 -if (CMAKE_Fortran_COMPILER_ID MATCHES "IntelLLVM") - if (WIN32) - # The $<$ part is so that this only applies when compiling Fortran - # files and *not* when compiling prima.c (because 'assume' is not a valid compiler flag when - # compiling C). - target_compile_options(primac PRIVATE $<$:/assume:norecursion>) - else () - target_compile_options(primac PRIVATE $<$:-assume norecursion>) - endif () -endif () target_include_directories (primac PUBLIC $ diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 23a769c187..e4a3748123 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -53,8 +53,10 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -procedure(COBJ), pointer :: obj_ptr -procedure(CCALLBACK), pointer :: cb_ptr +! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +procedure(COBJ), pointer :: obj_ptr => null() +procedure(CCALLBACK), pointer :: cb_ptr => null() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index 7f70e8bcbe..caa47bb514 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -70,8 +70,10 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP) :: xu_loc(n) real(RP) :: f0_loc real(RP) :: nlconstr0_loc(m_nlcon) -procedure(COBJCON), pointer :: objcon_ptr -procedure(CCALLBACK), pointer :: cb_ptr +! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +procedure(COBJCON), pointer :: objcon_ptr => null() +procedure(CCALLBACK), pointer :: cb_ptr => null() ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index 4e44f6f50d..d8c0e663b3 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -65,8 +65,10 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -procedure(COBJ), pointer :: obj_ptr -procedure(CCALLBACK), pointer :: cb_ptr +! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +procedure(COBJ), pointer :: obj_ptr => null() +procedure(CCALLBACK), pointer :: cb_ptr => null() ! Read the inputs and convert them to the Fortran side types ! Note that `transpose` is needed when reading 2D arrays, since they are stored in the row-major diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index bcba37edfc..a985f9d64c 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -47,8 +47,10 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) -procedure(COBJ), pointer :: obj_ptr -procedure(CCALLBACK), pointer :: cb_ptr +! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +procedure(COBJ), pointer :: obj_ptr => null() +procedure(CCALLBACK), pointer :: cb_ptr => null() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index ecb7ddae91..c87ac6380f 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -45,8 +45,10 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: rhoend_loc real(RP) :: ftarget_loc real(RP) :: x_loc(n) -procedure(COBJ), pointer :: obj_ptr -procedure(CCALLBACK), pointer :: cb_ptr +! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! See details here: https://fortran-lang.discourse.group/t/strange-issue-with-ifx-compiler-and-assume-recursion/7013 +procedure(COBJ), pointer :: obj_ptr => null() +procedure(CCALLBACK), pointer :: cb_ptr => null() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) From 65cbf5c2fd26233c7640048ce5f773c46d549a59 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 18:16:52 +0800 Subject: [PATCH 09/13] update comment --- c/bobyqa_c.f90 | 4 +++- c/cobyla_c.f90 | 4 +++- c/lincoa_c.f90 | 4 +++- c/newuoa_c.f90 | 4 +++- c/uobyqa_c.f90 | 4 +++- 5 files changed, 15 insertions(+), 5 deletions(-) diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index e4a3748123..069f351c14 100644 --- a/c/bobyqa_c.f90 +++ b/c/bobyqa_c.f90 @@ -53,8 +53,10 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, & real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! 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 +! initialization to null because it implies the 'save' attribute, which is undesirable. procedure(COBJ), pointer :: obj_ptr => null() procedure(CCALLBACK), pointer :: cb_ptr => null() diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90 index caa47bb514..961400bf1d 100644 --- a/c/cobyla_c.f90 +++ b/c/cobyla_c.f90 @@ -70,8 +70,10 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_ real(RP) :: xu_loc(n) real(RP) :: f0_loc real(RP) :: nlconstr0_loc(m_nlcon) -! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! 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 +! initialization to null because it implies the 'save' attribute, which is undesirable. procedure(COBJCON), pointer :: objcon_ptr => null() procedure(CCALLBACK), pointer :: cb_ptr => null() diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index d8c0e663b3..5e16b673c8 100644 --- a/c/lincoa_c.f90 +++ b/c/lincoa_c.f90 @@ -65,8 +65,10 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_ real(RP) :: x_loc(n) real(RP) :: xl_loc(n) real(RP) :: xu_loc(n) -! The initialization to null is necessary to avoid a bug with the newer Intel compiler. +! 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 +! initialization to null because it implies the 'save' attribute, which is undesirable. procedure(COBJ), pointer :: obj_ptr => null() procedure(CCALLBACK), pointer :: cb_ptr => null() diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index a985f9d64c..7585b55c30 100644 --- a/c/newuoa_c.f90 +++ b/c/newuoa_c.f90 @@ -47,8 +47,10 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: 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. +! 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 +! initialization to null because it implies the 'save' attribute, which is undesirable. procedure(COBJ), pointer :: obj_ptr => null() procedure(CCALLBACK), pointer :: cb_ptr => null() diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index c87ac6380f..3983b7dc2c 100644 --- a/c/uobyqa_c.f90 +++ b/c/uobyqa_c.f90 @@ -45,8 +45,10 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma real(RP) :: 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. +! 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 +! initialization to null because it implies the 'save' attribute, which is undesirable. procedure(COBJ), pointer :: obj_ptr => null() procedure(CCALLBACK), pointer :: cb_ptr => null() From 3eaf0950380ee632ddcf40b2d1a9a3c888e1d526 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 18:42:45 +0800 Subject: [PATCH 10/13] removed duplicate shell --- .github/workflows/cmake.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index f459f2042e..d672a1814f 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -122,7 +122,6 @@ jobs: cmake --build . --target install cmake --build . --target tests ctest --output-on-failure -V -E stress - shell: bash env: FC: ${{ steps.setup-fortran.outputs.fc }} shell: bash From 9ecef2f4e8abac3b88bc676e078bd37e2ebd65c7 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 18:51:39 +0800 Subject: [PATCH 11/13] remove unnecessary if statements and add comments to necessary ones --- c/CMakeLists.txt | 7 ++----- c/tests/CMakeLists.txt | 19 ++++++------------- fortran/CMakeLists.txt | 7 ++----- 3 files changed, 10 insertions(+), 23 deletions(-) diff --git a/c/CMakeLists.txt b/c/CMakeLists.txt index 960e131173..005aac3a70 100644 --- a/c/CMakeLists.txt +++ b/c/CMakeLists.txt @@ -40,17 +40,14 @@ macro (prima_add_c_test name) # Outside of CI we don't want to force people to run examples with gdb, so we test the executables by themselves. # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up # before we set up the examples for CI - if(APPLE OR UNIX) - add_test (NAME example_${name}_c COMMAND example_${name}_c_exe) - elseif(WIN32) - add_test (NAME example_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_c_exe.exe) - endif() + add_test (NAME example_${name}_c COMMAND example_${name}_c_exe) # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. # Of course this can be run locally as well if you define CI in your environment. if(NOT APPLE AND UNIX AND DEFINED ENV{CI}) # Apple security policy will not allow running gdb in CI add_test (NAME example_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb example_${name}_c_exe) elseif(WIN32 AND DEFINED ENV{CI}) + # For Windows we need to provide the full path to the executable since it is installed to a different directory add_test (NAME example_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_c_exe.exe) endif() diff --git a/c/tests/CMakeLists.txt b/c/tests/CMakeLists.txt index ed7f678116..0c9ff3d5b1 100644 --- a/c/tests/CMakeLists.txt +++ b/c/tests/CMakeLists.txt @@ -12,19 +12,11 @@ macro (prima_add_c_test_multi name) # Outside of CI we don't want to force people to run tests with gdb, so we test the executables by themselves. # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up # before we set up the tests for CI - if(APPLE OR UNIX) - add_test(NAME bobyqa_${name}_c COMMAND ${name}_c_exe bobyqa) - add_test(NAME cobyla_${name}_c COMMAND ${name}_c_exe cobyla) - add_test(NAME lincoa_${name}_c COMMAND ${name}_c_exe lincoa) - add_test(NAME newuoa_${name}_c COMMAND ${name}_c_exe newuoa) - add_test(NAME uobyqa_${name}_c COMMAND ${name}_c_exe uobyqa) - elseif(WIN32) - add_test(NAME bobyqa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe bobyqa) - add_test(NAME cobyla_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe cobyla) - add_test(NAME lincoa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe lincoa) - add_test(NAME newuoa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe newuoa) - add_test(NAME uobyqa_${name}_c COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe uobyqa) - endif() + add_test(NAME bobyqa_${name}_c COMMAND ${name}_c_exe bobyqa) + add_test(NAME cobyla_${name}_c COMMAND ${name}_c_exe cobyla) + add_test(NAME lincoa_${name}_c COMMAND ${name}_c_exe lincoa) + add_test(NAME newuoa_${name}_c COMMAND ${name}_c_exe newuoa) + add_test(NAME uobyqa_${name}_c COMMAND ${name}_c_exe uobyqa) # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. # Of course this can be run locally as well if you define CI in your environment. @@ -35,6 +27,7 @@ macro (prima_add_c_test_multi name) add_test(NAME newuoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe newuoa) add_test(NAME uobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${name}_c_exe uobyqa) elseif(WIN32 AND DEFINED ENV{CI}) + # For Windows we need to provide the full path to the executable since it is installed to a different directory add_test(NAME bobyqa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe bobyqa) add_test(NAME cobyla_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe cobyla) add_test(NAME lincoa_${name}_c_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb --args ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/${name}_c_exe.exe lincoa) diff --git a/fortran/CMakeLists.txt b/fortran/CMakeLists.txt index 3d1eb05dfe..80469a54e2 100644 --- a/fortran/CMakeLists.txt +++ b/fortran/CMakeLists.txt @@ -118,17 +118,14 @@ macro (prima_add_f_test name) # Outside of CI we don't want to force people to run examples with gdb, so we test the executables by themselves. # We want these to run in CI as well, because sometimes running with gdb masks an error, so we set them up # before we set up the examples for CI - if(APPLE OR UNIX) - add_test (NAME example_${name}_fortran COMMAND example_${name}_fortran_exe) - elseif(WIN32) - add_test (NAME example_${name}_fortran COMMAND ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_fortran_exe.exe) - endif() + add_test (NAME example_${name}_fortran COMMAND example_${name}_fortran_exe) # Within CI, we'd like to run with gdb so that if there's a segfault the logs will have a stacktrace we can use to investigate. # Of course this can be run locally as well if you define CI in your environment. if(NOT APPLE AND UNIX AND DEFINED ENV{CI}) # Apple security policy will not allow running gdb in CI add_test (NAME example_${name}_fortran_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb example_${name}_fortran_exe) elseif(WIN32 AND DEFINED ENV{CI}) + # For Windows we need to provide the full path to the executable since it is installed to a different directory add_test (NAME example_${name}_fortran_with_gdb COMMAND gdb -batch --command=${CMAKE_BINARY_DIR}/cmdfile.gdb ${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}/example_${name}_fortran_exe.exe) endif() From 77c6b787e859b35eae46631483b933efafbd8b2f Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 19:27:49 +0800 Subject: [PATCH 12/13] -changed "progress" to "best point so far" to make it clear what x, f, etc. represent -fixed bug in how cobyla reported nonlinear constraints -removed 'early' from the termination string -fixed comments in cobyla code related to the size of constraint arrays --- c/examples/bobyqa/bobyqa_example.c | 5 +++-- c/examples/cobyla/cobyla_example.c | 5 ++--- c/examples/lincoa/lincoa_example.c | 2 +- c/examples/newuoa/newuoa_example.c | 3 ++- c/examples/uobyqa/uobyqa_example.c | 5 +++-- c/prima.c | 2 +- fortran/cobyla/cobyla.f90 | 4 ++-- fortran/cobyla/cobylb.f90 | 4 ++-- fortran/common/message.f90 | 2 +- 9 files changed, 17 insertions(+), 15 deletions(-) diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index cd8adfae32..bd028a715f 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -15,10 +15,11 @@ static void fun(const double x[], double *f, const void *data) static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) { (void)n; - printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); - *terminate = 0; + (void)cstrv; (void)m_nlcon; (void)nlconstr; + printf("best point so far: x=[%g;%g] f=%g nf=%d tr=%d\n", x[0], x[1], f, nf, tr); + *terminate = 0; } int main(int argc, char * argv[]) diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c index ee0b0f6baa..25c6c46d46 100644 --- a/c/examples/cobyla/cobyla_example.c +++ b/c/examples/cobyla/cobyla_example.c @@ -20,10 +20,9 @@ static void fun(const double x[], double *f, double constr[], const void *data) static void callback(const int n, const double x[], const double f, const int nf, const int tr, const double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) { (void)n; - printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); - *terminate = 0; (void)m_nlcon; - (void)nlconstr; + printf("best point so far: x=[%g;%g] f=%g cstrv=%g nlconstr=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nlconstr[0], nf, tr); + *terminate = 0; } diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c index a1b4626ada..db19d93637 100644 --- a/c/examples/lincoa/lincoa_example.c +++ b/c/examples/lincoa/lincoa_example.c @@ -15,9 +15,9 @@ static void fun(const double x[], double *f, const void *data) static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate) { (void)n; - printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); (void)m_nlcon; (void)nlconstr; + printf("best point so far: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); *terminate = 0; } diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c index c05d370b45..096287815b 100644 --- a/c/examples/newuoa/newuoa_example.c +++ b/c/examples/newuoa/newuoa_example.c @@ -15,9 +15,10 @@ static void fun(const double x[], double *f, const void *data) static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate) { (void)n; - printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); + (void)cstrv; (void)m_nlcon; (void)nlconstr; + printf("best point so far: x=[%g;%g] f=%g nf=%d tr=%d\n", x[0], x[1], f, nf, tr); *terminate = 0; } diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c index 90f4e02c22..ddd1d36510 100644 --- a/c/examples/uobyqa/uobyqa_example.c +++ b/c/examples/uobyqa/uobyqa_example.c @@ -15,10 +15,11 @@ static void fun(const double x[], double *f, const void *data) static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate) { (void)n; - printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr); - *terminate = 0; + (void)cstrv; (void)m_nlcon; (void)nlconstr; + printf("best point so far: x=[%g;%g] f=%g nf=%d tr=%d\n", x[0], x[1], f, nf, tr); + *terminate = 0; } int main(int argc, char * argv[]) diff --git a/c/prima.c b/c/prima.c index c250c6d7ae..f6d9602f9f 100644 --- a/c/prima.c +++ b/c/prima.c @@ -262,7 +262,7 @@ const char *prima_get_rc_string(const prima_rc_t rc) case PRIMA_ZERO_LINEAR_CONSTRAINT: return "One of the linear constraints has a zero gradient"; case PRIMA_CALLBACK_TERMINATE: - return "Callback function requested early termination of optimization"; + return "Callback function requested termination of optimization"; case PRIMA_INVALID_INPUT: return "Invalid input"; case PRIMA_ASSERTION_FAILS: diff --git a/fortran/cobyla/cobyla.f90 b/fortran/cobyla/cobyla.f90 index bde122dc48..025cd23181 100644 --- a/fortran/cobyla/cobyla.f90 +++ b/fortran/cobyla/cobyla.f90 @@ -353,8 +353,8 @@ subroutine cobyla(calcfc, m_nlcon, x, & real(RP), allocatable :: bineq_loc(:) ! Bineq_LOC(Mineq) real(RP), allocatable :: bvec(:) ! BVEC(M_LCON) real(RP), allocatable :: chist_loc(:) ! CHIST_LOC(MAXCHIST) -real(RP), allocatable :: conhist_loc(:, :) ! CONHIST_LOC(M_NLCON, MAXCONHIST) -real(RP), allocatable :: constr_loc(:) ! CONSTR_LOC(M_NLCON) +real(RP), allocatable :: conhist_loc(:, :) ! CONHIST_LOC(M, MAXCONHIST) +real(RP), allocatable :: constr_loc(:) ! CONSTR_LOC(M) real(RP), allocatable :: fhist_loc(:) ! FHIST_LOC(MAXFHIST) real(RP), allocatable :: xhist_loc(:, :) ! XHIST_LOC(N, MAXXHIST) real(RP), allocatable :: xl_loc(:) ! XL_LOC(N) diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index 2aca8b1365..4953f2fc45 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -218,7 +218,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! Report the current best value, and check if user asks for early termination. terminate = .false. if (present(callback_fcn)) then - call callback_fcn(sim(:, n+1), fval(n+1), nf, 0, cval(n+1), conmat(:, n+1), terminate) + call callback_fcn(sim(:, n+1), fval(n+1), nf, 0, cval(n+1), conmat(m_lcon+1:m, n+1), terminate) if (terminate) then subinfo = CALLBACK_TERMINATE end if @@ -605,7 +605,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! Report the current best value, and check if user asks for early termination. if (present(callback_fcn)) then - call callback_fcn(sim(:, n+1), fval(n+1), nf, tr, cval(n+1), conmat(:, n+1), terminate) + call callback_fcn(sim(:, n+1), fval(n+1), nf, tr, cval(n+1), conmat(m_lcon+1:m, n+1), terminate) if (terminate) then info = CALLBACK_TERMINATE exit diff --git a/fortran/common/message.f90 b/fortran/common/message.f90 index 5f7b99d226..1a1346cc3e 100644 --- a/fortran/common/message.f90 +++ b/fortran/common/message.f90 @@ -125,7 +125,7 @@ subroutine retmsg(solver, info, iprint, nf, f, x, cstrv, constr) case (ZERO_LINEAR_CONSTRAINT) reason = 'one of the linear constraints has a zero gradient' case (CALLBACK_TERMINATE) - reason = 'callback function requested early termination of optimization' + reason = 'callback function requested termination of optimization' case default reason = 'UNKNOWN EXIT FLAG' end select From f191d6d9420d82d5ea0d31f9d02e85dd4189a98e Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Wed, 20 Dec 2023 19:46:52 +0800 Subject: [PATCH 13/13] keep problem definition together, add comments, add bounds to bobyqa example and make the cobyla bounds active --- c/examples/bobyqa/bobyqa_example.c | 10 +++++++++- c/examples/cobyla/cobyla_example.c | 20 ++++++++++++-------- c/examples/lincoa/lincoa_example.c | 16 ++++++++++------ c/examples/newuoa/newuoa_example.c | 4 ++++ c/examples/uobyqa/uobyqa_example.c | 4 ++++ 5 files changed, 39 insertions(+), 15 deletions(-) diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index bd028a715f..a4aa44eb4d 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -28,20 +28,28 @@ int main(int argc, char * argv[]) (void)argv; const int n = 2; double x0[2] = {0.0, 0.0}; + // set up the problem prima_problem_t problem; prima_init_problem(&problem, n); problem.x0 = x0; problem.calfun = &fun; + double xl[] = {0.0, 0.0}; + double xu[] = {1.0, 1.0}; + problem.xl = xl; + problem.xu = xu; + // set up the options prima_options_t options; prima_init_options(&options); options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; options.callback = &callback; + // initialize the result prima_result_t result; + // run the solver const int rc = prima_minimize(PRIMA_BOBYQA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); prima_free_problem(&problem); prima_free_result(&result); - return (fabs(result.x[0]-3)>2e-2 || fabs(result.x[1]-2)>2e-2); + return (fabs(result.x[0]-1)>2e-2 || fabs(result.x[1]-1)>2e-2); } diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c index 25c6c46d46..2b504b195d 100644 --- a/c/examples/cobyla/cobyla_example.c +++ b/c/examples/cobyla/cobyla_example.c @@ -12,7 +12,7 @@ static void fun(const double x[], double *f, double constr[], const void *data) const double x1 = x[0]; const double x2 = x[1]; *f = 5*(x1-3)*(x1-3)+7*(x2-2)*(x2-2)+0.1*(x1+x2)-10; - constr[0] = x1*x1 + x2*x2 - 13;// ||x||^2<=13 + constr[0] = x1*x1 + x2*x2 - 12;// ||x||^2<=12 (void)data; } @@ -32,16 +32,11 @@ int main(int argc, char * argv[]) (void)argv; const int n = 2; double x0[2] = {0.0, 0.0}; + // set up the problem prima_problem_t problem; prima_init_problem(&problem, n); problem.calcfc = &fun; problem.x0 = x0; - prima_options_t options; - prima_init_options(&options); - options.iprint = PRIMA_MSG_EXIT; - options.rhoend= 1e-3; - options.maxfun = 200*n; - options.callback = &callback; problem.m_nlcon = M_NLCON; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; @@ -57,10 +52,19 @@ int main(int argc, char * argv[]) double xu[2] = {6.0, 6.0}; problem.xl = xl; problem.xu = xu; + // set up the options + prima_options_t options; + prima_init_options(&options); + options.iprint = PRIMA_MSG_EXIT; + options.rhoend= 1e-3; + options.maxfun = 200*n; + options.callback = &callback; + // initialize the result prima_result_t result; + // run the solver const int rc = prima_minimize(PRIMA_COBYLA, &problem, &options, &result); printf("x*={%g, %g} f*=%g cstrv=%g nlconstr=%g rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], result.f, result.cstrv, result.nlconstr[0], rc, result.message, result.nf); prima_free_problem(&problem); prima_free_result(&result); - return (fabs(result.x[0]-3)>2e-2 || fabs(result.x[1]-2)>2e-2); + return (fabs(result.x[0]-2.86)>2e-2 || fabs(result.x[1]-1.94)>2e-2); } diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c index db19d93637..588c4c549e 100644 --- a/c/examples/lincoa/lincoa_example.c +++ b/c/examples/lincoa/lincoa_example.c @@ -27,16 +27,11 @@ int main(int argc, char * argv[]) (void)argv; const int n = 2; double x0[2] = {0.0, 0.0}; + // set up the problem prima_problem_t problem; prima_init_problem(&problem, n); problem.calfun = &fun; problem.x0 = x0; - prima_options_t options; - prima_init_options(&options); - options.iprint = PRIMA_MSG_EXIT; - options.rhoend= 1e-3; - options.maxfun = 200*n; - options.callback = &callback; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; double Aineq[3*2] = {1.0, 0.0, @@ -51,7 +46,16 @@ int main(int argc, char * argv[]) double xu[2] = {6.0, 6.0}; problem.xl = xl; problem.xu = xu; + // set up the options + prima_options_t options; + prima_init_options(&options); + options.iprint = PRIMA_MSG_EXIT; + options.rhoend= 1e-3; + options.maxfun = 200*n; + options.callback = &callback; + // initialize the result prima_result_t result; + // run the solver const int rc = prima_minimize(PRIMA_LINCOA, &problem, &options, &result); printf("x*={%g, %g} f*=%g cstrv=%g rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], result.f, result.cstrv, rc, result.message, result.nf); prima_free_problem(&problem); diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c index 096287815b..428c32b74d 100644 --- a/c/examples/newuoa/newuoa_example.c +++ b/c/examples/newuoa/newuoa_example.c @@ -28,17 +28,21 @@ int main(int argc, char * argv[]) (void)argv; const int n = 2; double x0[2] = {0.0, 0.0}; + // set up the problem prima_problem_t problem; prima_init_problem(&problem, n); problem.calfun = &fun; problem.x0 = x0; + // set up the options prima_options_t options; prima_init_options(&options); options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; options.callback = &callback; + // initialize the result prima_result_t result; + // run the solver const int rc = prima_minimize(PRIMA_NEWUOA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); prima_free_problem(&problem); diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c index ddd1d36510..166c424111 100644 --- a/c/examples/uobyqa/uobyqa_example.c +++ b/c/examples/uobyqa/uobyqa_example.c @@ -28,17 +28,21 @@ int main(int argc, char * argv[]) (void)argv; const int n = 2; double x0[2] = {0.0, 0.0}; + // set up the problem prima_problem_t problem; prima_init_problem(&problem, n); problem.calfun = &fun; problem.x0 = x0; + // set up the options prima_options_t options; prima_init_options(&options); options.iprint = PRIMA_MSG_EXIT; options.rhoend= 1e-3; options.maxfun = 200*n; options.callback = &callback; + // initialize the result prima_result_t result; + // run the solver const int rc = prima_minimize(PRIMA_UOBYQA, &problem, &options, &result); printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf); prima_free_problem(&problem);