diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt index c6b95fe345..57e1859dc8 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,6 @@ orthtol nouninit libgfortran chocolatey +fcn +BINDIR +cmdfile diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index adfe78991a..d672a1814f 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,9 @@ 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 env: FC: ${{ steps.setup-fortran.outputs.fc }} shell: bash @@ -126,7 +129,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 +160,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 +178,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..005aac3a70 100644 --- a/c/CMakeLists.txt +++ b/c/CMakeLists.txt @@ -6,6 +6,7 @@ if (WIN32) set_target_properties(primac PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin) endif() + target_include_directories (primac PUBLIC $ $ @@ -35,7 +36,21 @@ macro (prima_add_c_test name) if (WIN32) set_target_properties(example_${name}_c_exe PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${PROJECT_BINARY_DIR}/bin) endif() + + # 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 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() + add_dependencies(examples example_${name}_c_exe) endmacro () diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90 index 7cb90a4466..069f351c14 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,12 @@ 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) +! 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() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -62,10 +70,20 @@ 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 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 + 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)) @@ -79,16 +97,99 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta ! 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 use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj 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 -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 + +!--------------------------------------------------------------------------------------------------! +! 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, 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 +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), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: 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_sub_loc +integer(C_INT) :: m_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 +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)) + +! 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 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 +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 + 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..961400bf1d 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,12 @@ 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 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() ! 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 +95,24 @@ 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 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, & + & 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)) @@ -110,17 +128,103 @@ 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 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 use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobjcon 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(:) -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 + +!--------------------------------------------------------------------------------------------------! +! 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, 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 +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), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: 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_sub_loc +integer(C_INT) :: m_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 +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)) + +! 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 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 +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 + + end subroutine cobyla_c diff --git a/c/examples/bobyqa/bobyqa_example.c b/c/examples/bobyqa/bobyqa_example.c index d871944961..a4aa44eb4d 100644 --- a/c/examples/bobyqa/bobyqa_example.c +++ b/c/examples/bobyqa/bobyqa_example.c @@ -12,25 +12,44 @@ 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; + (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[]) { (void)argc; (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 8053eabc2b..2b504b195d 100644 --- a/c/examples/cobyla/cobyla_example.c +++ b/c/examples/cobyla/cobyla_example.c @@ -12,25 +12,31 @@ 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; } + +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; + (void)m_nlcon; + 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; +} + + int main(int argc, char * argv[]) { (void)argc; (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; problem.m_nlcon = M_NLCON; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; @@ -46,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 e5aa89d05e..588c4c549e 100644 --- a/c/examples/lincoa/lincoa_example.c +++ b/c/examples/lincoa/lincoa_example.c @@ -12,21 +12,26 @@ 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; + (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; +} + int main(int argc, char * argv[]) { (void)argc; (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; // x1<=4, x2<=3, x1+x2<=10 problem.m_ineq = 3; double Aineq[3*2] = {1.0, 0.0, @@ -41,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 00c4aa709b..428c32b74d 100644 --- a/c/examples/newuoa/newuoa_example.c +++ b/c/examples/newuoa/newuoa_example.c @@ -12,22 +12,37 @@ 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; + (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[]) { (void)argc; (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 be53e87c21..166c424111 100644 --- a/c/examples/uobyqa/uobyqa_example.c +++ b/c/examples/uobyqa/uobyqa_example.c @@ -12,22 +12,37 @@ 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; + (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[]) { (void)argc; (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); diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h index bb38c17961..e2d63df0ea 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_CALLBACK_TERMINATE = 30, PRIMA_INVALID_INPUT = 100, PRIMA_ASSERTION_FAILS = 101, PRIMA_VALIDATION_FAILS = 102, @@ -67,12 +70,12 @@ 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 * 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 @@ -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 { @@ -102,9 +120,12 @@ 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) + prima_callback_t callback; + } prima_options_t; /* Initialize problem */ diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90 index b41746394d..5e16b673c8 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,12 @@ 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 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() ! 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 +88,25 @@ 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 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, & + & 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)) @@ -102,16 +121,99 @@ 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 use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj 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 -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 + +!--------------------------------------------------------------------------------------------------! +! 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, 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 +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), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: 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_sub_loc +integer(C_INT) :: m_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 +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)) + +! 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 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 +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 + end subroutine lincoa_c diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90 index 5177d1970c..7585b55c30 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,12 @@ 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 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() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -55,10 +62,20 @@ 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 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 + 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)) @@ -72,16 +89,99 @@ 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 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj 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 -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 + +!--------------------------------------------------------------------------------------------------! +! 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, 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 +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), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: 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_sub_loc +integer(C_INT) :: m_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 +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)) + +! 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 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 +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 + end subroutine newuoa_c diff --git a/c/prima.c b/c/prima.c index 38d2d96127..f6d9602f9f 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_CALLBACK_TERMINATE: + return "Callback function requested termination of optimization"; case PRIMA_INVALID_INPUT: return "Invalid input"; case PRIMA_ASSERTION_FAILS: diff --git a/c/tests/CMakeLists.txt b/c/tests/CMakeLists.txt index 6fc0d612df..0c9ff3d5b1 100644 --- a/c/tests/CMakeLists.txt +++ b/c/tests/CMakeLists.txt @@ -9,11 +9,32 @@ 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) + # 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 + 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. + 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}) + # 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) + 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/c/uobyqa_c.f90 b/c/uobyqa_c.f90 index 07d920dc13..3983b7dc2c 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 @@ -32,6 +32,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 +45,12 @@ 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 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() ! Read the inputs and convert them to the Fortran side types x_loc = real(x, kind(x_loc)) @@ -52,10 +59,20 @@ 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 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 + 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)) @@ -69,16 +86,99 @@ 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 out of 5 algorithms; COBYLA requires a slightly different +! signature. !--------------------------------------------------------------------------------------------------! subroutine calfun(x_sub, f_sub) +use, intrinsic :: iso_c_binding, only : C_DOUBLE use, non_intrinsic :: consts_mod, only : RP -use, non_intrinsic :: cintrf_mod, only : evalcobj 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 -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 + +!--------------------------------------------------------------------------------------------------! +! 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, 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 +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), optional :: cstrv_sub +real(RP), intent(in), optional :: nlconstr_sub(:) +logical, intent(out), optional :: 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_sub_loc +integer(C_INT) :: m_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 +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)) + +! 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 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 +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 + end subroutine uobyqa_c diff --git a/fortran/CMakeLists.txt b/fortran/CMakeLists.txt index 1a926fd4b4..80469a54e2 100644 --- a/fortran/CMakeLists.txt +++ b/fortran/CMakeLists.txt @@ -114,7 +114,21 @@ macro (prima_add_f_test name) if (PRIMA_ENABLE_EXAMPLES) set_target_properties (example_${name}_fortran_exe PROPERTIES EXCLUDE_FROM_ALL FALSE) endif () + + # 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 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() + add_dependencies(examples example_${name}_fortran_exe) endmacro () diff --git a/fortran/bobyqa/bobyqa.f90 b/fortran/bobyqa/bobyqa.f90 index 24e7f71bc5..8092eae92d 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 optionally 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..2a3a8aabac 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, 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 +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,15 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) + if (terminate) then + subinfo = CALLBACK_TERMINATE + 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 +582,16 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) + if (terminate) then + info = CALLBACK_TERMINATE + 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..025cd23181 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 optionally 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 @@ -349,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) @@ -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..4953f2fc45 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, 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 +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(m_lcon+1:m, n+1), terminate) + if (terminate) then + subinfo = CALLBACK_TERMINATE + 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(m_lcon+1:m, n+1), terminate) + if (terminate) then + info = CALLBACK_TERMINATE + 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..e474a5b98b 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 :: CALLBACK_TERMINATE 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 :: 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 779259e6a6..1a1346cc3e 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, CALLBACK_TERMINATE 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 (CALLBACK_TERMINATE) + reason = 'callback function requested termination of optimization' case default reason = 'UNKNOWN EXIT FLAG' end select diff --git a/fortran/common/pintrf.f90 b/fortran/common/pintrf.f90 index 8a610e0c60..cbc1f1490c 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), optional :: cstrv + real(RP), intent(in), optional :: nlconstr(:) + logical, intent(out), optional :: terminate + end subroutine CALLBACK +end interface end module pintrf_mod diff --git a/fortran/lincoa/lincoa.f90 b/fortran/lincoa/lincoa.f90 index c6aac8a3c7..8cd7a54f8b 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 optionally 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..35ce812e44 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, 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 -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,15 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), terminate=terminate) + if (terminate) then + subinfo = CALLBACK_TERMINATE + 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 +646,16 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), terminate=terminate) + if (terminate) then + info = CALLBACK_TERMINATE + 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..05318d4194 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 optionally 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..8b948f287f 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, 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 +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,15 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) + if (terminate) then + subinfo = CALLBACK_TERMINATE + end if +end if + ! Initialize X and F according to KOPT. x = xbase + xpt(:, kopt) f = fval(kopt) @@ -574,6 +585,16 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) + if (terminate) then + info = CALLBACK_TERMINATE + 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/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/uobyqa.f90 b/fortran/uobyqa/uobyqa.f90 index 55261e3398..546c8d84d6 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 optionally 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..244782ab5e 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, 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 -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,16 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, terminate=terminate) + if (terminate) then + subinfo = CALLBACK_TERMINATE + end if +end if + ! Initialize X and F according to KOPT. x = xbase + xpt(:, kopt) f = fval(kopt) @@ -458,6 +470,15 @@ 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 + call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, terminate=terminate) + if (terminate) then + info = CALLBACK_TERMINATE + 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