From c1297dd6da2d74ea2a64d91542737df1ef40c4f6 Mon Sep 17 00:00:00 2001
From: Nickolai Belakovski <nbelakovski@gmail.com>
Date: Thu, 23 Nov 2023 18:10:48 +0800
Subject: [PATCH] Manually applied jschueller's closed callback PR

---
 .github/actions/spelling/allow.txt |  4 ++
 c/bobyqa_c.f90                     | 26 ++++++++++++-
 c/cintrf.f90                       | 62 +++++++++++++++++++++++++++++-
 c/cobyla_c.f90                     | 27 +++++++++++--
 c/examples/bobyqa/bobyqa_example.c | 10 +++++
 c/examples/cobyla/cobyla_example.c | 12 ++++++
 c/examples/lincoa/lincoa_example.c | 10 +++++
 c/examples/newuoa/newuoa_example.c | 10 +++++
 c/examples/uobyqa/uobyqa_example.c | 10 +++++
 c/include/prima/prima.h            | 23 ++++++++++-
 c/lincoa_c.f90                     | 26 ++++++++++++-
 c/newuoa_c.f90                     | 24 +++++++++++-
 c/prima.c                          | 22 ++++++-----
 c/uobyqa_c.f90                     | 26 ++++++++++++-
 fortran/bobyqa/bobyqa.f90          | 14 +++++--
 fortran/bobyqa/bobyqb.f90          | 22 +++++++++--
 fortran/bobyqa/initialize.f90      | 18 +++++++--
 fortran/cobyla/cobyla.f90          | 18 +++++++--
 fortran/cobyla/cobylb.f90          | 21 +++++++---
 fortran/cobyla/initialize.f90      | 15 ++++++--
 fortran/common/infos.f90           |  3 +-
 fortran/common/pintrf.f90          | 32 ++++++++++++++-
 fortran/lincoa/initialize.f90      | 17 ++++++--
 fortran/lincoa/lincoa.f90          | 18 +++++++--
 fortran/lincoa/lincob.f90          | 21 ++++++++--
 fortran/newuoa/initialize.f90      | 18 +++++++--
 fortran/newuoa/newuoa.f90          | 17 ++++++--
 fortran/newuoa/newuob.f90          | 21 ++++++++--
 fortran/uobyqa/initialize.f90      | 17 ++++++--
 fortran/uobyqa/uobyqa.f90          | 18 +++++++--
 fortran/uobyqa/uobyqb.f90          | 20 ++++++++--
 31 files changed, 526 insertions(+), 76 deletions(-)

diff --git a/.github/actions/spelling/allow.txt b/.github/actions/spelling/allow.txt
index f7ec8ee176..8f0a61e5bc 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,6 +2049,7 @@ COBJCON
 cobjfun
 cobjfuncon
 constrc
+evalcallback
 evalcobj
 evalcobjcon
 execstack
@@ -2133,3 +2135,5 @@ TWOBARS
 lang
 archnorma
 orthtol
+fcb
+fcn
diff --git a/c/bobyqa_c.f90 b/c/bobyqa_c.f90
index 7cb90a4466..eeb83394a6 100644
--- a/c/bobyqa_c.f90
+++ b/c/bobyqa_c.f90
@@ -13,7 +13,8 @@ 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)
+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
 use, non_intrinsic :: cintrf_mod, only : COBJ
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -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
@@ -65,7 +67,7 @@ subroutine bobyqa_c(cobj_ptr, data_ptr, n, x, f, xl, xu, nf, rhobeg, rhoend, fta
 
 ! 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)
+    & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -89,6 +91,26 @@ subroutine calfun(x_sub, f_sub)
 call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub)
 end subroutine calfun
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
+
 end subroutine bobyqa_c
 
 
diff --git a/c/cintrf.f90 b/c/cintrf.f90
index f837dcc9f7..08c4ec86c9 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, evalcobj, evalcobjcon, evalcallback
 
 
 abstract interface
@@ -31,6 +31,22 @@ subroutine COBJCON(x, f, constr, data_ptr) bind(c)
     type(C_PTR), intent(in), value :: data_ptr
     end subroutine COBJCON
 
+    subroutine CCALLBACK(n, x, f, nf, tr, cstrv, m_nlcon, nlconstr, terminate) bind(c)
+    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
+
+
 end interface
 
 
@@ -93,4 +109,48 @@ subroutine evalcobjcon(cobjcon_ptr, data_ptr, x, f, constr)
 
 end subroutine evalcobjcon
 
+
+subroutine evalcallback(fcb_ptr, x, f, nf, tr, cstrv, nlconstr, terminate)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_FUNPTR, C_F_PROCPOINTER, C_BOOL, C_INT
+implicit none
+type(C_FUNPTR), intent(in) :: fcb_ptr
+real(RP), intent(in) :: x(:)
+real(RP), intent(in) :: f
+integer(IK), intent(in) :: nf
+integer(IK), intent(in) :: tr
+real(RP), intent(in) :: cstrv
+real(RP), intent(in) :: nlconstr(:)
+logical, intent(out) :: terminate
+
+! Local variables
+procedure(CCALLBACK), pointer :: cb_ptr
+integer(C_INT) :: n_loc
+real(C_DOUBLE) :: x_loc(size(x))
+real(C_DOUBLE) :: f_loc
+integer(C_INT) :: nf_loc
+integer(C_INT) :: tr_loc
+real(C_DOUBLE) :: cstrv_loc
+real(C_DOUBLE) :: nlconstr_loc(size(nlconstr))
+logical(C_BOOL) :: terminate_loc
+
+! Read the inputs and convert them to the types specified in CCALLBACK
+n_loc = size(x)
+x_loc = real(x, kind(x_loc))
+f_loc = real(f, kind(f_loc))
+nf_loc = int(nf, kind(nf_loc))
+tr_loc = int(tr, kind(tr_loc))
+cstrv_loc = real(cstrv, kind(cstrv_loc))
+nlconstr_loc = real(nlconstr, kind(nlconstr_loc))
+call C_F_PROCPOINTER(fcb_ptr, cb_ptr)
+
+! Call the C objective function
+call cb_ptr(n_loc, x_loc, f_loc, nf_loc, tr_loc, cstrv_loc, size(nlconstr_loc), nlconstr_loc, terminate_loc)
+
+! Write the output
+terminate = logical(terminate_loc, kind(terminate))
+
+end subroutine evalcallback
+
+
 end module cintrf_mod
diff --git a/c/cobyla_c.f90 b/c/cobyla_c.f90
index 35714a1884..083b26e470 100644
--- a/c/cobyla_c.f90
+++ b/c/cobyla_c.f90
@@ -12,8 +12,8 @@ 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)
+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
 use, non_intrinsic :: cintrf_mod, only : COBJCON
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -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
@@ -94,7 +95,7 @@ subroutine cobyla_c(m_nlcon, cobjcon_ptr, data_ptr, n, x, f, cstrv, nlconstr, m_
     & 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)
+    & iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -121,6 +122,26 @@ subroutine calcfc(x_sub, f_sub, constr_sub)
 call evalcobjcon(cobjcon_ptr, data_ptr, x_sub, f_sub, constr_sub)
 end subroutine calcfc
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+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..cd8adfae32 100644
--- a/c/examples/bobyqa/bobyqa_example.c
+++ b/c/examples/bobyqa/bobyqa_example.c
@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_BOBYQA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/examples/cobyla/cobyla_example.c b/c/examples/cobyla/cobyla_example.c
index 8053eabc2b..ee0b0f6baa 100644
--- a/c/examples/cobyla/cobyla_example.c
+++ b/c/examples/cobyla/cobyla_example.c
@@ -16,6 +16,17 @@ static void fun(const double x[], double *f, double constr[], const void *data)
   (void)data;
 }
 
+
+static void callback(const int n, const double x[], const double f, const int nf, const int tr, const double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -31,6 +42,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   problem.m_nlcon = M_NLCON;
   // x1<=4, x2<=3, x1+x2<=10
   problem.m_ineq = 3;
diff --git a/c/examples/lincoa/lincoa_example.c b/c/examples/lincoa/lincoa_example.c
index e5aa89d05e..a1b4626ada 100644
--- a/c/examples/lincoa/lincoa_example.c
+++ b/c/examples/lincoa/lincoa_example.c
@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  (void)m_nlcon;
+  (void)nlconstr;
+  *terminate = 0;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   // x1<=4, x2<=3, x1+x2<=10
   problem.m_ineq = 3;
   double Aineq[3*2] = {1.0, 0.0,
diff --git a/c/examples/newuoa/newuoa_example.c b/c/examples/newuoa/newuoa_example.c
index 00c4aa709b..c05d370b45 100644
--- a/c/examples/newuoa/newuoa_example.c
+++ b/c/examples/newuoa/newuoa_example.c
@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  (void)m_nlcon;
+  (void)nlconstr;
+  *terminate = 0;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_NEWUOA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/examples/uobyqa/uobyqa_example.c b/c/examples/uobyqa/uobyqa_example.c
index be53e87c21..90f4e02c22 100644
--- a/c/examples/uobyqa/uobyqa_example.c
+++ b/c/examples/uobyqa/uobyqa_example.c
@@ -12,6 +12,15 @@ static void fun(const double x[], double *f, const void *data)
   (void)data;
 }
 
+static void callback(int n, const double x[], double f, int nf, int tr, double cstrv, const int m_nlcon, const double nlconstr[], bool *terminate)
+{
+  (void)n;
+  printf("progress: x=[%g;%g] f=%g cstrv=%g nf=%d tr=%d\n", x[0], x[1], f, cstrv, nf, tr);
+  *terminate = 0;
+  (void)m_nlcon;
+  (void)nlconstr;
+}
+
 int main(int argc, char * argv[])
 {
   (void)argc;
@@ -27,6 +36,7 @@ int main(int argc, char * argv[])
   options.iprint = PRIMA_MSG_EXIT;
   options.rhoend= 1e-3;
   options.maxfun = 200*n;
+  options.callback = &callback;
   prima_result_t result;
   const int rc = prima_minimize(PRIMA_UOBYQA, &problem, &options, &result);
   printf("x*={%g, %g} rc=%d msg='%s' evals=%d\n", result.x[0], result.x[1], rc, result.message, result.nf);
diff --git a/c/include/prima/prima.h b/c/include/prima/prima.h
index bb38c17961..babc335b2f 100644
--- a/c/include/prima/prima.h
+++ b/c/include/prima/prima.h
@@ -3,6 +3,8 @@
 #ifndef PRIMA_H
 #define PRIMA_H
 
+#include <stdbool.h>
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -49,6 +51,7 @@ typedef enum
   PRIMA_NO_SPACE_BETWEEN_BOUNDS = 6,
   PRIMA_DAMAGING_ROUNDING = 7,
   PRIMA_ZERO_LINEAR_CONSTRAINT = 8,
+  PRIMA_USER_STOP = 9,
   PRIMA_INVALID_INPUT = 100,
   PRIMA_ASSERTION_FAILS = 101,
   PRIMA_VALIDATION_FAILS = 102,
@@ -67,7 +70,7 @@ PRIMAC_API
 const char *prima_get_rc_string(const prima_rc_t rc);
 
 /*
- * A function as required by solvers
+ * The objective function as required by solvers
  *
  * x     : on input, then vector of variables (should not be modified)
  * f     : on output, the value of the function
@@ -80,6 +83,21 @@ const char *prima_get_rc_string(const prima_rc_t rc);
 typedef void (*prima_obj_t)(const double x[], double *f, const void *data);
 typedef void (*prima_objcon_t)(const double x[], double *f, double constr[], const void *data);
 
+/* An optional callback function to report algorithm progress
+ *
+ * n     : number of variables
+ * x     : the current best point
+ * f     : the function value of the best point
+ * nf    : number of objective function calls
+ * tr    : iteration number
+ * cstrv : the constraint value verified by the current best point
+ * m_nlcon : number of non-linear constraints (cobyla only)
+ * nlconstr : non-linear constraints values verified by the current best point (cobyla only)
+ * terminate : a boolean to ask from early optimization exit
+*/
+typedef void (*prima_callback_t)(const int n, const double x[], const double f, int nf, int tr,
+                               const double cstrv, int m_nlcon, const double nlconstr[], bool *terminate);
+
 
 typedef struct {
   
@@ -105,6 +123,9 @@ typedef struct {
   // user-data, will be passed through the objective function callback
   void *data;
 
+  // callback function to report algorithm progress (default=NULL)
+  prima_callback_t callback;
+
 } prima_options_t;
 
 /* Initialize problem */
diff --git a/c/lincoa_c.f90 b/c/lincoa_c.f90
index b41746394d..0140decdf3 100644
--- a/c/lincoa_c.f90
+++ b/c/lincoa_c.f90
@@ -14,7 +14,7 @@ module lincoa_c_mod
 
 
 subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_eq, Aeq, beq, xl, xu, &
-    & nf, rhobeg, rhoend, ftarget, maxfun, npt, iprint, info) bind(C)
+    & 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
 use, non_intrinsic :: cintrf_mod, only : COBJ
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -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
@@ -87,7 +88,8 @@ subroutine lincoa_c(cobj_ptr, data_ptr, n, x, f, cstrv, m_ineq, Aineq, bineq, m_
     & 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)
+    & ftarget=ftarget_loc, maxfun=maxfun_loc, npt=npt_loc, &
+    & iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -112,6 +114,26 @@ subroutine calfun(x_sub, f_sub)
 call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub)
 end subroutine calfun
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
+
 end subroutine lincoa_c
 
 
diff --git a/c/newuoa_c.f90 b/c/newuoa_c.f90
index 5177d1970c..a579584d95 100644
--- a/c/newuoa_c.f90
+++ b/c/newuoa_c.f90
@@ -12,7 +12,7 @@ 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)
+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
 use, non_intrinsic :: cintrf_mod, only : COBJ
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -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
@@ -58,7 +59,7 @@ subroutine newuoa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 
 ! 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)
+    & maxfun=maxfun_loc, npt=npt_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -82,6 +83,25 @@ subroutine calfun(x_sub, f_sub)
 call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub)
 end subroutine calfun
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
 end subroutine newuoa_c
 
 
diff --git a/c/prima.c b/c/prima.c
index 38d2d96127..a7ecfa75b2 100644
--- a/c/prima.c
+++ b/c/prima.c
@@ -44,19 +44,19 @@ int cobyla_c(const int m_nlcon, const prima_objcon_t calcfc, const void *data, c
              const int m_eq, const double Aeq[], const double beq[],
              const double xl[], const double xu[],
              double f0, const double nlconstr0[],
-             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, int *info);
+             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, const prima_callback_t callback, int *info);
 int bobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f, const double xl[], const double xu[],
-             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, int *info);
+             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *info);
 int newuoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f,
-             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, int *info);
+             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *info);
 int uobyqa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f,
-             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, int *info);
+             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int iprint, const prima_callback_t callback, int *info);
 int lincoa_c(prima_obj_t calfun, const void *data, const int n, double x[], double *f,
              double *cstrv,
              const int m_ineq, const double Aineq[], const double bineq[],
              const int m_eq, const double Aeq[], const double beq[],
              const double xl[], const double xu[],
-             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, int *info);
+             int *nf, const double rhobeg, const double rhoend, const double ftarget, const int maxfun, const int npt, const int iprint, const prima_callback_t callback, int *info);
 
 int prima_check_problem(prima_problem_t *problem, prima_options_t *options, int alloc_bounds, int use_constr)
 {
@@ -203,27 +203,27 @@ int prima_minimize(const prima_algorithm_t algorithm, prima_problem_t *problem,
     switch (algorithm)
     {
       case PRIMA_BOBYQA:
-        bobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, &info);
+        bobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info);
       break;
 
       case PRIMA_COBYLA:
         cobyla_c(problem->m_nlcon, problem->calcfc, options->data, problem->n, result->x, &(result->f), &(result->cstrv), result->nlconstr,
               problem->m_ineq, problem->Aineq, problem->bineq, problem->m_eq, problem->Aeq, problem->beq,
-              problem->xl, problem->xu, problem->f0, problem->nlconstr0, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, &info);
+              problem->xl, problem->xu, problem->f0, problem->nlconstr0, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, options->callback, &info);
       break;
 
       case PRIMA_LINCOA:
         lincoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->cstrv),
             problem->m_ineq, problem->Aineq, problem->bineq, problem->m_eq, problem->Aeq, problem->beq,
-            problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, &info);
+            problem->xl, problem->xu, &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info);
       break;
 
       case PRIMA_NEWUOA:
-        newuoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, &info);
+        newuoa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->npt, options->iprint, options->callback, &info);
       break;
 
       case PRIMA_UOBYQA:
-        uobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, &info);
+        uobyqa_c(problem->calfun, options->data, problem->n, result->x, &(result->f), &(result->nf), options->rhobeg, options->rhoend, options->ftarget, options->maxfun, options->iprint, options->callback, &info);
       break;
 
       default:
@@ -261,6 +261,8 @@ const char *prima_get_rc_string(const prima_rc_t rc)
       return "Rounding errors are becoming damaging";
     case PRIMA_ZERO_LINEAR_CONSTRAINT:
       return "One of the linear constraints has a zero gradient";
+    case PRIMA_USER_STOP:
+      return "User requested end of optimization";
     case PRIMA_INVALID_INPUT:
       return "Invalid input";
     case PRIMA_ASSERTION_FAILS:
diff --git a/c/uobyqa_c.f90 b/c/uobyqa_c.f90
index 07d920dc13..ff1335b669 100644
--- a/c/uobyqa_c.f90
+++ b/c/uobyqa_c.f90
@@ -12,7 +12,7 @@ module uobyqa_c_mod
 contains
 
 
-subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, maxfun, iprint, info) bind(C)
+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
 use, non_intrinsic :: cintrf_mod, only : COBJ
 use, non_intrinsic :: consts_mod, only : RP, IK
@@ -22,6 +22,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 ! Compulsory arguments
 type(C_FUNPTR), intent(IN), value :: cobj_ptr
 type(C_PTR), intent(in), value :: data_ptr
+
 integer(C_INT), intent(in), value :: n
 ! We cannot use assumed-shape arrays for C interoperability
 real(C_DOUBLE), intent(inout) :: x(n)
@@ -32,6 +33,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 real(C_DOUBLE), intent(in), value :: ftarget
 integer(C_INT), intent(in), value :: maxfun
 integer(C_INT), intent(in), value :: iprint
+type(C_FUNPTR), intent(in), value :: callback_ptr
 integer(C_INT), intent(out) :: info
 
 ! Local variables
@@ -55,7 +57,7 @@ subroutine uobyqa_c(cobj_ptr, data_ptr, n, x, f, nf, rhobeg, rhoend, ftarget, ma
 
 ! 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)
+    & maxfun=maxfun_loc, iprint=iprint_loc, callback_fcn=callback_fcn, info=info_loc)
 
 ! Write the outputs
 x = real(x_loc, kind(x))
@@ -79,6 +81,26 @@ subroutine calfun(x_sub, f_sub)
 call evalcobj(cobj_ptr, data_ptr, x_sub, f_sub)
 end subroutine calfun
 
+
+subroutine callback_fcn(x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+use, non_intrinsic :: consts_mod, only : RP, IK
+use, non_intrinsic :: cintrf_mod, only : evalcallback
+use, intrinsic :: iso_c_binding, only : C_DOUBLE, C_INT, C_ASSOCIATED
+implicit none
+real(RP), intent(in) :: x_sub(:)
+real(RP), intent(in) :: f_sub
+integer(IK), intent(in) :: nf_sub
+integer(IK), intent(in) :: tr_sub
+real(RP), intent(in) :: cstrv_sub
+real(RP), intent(in) :: nlconstr_sub(:)
+logical, intent(out) :: terminate_sub
+terminate_sub = .false.
+if (C_ASSOCIATED(callback_ptr)) then
+    call evalcallback(callback_ptr, x_sub, f_sub, nf_sub, tr_sub, cstrv_sub, nlconstr_sub, terminate_sub)
+end if
+end subroutine callback_fcn
+
+
 end subroutine uobyqa_c
 
 
diff --git a/fortran/bobyqa/bobyqa.f90 b/fortran/bobyqa/bobyqa.f90
index 24e7f71bc5..e64d07645b 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
@@ -198,7 +198,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, NO_OP_CALLBACK
 use, non_intrinsic :: preproc_mod, only : preproc
 use, non_intrinsic :: string_mod, only : num2str
 
@@ -212,6 +212,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
@@ -237,6 +238,7 @@ subroutine bobyqa(calfun, x, &
 ! Local variables
 character(len=*), parameter :: solver = 'BOBYQA'
 character(len=*), parameter :: srname = 'BOBYQA'
+procedure(CALLBACK), pointer :: callback_fcn_loc
 integer(IK) :: info_loc
 integer(IK) :: iprint_loc
 integer(IK) :: k
@@ -402,6 +404,12 @@ subroutine bobyqa(calfun, x, &
     honour_x0_loc = (.not. has_rhobeg)
 end if
 
+if (present(callback_fcn)) then
+    callback_fcn_loc => callback_fcn
+else
+    callback_fcn_loc => NO_OP_CALLBACK
+end if
+
 ! Preprocess the inputs in case some of them are invalid. It does nothing if all inputs are valid.
 call preproc(solver, n, iprint_loc, maxfun_loc, maxhist_loc, ftarget_loc, rhobeg_loc, rhoend_loc, &
     & npt=npt_loc, eta1=eta1_loc, eta2=eta2_loc, gamma1=gamma1_loc, gamma2=gamma2_loc, &
@@ -417,7 +425,7 @@ 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)
+    & fhist_loc, xhist_loc, callback_fcn_loc, info_loc)
 !--------------------------------------------------------------------------------------------------!
 
 ! Write the outputs.
diff --git a/fortran/bobyqa/bobyqb.f90 b/fortran/bobyqa/bobyqb.f90
index 57fb718116..a3fd9fdfe8 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, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine performs the major calculations of BOBYQA.
 !
@@ -82,10 +82,10 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist, rangehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_finite, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED
+use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP
 use, non_intrinsic :: linalg_mod, only : norm
 use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : quadinc, calden, calvlag!, errquad
 use, non_intrinsic :: ratio_mod, only : redrat
 use, non_intrinsic :: redrho_mod, only : redrho
@@ -103,6 +103,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 integer(IK), intent(in) :: npt
@@ -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))
@@ -214,7 +216,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 
 ! Initialize XBASE, XPT, SL, SU, FVAL, and KOPT, together with the history, NF, and IJ.
 call initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x, ij, kopt, nf, fhist, fval, &
-    & sl, su, xbase, xhist, xpt, subinfo)
+    & sl, su, xbase, xhist, xpt, callback_fcn, subinfo)
 
 ! Initialize X and F according to KOPT.
 x = xinbd(xbase, xpt(:, kopt), xl, xu, sl, su)  ! In precise arithmetic, X = XBASE + XOPT.
@@ -272,6 +274,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 ! model but the geometry step is not invoked. Thus the following MAXTR is unlikely to be reached.
 maxtr = max(maxfun, 2_IK * maxfun)  ! MAX: precaution against overflow, which will make 2*MAXFUN < 0.
 info = MAXTR_REACHED
+terminate = .false.
 
 ! Begin the iterative procedure.
 ! After solving a trust-region subproblem, we use three boolean variables to control the workflow.
@@ -571,6 +574,17 @@ 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.
+    ! Since this algorithm is a feasible method, the constraint violation is always 0.
+    ! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set
+    ! to an empty array.
+    call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate)  
+    if (terminate) then
+        info = USER_STOP
+        exit
+    end if
+    
 end 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/bobyqa/initialize.f90 b/fortran/bobyqa/initialize.f90
index 4f689c8f43..b05ef1ea0d 100644
--- a/fortran/bobyqa/initialize.f90
+++ b/fortran/bobyqa/initialize.f90
@@ -20,7 +20,7 @@ module initialize_bobyqa_mod
 
 
 subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x0, ij, kopt, nf, fhist, fval, &
-    & sl, su, xbase, xhist, xpt, info)
+    & sl, su, xbase, xhist, xpt, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization about the interpolation points & their function values.
 !
@@ -46,9 +46,9 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x0, ij, kopt,
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf, is_finite
-use, non_intrinsic :: infos_mod, only : INFO_DFT
+use, non_intrinsic :: infos_mod, only : INFO_DFT, USER_STOP
 use, non_intrinsic :: message_mod, only : fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : setij
 use, non_intrinsic :: xinbd_mod, only : xinbd
 
@@ -56,6 +56,7 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x0, ij, kopt,
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: ftarget
@@ -90,6 +91,7 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x0, ij, kopt,
 integer(IK) :: npt
 integer(IK) :: subinfo
 logical :: evaluated(size(xpt, 2))
+logical :: terminate
 real(RP) :: f
 real(RP) :: x(size(xpt, 1))
 
@@ -271,6 +273,16 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, xl, xu, x0, ij, kopt,
 kopt = int(minloc(fval, mask=evaluated, dim=1), kind(kopt))
 !!MATLAB: fopt = min(fval(evaluated)); kopt = find(evaluated & ~(fval > fopt), 1, 'first');
 
+! Report the current best value, and check if user asks for early termination.
+! Since this algorithm is a feasible method, the constraint violation is always 0.
+! Secondly, since we only solve bound constrained problems, the nonlinear constraint is set
+! to an empty array.
+terminate = .false.
+call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate)  
+if (terminate) then
+    info = USER_STOP
+end if
+
 !====================!
 !  Calculation ends  !
 !====================!
diff --git a/fortran/cobyla/cobyla.f90 b/fortran/cobyla/cobyla.f90
index e9dca05b5d..af31216731 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
@@ -198,6 +198,10 @@ subroutine cobyla(calcfc, m_nlcon, x, &
 !      be appended to the end of this file if it already exists.
 !   Note that IPRINT = +/-3 can be costly in terms of time and/or space.
 !
+! CALLBACK
+!   Input, function to report progress and request termination.
+!
+!
 ! ETA1, ETA2, GAMMA1, GAMMA2
 !   Input, REAL(RP) scalars, default: ETA1 = 0.1, ETA2 = 0.7, GAMMA1 = 0.5, and GAMMA2 = 2.
 !   ETA1, ETA2, GAMMA1, and GAMMA2 are parameters in the updating scheme of the trust-region radius
@@ -264,7 +268,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, NO_OP_CALLBACK
 use, non_intrinsic :: selectx_mod, only : isbetter
 use, non_intrinsic :: preproc_mod, only : preproc
 use, non_intrinsic :: string_mod, only : num2str
@@ -280,6 +284,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
@@ -316,6 +321,7 @@ subroutine cobyla(calcfc, m_nlcon, x, &
 ! Local variables
 character(len=*), parameter :: solver = 'COBYLA'
 character(len=*), parameter :: srname = 'COBYLA'
+procedure(CALLBACK), pointer :: callback_fcn_loc
 integer(IK) :: info_loc
 integer(IK) :: iprint_loc
 integer(IK) :: m
@@ -587,6 +593,12 @@ subroutine cobyla(calcfc, m_nlcon, x, &
     maxfilt_loc = MAXFILT_DFT
 end if
 
+if (present(callback_fcn)) then
+    callback_fcn_loc => callback_fcn
+else
+    callback_fcn_loc => NO_OP_CALLBACK
+end if
+
 ! Preprocess the inputs in case some of them are invalid. It does nothing if all inputs are valid.
 call preproc(solver, n, iprint_loc, maxfun_loc, maxhist_loc, ftarget_loc, rhobeg_loc, rhoend_loc, &
     & m=m, is_constrained=(m > 0), ctol=ctol_loc, cweight=cweight_loc, eta1=eta1_loc, &
@@ -604,7 +616,7 @@ subroutine cobyla(calcfc, m_nlcon, x, &
 !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, &
     & 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)
+    & f_loc, x, nf_loc, chist_loc, conhist_loc, cstrv_loc, fhist_loc, xhist_loc, callback_fcn_loc, info_loc)
 !--------------------------------------------------------------------------------------------------!
 
 ! 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..6b4a2cd299 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, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine performs the actual calculations of COBYLA.
 !
 ! IPRINT, MAXFILT, MAXFUN, MAXHIST, CTOL, CWEIGHT, ETA1, ETA2, FTARGET, GAMMA1, GAMMA2, RHOBEG,
-! RHOEND, X, NF, F, XHIST, FHIST, CHIST, CONHIST, CSTRV and INFO are identical to the corresponding
+! RHOEND, X, NF, F, XHIST, FHIST, CHIST, CONHIST, CSTRV, INFO and CALLBACK are identical to the corresponding
 ! arguments in subroutine COBYLA.
 !--------------------------------------------------------------------------------------------------!
 
@@ -46,10 +46,10 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist, rangehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING
+use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, DAMAGING_ROUNDING, USER_STOP
 use, non_intrinsic :: linalg_mod, only : inprod, matprod, norm
 use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJCON
+use, non_intrinsic :: pintrf_mod, only : OBJCON, CALLBACK
 use, non_intrinsic :: ratio_mod, only : redrat
 use, non_intrinsic :: redrho_mod, only : redrho
 use, non_intrinsic :: selectx_mod, only : savefilt, selectx, isbetter
@@ -69,6 +69,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: amat(:, :) ! AMAT(N, M_LCON)
 real(RP), intent(in) :: bvec(:) ! BVEC(M_LCON)
+procedure(CALLBACK) :: callback_fcn
 real(RP), intent(in) :: ctol
 real(RP), intent(in) :: cweight
 real(RP), intent(in) :: eta1
@@ -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
@@ -211,7 +213,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 ! from the other vertices to SIM(:, N+1). FVAL, CONMAT, and CVAL hold the function values,
 ! constraint values, and constraint violations on the vertices in the order corresponding to SIM.
 call initxfc(calcfc_internal, iprint, maxfun, constr, ctol, f, ftarget, rhobeg, x, nf, chist, conhist, &
-   & conmat, cval, fhist, fval, sim, simi, xhist, evaluated, subinfo)
+& conmat, cval, fhist, fval, sim, simi, xhist, evaluated, callback_fcn, subinfo)
 
 ! 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
@@ -295,6 +297,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 ! fails but the geometry step is not invoked. Thus the following MAXTR is unlikely to be reached.
 maxtr = max(maxfun, 2_IK * maxfun)  ! MAX: precaution against overflow, which will make 2*MAXFUN < 0.
 info = MAXTR_REACHED
+terminate = .false.
 
 ! Begin the iterative procedure.
 ! After solving a trust-region subproblem, we use three boolean variables to control the workflow.
@@ -304,6 +307,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 ! REDUCE_RHO - Will we reduce rho after the trust-region iteration?
 ! COBYLA never sets IMPROVE_GEO and REDUCE_RHO to TRUE simultaneously.
 do tr = 1, maxtr
+
     ! Increase the penalty parameter CPEN, if needed, so that PREREM = PREREF + CPEN * PREREC > 0.
     ! This is the first (out of two) update of CPEN, where CPEN increases or remains the same.
     ! N.B.: CPEN and the merit function PHI = FVAL + CPEN*CVAL are used at three places only.
@@ -592,6 +596,13 @@ 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.
+    call callback_fcn(sim(:, n+1), fval(n+1), nf, tr, cval(n+1), conmat(:, n+1), terminate)
+    if (terminate) then
+        info = USER_STOP
+        exit
+    end if
+
 end 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/cobyla/initialize.f90 b/fortran/cobyla/initialize.f90
index a3779342f3..9c1384204d 100644
--- a/fortran/cobyla/initialize.f90
+++ b/fortran/cobyla/initialize.f90
@@ -20,7 +20,7 @@ module initialize_cobyla_mod
 
 
 subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x0, nf, chist, &
-    & conhist, conmat, cval, fhist, fval, sim, simi, xhist, evaluated, info)
+    & conhist, conmat, cval, fhist, fval, sim, simi, xhist, evaluated, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization concerning X, function values, and constraints.
 !--------------------------------------------------------------------------------------------------!
@@ -32,15 +32,16 @@ subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x
 use, non_intrinsic :: evaluate_mod, only : evaluate, moderatef, moderatec
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf, is_finite
-use, non_intrinsic :: infos_mod, only : INFO_DFT
+use, non_intrinsic :: infos_mod, only : INFO_DFT, USER_STOP
 use, non_intrinsic :: linalg_mod, only : eye, inv, isinv
 use, non_intrinsic :: message_mod, only : fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJCON
+use, non_intrinsic :: pintrf_mod, only : OBJCON, CALLBACK
 
 implicit none
 
 ! Inputs
 procedure(OBJCON) :: calcfc ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: constr0(:)  ! CONSTR0(M)
@@ -77,6 +78,7 @@ subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x
 integer(IK) :: maxxhist
 integer(IK) :: n
 integer(IK) :: subinfo
+logical :: terminate
 real(RP) :: constr(size(conmat, 1))
 real(RP) :: cstrv
 real(RP) :: f
@@ -195,6 +197,13 @@ subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x
     simi = inv(sim(:, 1:n))
 end if
 
+! Report the current best value, and check if user asks for early termination.
+terminate = .false.
+call callback_fcn(sim(:, n+1), fval(n+1), nf, 0, cval(n+1), conmat(:, n+1), terminate)
+if (terminate) then
+    info = USER_STOP
+end if
+
 !====================!
 !  Calculation ends  !
 !====================!
diff --git a/fortran/common/infos.f90 b/fortran/common/infos.f90
index 72698fd48f..95344a4664 100644
--- a/fortran/common/infos.f90
+++ b/fortran/common/infos.f90
@@ -25,12 +25,12 @@ module infos_mod
 public :: DAMAGING_ROUNDING
 public :: NO_SPACE_BETWEEN_BOUNDS
 public :: ZERO_LINEAR_CONSTRAINT
+public :: USER_STOP
 public :: INVALID_INPUT
 public :: ASSERTION_FAILS
 public :: VALIDATION_FAILS
 public :: MEMORY_ALLOCATION_FAILS
 
-
 integer(IK), parameter :: INFO_DFT = 0
 integer(IK), parameter :: SMALL_TR_RADIUS = 0
 integer(IK), parameter :: FTARGET_ACHIEVED = 1
@@ -43,6 +43,7 @@ module infos_mod
 integer(IK), parameter :: NO_SPACE_BETWEEN_BOUNDS = 6
 integer(IK), parameter :: DAMAGING_ROUNDING = 7
 integer(IK), parameter :: ZERO_LINEAR_CONSTRAINT = 8
+integer(IK), parameter :: USER_STOP = 9
 
 ! Stop-codes.
 ! The following codes are used by ERROR STOP as stop-codes, which should be default integers.
diff --git a/fortran/common/pintrf.f90 b/fortran/common/pintrf.f90
index 8a610e0c60..e2d3c0f4da 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, NO_OP_CALLBACK
 
 
 abstract interface
@@ -36,7 +36,37 @@ subroutine OBJCON(x, f, constr)
     real(RP), intent(out) :: constr(:)
     end subroutine OBJCON
 
+    subroutine CALLBACK(x, f, nf, tr, cstrv, nlconstr, terminate)
+    use consts_mod, only : RP, IK
+    implicit none
+    real(RP), intent(in) :: x(:)
+    real(RP), intent(in) :: f
+    integer(IK), intent(in) :: nf
+    integer(IK), intent(in) :: tr
+    real(RP), intent(in) :: cstrv
+    real(RP), intent(in) :: nlconstr(:)
+    logical, intent(out) :: terminate
+    end subroutine CALLBACK
+
 end interface
 
+contains
+! A dummy subroutine for the callback function.
+subroutine NO_OP_CALLBACK(x, f, nf, tr, cstrv, nlconstr, terminate)
+use consts_mod, only : RP, IK
+implicit none
+real(RP), intent(in) :: x(:)
+real(RP), intent(in) :: f
+integer(IK), intent(in) :: nf
+integer(IK), intent(in) :: tr
+real(RP), intent(in) :: cstrv
+real(RP), intent(in) :: nlconstr(:)
+logical, intent(out) :: terminate
+terminate = .false.
+! We use the variables in a dummy statement to avoid compiler warning about unused variables.
+! The fact that nlconstr is empty for algorithms other than COBYLA doesn't matter since the statement
+! never gets executed.
+if (.false.) terminate = x(1) < 0 .and. f < 0 .and. nf < 0 .and. tr < 0 .and. cstrv < 0 .and. nlconstr(1) < 0
+end subroutine NO_OP_CALLBACK
 
 end module pintrf_mod
diff --git a/fortran/lincoa/initialize.f90 b/fortran/lincoa/initialize.f90
index 09c24fb7b9..32fae2ad13 100644
--- a/fortran/lincoa/initialize.f90
+++ b/fortran/lincoa/initialize.f90
@@ -20,7 +20,7 @@ module initialize_lincoa_mod
 
 
 subroutine initxf(calfun, iprint, maxfun, Aeq, Aineq, amat, beq, bineq, ctol, ftarget, rhobeg, xl, xu,&
-    & x0, b, ij, kopt, nf, chist, cval, fhist, fval, xbase, xhist, xpt, evaluated, info)
+    & x0, b, ij, kopt, nf, chist, cval, fhist, fval, xbase, xhist, xpt, evaluated, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization about the interpolation points & their function values.
 !
@@ -46,17 +46,18 @@ subroutine initxf(calfun, iprint, maxfun, Aeq, Aineq, amat, beq, bineq, ctol, ft
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf, is_finite
-use, non_intrinsic :: infos_mod, only : INFO_DFT
+use, non_intrinsic :: infos_mod, only : INFO_DFT, USER_STOP
 use, non_intrinsic :: linalg_mod, only : matprod, maximum, eye, trueloc
 use, non_intrinsic :: memory_mod, only : safealloc
 use, non_intrinsic :: message_mod, only : fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : setij
 
 implicit none
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: Aeq(:, :)  ! AMAT(Meq, N)
@@ -103,6 +104,7 @@ subroutine initxf(calfun, iprint, maxfun, Aeq, Aineq, amat, beq, bineq, ctol, ft
 integer(IK), allocatable :: ixl(:)
 integer(IK), allocatable :: ixu(:)
 logical :: feasible(size(xpt, 2))
+logical :: terminate
 real(RP) :: constr(count(xl > -REALMAX) + count(xu < REALMAX) + 2 * size(beq) + size(bineq))
 real(RP) :: constr_leq(size(beq))
 real(RP) :: cstrv
@@ -253,6 +255,15 @@ subroutine initxf(calfun, iprint, maxfun, Aeq, Aineq, amat, beq, bineq, ctol, ft
 !!fopt = min(fval(evaluated & feasible));
 !!kopt = find(evaluated & feasible & ~(fval > fopt), 1, 'first');
 
+! Report the current best value, and check if user asks for early termination.
+! Since we only solve linearly constrained problems, the nonlinear constraint is set
+! to an empty array.
+terminate = .false.
+call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, cval(kopt), [real(RP) ::], terminate)
+if (terminate) then
+    info = USER_STOP
+end if
+
 !====================!
 !  Calculation ends  !
 !====================!
diff --git a/fortran/lincoa/lincoa.f90 b/fortran/lincoa/lincoa.f90
index c6aac8a3c7..9fa7916ebb 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
@@ -168,6 +168,10 @@ subroutine lincoa(calfun, x, &
 !      be appended to the end of this file if it already exists.
 !   Note that IPRINT = +/-3 can be costly in terms of time and/or space.
 !
+! CALLBACK
+!   Input, function to report progress and request termination.
+!
+!
 ! ETA1, ETA2, GAMMA1, GAMMA2
 !   Input, REAL(RP) scalars, default: ETA1 = 0.1, ETA2 = 0.7, GAMMA1 = 0.5, and GAMMA2 = 2.
 !   ETA1, ETA2, GAMMA1, and GAMMA2 are parameters in the updating scheme of the trust-region radius
@@ -224,7 +228,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, NO_OP_CALLBACK
 use, non_intrinsic :: preproc_mod, only : preproc
 use, non_intrinsic :: selectx_mod, only : isbetter
 use, non_intrinsic :: string_mod, only : num2str
@@ -239,6 +243,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
@@ -272,6 +277,7 @@ subroutine lincoa(calfun, x, &
 ! Local variables
 character(len=*), parameter :: solver = 'LINCOA'
 character(len=*), parameter :: srname = 'LINCOA'
+procedure(CALLBACK), pointer :: callback_fcn_loc
 integer(IK) :: info_loc
 integer(IK) :: iprint_loc
 integer(IK) :: maxfilt_loc
@@ -494,6 +500,12 @@ subroutine lincoa(calfun, x, &
     maxfilt_loc = MAXFILT_DFT
 end if
 
+if (present(callback_fcn)) then
+    callback_fcn_loc => callback_fcn
+else
+    callback_fcn_loc => NO_OP_CALLBACK
+end if
+
 ! Preprocess the inputs in case some of them are invalid. It does nothing if all inputs are valid.
 call preproc(solver, n, iprint_loc, maxfun_loc, maxhist_loc, ftarget_loc, rhobeg_loc, rhoend_loc, &
     & npt=npt_loc, ctol=ctol_loc, cweight=cweight_loc, eta1=eta1_loc, eta2=eta2_loc, gamma1=gamma1_loc, &
@@ -512,7 +524,7 @@ subroutine lincoa(calfun, x, &
 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)
+    & f_loc, fhist_loc, xhist_loc, callback_fcn_loc, info_loc)
 !--------------------------------------------------------------------------------------------------!
 
 ! 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..5f8b056fd1 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, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine performs the actual calculations of LINCOA.
 !
@@ -78,11 +78,11 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist, rangehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS
+use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP
 use, non_intrinsic :: linalg_mod, only : matprod, maximum, eye, trueloc, linspace, norm, trueloc
 use, non_intrinsic :: memory_mod, only : safealloc
 use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : quadinc, omega_mul, hess_mul, updateh
 use, non_intrinsic :: ratio_mod, only : redrat
 use, non_intrinsic :: redrho_mod, only : redrho
@@ -99,6 +99,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfilt
 integer(IK), intent(in) :: maxfun
@@ -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))
@@ -251,7 +253,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b
 ! Initialize B, XBASE, XPT, FVAL, CVAL, and KOPT, together with the history, NF, IJ, and EVALUATED.
 b = bvec
 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)
+    & x, b, ij, kopt, nf, chist, cval, fhist, fval, xbase, xhist, xpt, evaluated, callback_fcn, subinfo)
 
 ! 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
@@ -342,6 +344,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b
 ! model but the geometry step is not invoked. Thus the following MAXTR is unlikely to be reached.
 maxtr = max(maxfun, 2_IK * maxfun)  ! MAX: precaution against overflow, which will make 2*MAXFUN < 0.
 info = MAXTR_REACHED
+terminate = .false.
 
 ! Begin the iterative procedure.
 ! After solving a trust-region subproblem, we use three boolean variables to control the workflow.
@@ -635,6 +638,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.
+    ! Since we only solve linearly constrained problems, the nonlinear constraint is set
+    ! to an empty array.
+    call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, cval(kopt), [real(RP) ::], terminate)
+    if (terminate) then
+        info = USER_STOP
+        exit
+    end if
+
 end 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/initialize.f90 b/fortran/newuoa/initialize.f90
index c56ec1d382..74a923592d 100644
--- a/fortran/newuoa/initialize.f90
+++ b/fortran/newuoa/initialize.f90
@@ -20,7 +20,7 @@ module initialize_newuoa_mod
 
 
 subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, ij, kopt, nf, fhist, fval, xbase, &
-    & xhist, xpt, info)
+    & xhist, xpt, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization about the interpolation points & their function values.
 !
@@ -49,16 +49,17 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, ij, kopt, nf, fhi
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_finite, is_nan, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT
+use, non_intrinsic :: infos_mod, only : INFO_DFT, USER_STOP
 use, non_intrinsic :: linalg_mod, only : eye
 use, non_intrinsic :: message_mod, only : fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : setij
 
 implicit none
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: ftarget
@@ -87,6 +88,7 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, ij, kopt, nf, fhi
 integer(IK) :: npt
 integer(IK) :: subinfo
 logical :: evaluated(size(fval))
+logical :: terminate
 real(RP) :: f
 real(RP) :: x(size(x0))
 
@@ -224,6 +226,16 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, ij, kopt, nf, fhi
 kopt = int(minloc(fval, mask=evaluated, dim=1), kind(kopt))
 !!MATLAB: fopt = min(fval(evaluated)); kopt = find(evaluated & ~(fval > fopt), 1, 'first')
 
+! Report the current best value, and check if user asks for early termination.
+! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear
+! constraint is set to an empty array.
+terminate = .false.
+call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate)
+if (terminate) then
+    info = USER_STOP
+end if
+
+
 !====================!
 !  Calculation ends  !
 !====================!
diff --git a/fortran/newuoa/newuoa.f90 b/fortran/newuoa/newuoa.f90
index f2c7b2d16f..c50aaf3a94 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
@@ -117,6 +117,9 @@ subroutine newuoa(calfun, x, &
 !      be appended to the end of this file if it already exists.
 !   Note that IPRINT = +/-3 can be costly in terms of time and/or space.
 !
+! CALLBACK
+! Input, function to report progress and request termination.
+!
 ! ETA1, ETA2, GAMMA1, GAMMA2
 !   Input, REAL(RP) scalars, default: ETA1 = 0.1, ETA2 = 0.7, GAMMA1 = 0.5, and GAMMA2 = 2.
 !   ETA1, ETA2, GAMMA1, and GAMMA2 are parameters in the updating scheme of the trust-region radius
@@ -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, NO_OP_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
@@ -205,6 +209,7 @@ subroutine newuoa(calfun, x, &
 ! Local variables
 character(len=*), parameter :: solver = 'NEWUOA'
 character(len=*), parameter :: srname = 'NEWUOA'
+procedure(CALLBACK), pointer :: callback_fcn_loc
 integer(IK) :: info_loc
 integer(IK) :: iprint_loc
 integer(IK) :: maxfun_loc
@@ -321,6 +326,12 @@ subroutine newuoa(calfun, x, &
     maxhist_loc = maxval([maxfun_loc, n + 3_IK, MAXFUN_DIM_DFT * n])
 end if
 
+if (present(callback_fcn)) then
+    callback_fcn_loc => callback_fcn
+else
+    callback_fcn_loc => NO_OP_CALLBACK
+end if
+
 ! Preprocess the inputs in case some of them are invalid.
 call preproc(solver, n, iprint_loc, maxfun_loc, maxhist_loc, ftarget_loc, rhobeg_loc, rhoend_loc, &
     & npt=npt_loc, eta1=eta1_loc, eta2=eta2_loc, gamma1=gamma1_loc, gamma2=gamma2_loc)
@@ -334,7 +345,7 @@ 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)
+    & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, callback_fcn_loc, info_loc)
 !--------------------------------------------------------------------------------------------------!
 
 
diff --git a/fortran/newuoa/newuob.f90 b/fortran/newuoa/newuob.f90
index 742bf787ef..0bdb7d5b39 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, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine performs the actual calculations of NEWUOA.
 !
@@ -57,10 +57,10 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist, rangehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS
+use, non_intrinsic :: infos_mod, only : INFO_DFT, MAXTR_REACHED, SMALL_TR_RADIUS, USER_STOP
 use, non_intrinsic :: linalg_mod, only : norm
 use, non_intrinsic :: message_mod, only : retmsg, rhomsg, fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : quadinc, updateh
 use, non_intrinsic :: ratio_mod, only : redrat
 use, non_intrinsic :: redrho_mod, only : redrho
@@ -76,6 +76,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 integer(IK), intent(in) :: npt
@@ -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))
@@ -173,7 +175,7 @@ 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)
+call initxf(calfun, iprint, maxfun, ftarget, rhobeg, x, ij, kopt, nf, fhist, fval, xbase, xhist, xpt, callback_fcn, subinfo)
 
 ! Initialize X and F according to KOPT.
 x = xbase + xpt(:, kopt)
@@ -243,6 +245,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm
 ! model but the geometry step is not invoked. Thus the following MAXTR is unlikely to be reached.
 maxtr = max(maxfun, 2_IK * maxfun)  ! MAX: precaution against overflow, which will make 2*MAXFUN < 0.
 info = MAXTR_REACHED
+terminate = .false.
 
 ! Begin the iterative procedure.
 ! After solving a trust-region subproblem, we use three boolean variables to control the workflow.
@@ -574,6 +577,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.
+    ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear
+    ! constraint is set to an empty array.
+    call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate)
+    if (terminate) then
+        info = USER_STOP
+        exit
+    end if
+
 end do  ! End of DO TR = 1, MAXTR. The iterative procedure ends.
 
 ! Return from the calculation, after trying the Newton-Raphson step if it has not been tried yet.
diff --git a/fortran/uobyqa/initialize.f90 b/fortran/uobyqa/initialize.f90
index f764ab5436..619ec6b774 100644
--- a/fortran/uobyqa/initialize.f90
+++ b/fortran/uobyqa/initialize.f90
@@ -20,7 +20,7 @@ module initialize_uobyqa_mod
 
 
 subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, kopt, nf, fhist, fval, xbase, &
-    & xhist, xpt, info)
+    & xhist, xpt, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization about the interpolation points & their function values.
 ! See Section 4 of the UOBYQA paper.
@@ -33,15 +33,16 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, kopt, nf, fhist,
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_finite, is_posinf, is_nan
-use, non_intrinsic :: infos_mod, only : INFO_DFT
+use, non_intrinsic :: infos_mod, only : INFO_DFT, USER_STOP
 use, non_intrinsic :: linalg_mod, only : eye, trueloc, linspace
 use, non_intrinsic :: message_mod, only : fmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 
 implicit none
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: ftarget
@@ -72,6 +73,7 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, kopt, nf, fhist,
 integer(IK) :: npt
 integer(IK) :: subinfo
 logical :: evaluated(size(xpt, 2))
+logical :: terminate
 real(RP) :: f
 real(RP) :: x(size(x0))
 real(RP) :: xw(size(x0))
@@ -200,6 +202,15 @@ subroutine initxf(calfun, iprint, maxfun, ftarget, rhobeg, x0, kopt, nf, fhist,
 kopt = int(minloc(fval, mask=evaluated, dim=1), kind(kopt))
 !!MATLAB: fopt = min(fval(evaluated)); kopt = find(evaluated & ~(fval > fopt), 1, 'first')
 
+! Report the current best value, and check if user asks for early termination.
+! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear
+! constraint is set to an empty array.
+terminate = .false.
+call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, 0, 0.0_RP, [real(RP) ::], terminate)
+if (terminate) then
+    info = USER_STOP
+end if
+
 !====================!
 !  Calculation ends  !
 !====================!
diff --git a/fortran/uobyqa/uobyqa.f90 b/fortran/uobyqa/uobyqa.f90
index 55261e3398..f4ad843a93 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
@@ -109,6 +109,10 @@ subroutine uobyqa(calfun, x, &
 !      be appended to the end of this file if it already exists.
 !   Note that IPRINT = +/-3 can be costly in terms of time and/or space.
 !
+! CALLBACK
+!   Input, function to report progress and request termination.
+!
+!
 ! ETA1, ETA2, GAMMA1, GAMMA2
 !   Input, REAL(RP) scalars, default: ETA1 = 0.1, ETA2 = 0.7, GAMMA1 = 0.5, and GAMMA2 = 2.
 !   ETA1, ETA2, GAMMA1, and GAMMA2 are parameters in the updating scheme of the trust-region radius
@@ -161,7 +165,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, NO_OP_CALLBACK
 use, non_intrinsic :: preproc_mod, only : preproc
 use, non_intrinsic :: string_mod, only : num2str
 
@@ -175,6 +179,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
@@ -196,6 +201,7 @@ subroutine uobyqa(calfun, x, &
 ! Local variables
 character(len=*), parameter :: solver = 'UOBYQA'
 character(len=*), parameter :: srname = 'UOBYQA'
+procedure(CALLBACK), pointer :: callback_fcn_loc
 integer(IK) :: info_loc
 integer(IK) :: iprint_loc
 integer(IK) :: maxfun_loc
@@ -304,6 +310,12 @@ subroutine uobyqa(calfun, x, &
     maxhist_loc = maxval([maxfun_loc, 1_IK + (n + 1_IK) * (n + 2_IK) / 2_IK, MAXFUN_DIM_DFT * n])
 end if
 
+if (present(callback_fcn)) then
+    callback_fcn_loc => callback_fcn
+else
+    callback_fcn_loc => NO_OP_CALLBACK
+end if
+
 ! Preprocess the inputs in case some of them are invalid.
 call preproc(solver, n, iprint_loc, maxfun_loc, maxhist_loc, ftarget_loc, rhobeg_loc, rhoend_loc, &
     & eta1=eta1_loc, eta2=eta2_loc, gamma1=gamma1_loc, gamma2=gamma2_loc)
@@ -317,7 +329,7 @@ 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)
+    & gamma2_loc, rhobeg_loc, rhoend_loc, x, nf_loc, f_loc, fhist_loc, xhist_loc, callback_fcn_loc, info_loc)
 !--------------------------------------------------------------------------------------------------!
 
 
diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90
index 0f9c36015c..49a394f0bb 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, callback_fcn, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine performs the major calculations of UOBYQA.
 !
@@ -50,11 +50,11 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 use, non_intrinsic :: evaluate_mod, only : evaluate
 use, non_intrinsic :: history_mod, only : savehist, rangehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
-use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED
+use, non_intrinsic :: infos_mod, only : INFO_DFT, SMALL_TR_RADIUS, MAXTR_REACHED, USER_STOP
 use, non_intrinsic :: linalg_mod, only : vec2smat, smat_mul_vec, norm
 use, non_intrinsic :: memory_mod, only : safealloc
 use, non_intrinsic :: message_mod, only : fmsg, rhomsg, retmsg
-use, non_intrinsic :: pintrf_mod, only : OBJ
+use, non_intrinsic :: pintrf_mod, only : OBJ, CALLBACK
 use, non_intrinsic :: powalg_mod, only : quadinc
 use, non_intrinsic :: ratio_mod, only : redrat
 use, non_intrinsic :: redrho_mod, only : redrho
@@ -70,6 +70,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 
 ! Inputs
 procedure(OBJ) :: calfun  ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+procedure(CALLBACK) :: callback_fcn
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: eta1
@@ -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
@@ -163,7 +165,7 @@ 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)
+call initxf(calfun, iprint, maxfun, ftarget, rhobeg, x, kopt, nf, fhist, fval, xbase, xhist, xpt, callback_fcn, subinfo)
 
 ! Initialize X and F according to KOPT.
 x = xbase + xpt(:, kopt)
@@ -220,6 +222,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r
 ! model but the geometry step is not invoked. Thus the following MAXTR is unlikely to be reached.
 maxtr = max(maxfun, 2_IK * maxfun)  ! MAX: precaution against overflow, which will make 2*MAXFUN < 0.
 info = MAXTR_REACHED
+terminate = .false.
 
 ! Begin the iterative procedure.
 ! After solving a trust-region subproblem, we use three boolean variables to control the workflow.
@@ -458,6 +461,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.
+    ! Since this algorithm is unconstrained, the constraint violation is set to 0 and the nonlinear
+    ! constraint is set to an empty array.
+    call callback_fcn(xbase + xpt(:, kopt), fval(kopt), nf, tr, 0.0_RP, [real(RP) ::], terminate)
+    if (terminate) then
+        info = USER_STOP
+        exit
+    end if
 end do  ! End of DO TR = 1, MAXTR. The iterative procedure ends.
 
 ! Deallocate PL. F2003 automatically deallocate local ALLOCATABLE variables at exit, yet we prefer