libprima · nbelakovski · Aug 25, 2024
@@ -499,13 +499,14 @@ subroutine cobyla(calcfc, m_nlcon, x, &
 ! If NLCONSTR0 is present, then F0 must be present, and we assume that F(X0) = F0 even if F0 is NaN.
 ! If NLCONSTR0 is absent, then F0 must be either absent or NaN, both of which will be interpreted as
 ! F(X0) is not provided and we have to evaluate F(X0) and NLCONSTR(X0) now.
-constr_loc(1:m - m_nlcon) = moderatec(matprod(x, amat) - bvec)
 if (present(f0) .and. present(nlconstr0) .and. all(is_finite(x))) then
     f_loc = moderatef(f0)
+    constr_loc(1:m - m_nlcon) = moderatec(matprod(x, amat) - bvec)
     constr_loc(m - m_nlcon + 1:m) = moderatec(nlconstr0)
 else
     x = moderatex(x)
-    call evaluate(calcfc, x, f_loc, constr_loc(m - m_nlcon + 1:m))
+    call evaluate(calcfc, x, f_loc, constr_loc, amat, bvec)
+    constr_loc(1:m - m_nlcon) = moderatec(constr_loc(1:m - m_nlcon))
     ! N.B.: Do NOT call FMSG, SAVEHIST, or SAVEFILT for the function/constraint evaluation at X0.
     ! They will be called during the initialization, which will read the function/constraint at X0.
 end if

@@ -202,7 +202,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 ! function value (regardless of the constraint violation), and SIM(:, 1:N) holds the displacements
 ! 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, &
+call initxfc(calcfc, iprint, maxfun, constr, amat, bvec, 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.
@@ -393,7 +393,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
             cstrv = cval(j)
         else
             ! Evaluate the objective and constraints at X, taking care of possible Inf/NaN values.
-            call evaluate(calcfc_internal, x, f, constr)
+            call evaluate(calcfc, x, f, constr, amat, bvec)
             cstrv = maximum([ZERO, constr])
             nf = nf + 1_IK
             ! Save X, F, CONSTR, CSTRV into the history.
@@ -585,7 +585,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
             cstrv = cval(j)
         else
             ! Evaluate the objective and constraints at X, taking care of possible Inf/NaN values.
-            call evaluate(calcfc_internal, x, f, constr)
+            call evaluate(calcfc, x, f, constr, amat, bvec)
             cstrv = maximum([ZERO, constr])
             nf = nf + 1_IK
             ! Save X, F, CONSTR, CSTRV into the history.
@@ -650,7 +650,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
 ! Ensure that D has not been updated after SHORTD == TRUE occurred, or the code below is incorrect.
 x = sim(:, n + 1) + d
 if (info == SMALL_TR_RADIUS .and. shortd .and. norm(x - sim(:, n + 1)) > 1.0E-3_RP * rhoend .and. nf < maxfun) then
-    call evaluate(calcfc_internal, x, f, constr)
+    call evaluate(calcfc, x, f, constr, amat, bvec)
     cstrv = maximum([ZERO, constr])
     nf = nf + 1_IK
     ! Save X, F, CONSTR, CSTRV into the history.
@@ -702,25 +702,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et
         & 'No point in the history is better than X', srname)
 end if
 
-
-contains
-
-
-subroutine calcfc_internal(x_internal, f_internal, constr_internal)
-!--------------------------------------------------------------------------------------------------!
-! This internal subroutine evaluates the objective function and ALL the constraints.
-! In MATLAB/Python/R/Julia, this can be implemented as a lambda function / anonymous function.
-!--------------------------------------------------------------------------------------------------!b
-implicit none
-! Inputs
-real(RP), intent(in) :: x_internal(:)
-! Outputs
-real(RP), intent(out) :: f_internal
-real(RP), intent(out) :: constr_internal(:)
-constr_internal(1:m_lcon) = matprod(x_internal, amat) - bvec
-call calcfc(x_internal, f_internal, constr_internal(m_lcon + 1:m))
-end subroutine calcfc_internal
-
 end subroutine cobylb
 
 

@@ -19,7 +19,7 @@ module initialize_cobyla_mod
 contains
 
 
-subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x0, nf, chist, &
+subroutine initxfc(calcfc, iprint, maxfun, constr0, amat, bvec, ctol, f0, ftarget, rhobeg, x0, nf, chist, &
     & conhist, conmat, cval, fhist, fval, sim, simi, xhist, evaluated, info)
 !--------------------------------------------------------------------------------------------------!
 ! This subroutine does the initialization concerning X, function values, and constraints.
@@ -44,6 +44,8 @@ subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x
 integer(IK), intent(in) :: iprint
 integer(IK), intent(in) :: maxfun
 real(RP), intent(in) :: constr0(:)  ! CONSTR0(M)
+real(RP), intent(in) :: amat(:, :)
+real(RP), intent(in) :: bvec(:)
 real(RP), intent(in) :: ctol
 real(RP), intent(in) :: f0
 real(RP), intent(in) :: ftarget
@@ -156,7 +158,7 @@ subroutine initxfc(calcfc, iprint, maxfun, constr0, ctol, f0, ftarget, rhobeg, x
     else
         j = k - 1_IK
         x(j) = x(j) + rhobeg
-        call evaluate(calcfc, x, f, constr)
+        call evaluate(calcfc, x, f, constr, amat, bvec)
     end if
     cstrv = maximum([ZERO, constr])
 

@@ -148,20 +148,23 @@ subroutine evaluatef(calfun, x, f)
 end subroutine evaluatef
 
 
-subroutine evaluatefc(calcfc, x, f, constr)
+subroutine evaluatefc(calcfc, x, f, constr, amat, bvec)
 !--------------------------------------------------------------------------------------------------!
 ! This function evaluates CALCFC at X, setting F to the objective function value and CONSTR to the
 ! constraint value. Nan/Inf are handled by a moderated extreme barrier.
 !--------------------------------------------------------------------------------------------------!
 ! Common modules
-use, non_intrinsic :: consts_mod, only : RP, DEBUGGING
+use, non_intrinsic :: consts_mod, only : DEBUGGING, IK, RP
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
+use, non_intrinsic :: linalg_mod, only : matprod
 use, non_intrinsic :: pintrf_mod, only : OBJCON
 implicit none
 
 ! Inputs
 procedure(OBJCON) :: calcfc ! N.B.: INTENT cannot be specified if a dummy procedure is not a POINTER
+real(RP), intent(in) :: amat(:, :)
+real(RP), intent(in) :: bvec(:)
 real(RP), intent(in) :: x(:)
 
 ! Outputs
@@ -170,6 +173,13 @@ subroutine evaluatefc(calcfc, x, f, constr)
 
 ! Local variables
 character(len=*), parameter :: srname = 'EVALUATEFC'
+integer(IK) :: m
+integer(IK) :: m_lcon
+
+! Sizes
+m = int(size(constr), kind(m))
+m_lcon = int(size(bvec), kind(m_lcon))
+
 
 ! Preconditions
 if (DEBUGGING) then
@@ -182,18 +192,20 @@ subroutine evaluatefc(calcfc, x, f, constr)
 ! Calculation starts !
 !====================!
 
+constr(1:m_lcon) = matprod(x, amat) - bvec
+
 if (any(is_nan(x))) then
     ! Although this should not happen unless there is a bug, we include this case for robustness.
     ! Set F, CONSTR, and CSTRV to NaN.
     f = sum(x)
-    constr = f
+    constr(m_lcon + 1:m) = f
 else
-    call calcfc(moderatex(x), f, constr)  ! Evaluate F and CONSTR; We moderate X before doing so.
+    call calcfc(moderatex(x), f, constr(m_lcon + 1:m))  ! Evaluate F and CONSTR; We moderate X before doing so.
 
     ! Moderated extreme barrier: replace NaN/huge objective or constraint values with a large but
     ! finite value. This is naive, and better approaches surely exist.
     f = moderatef(f)
-    constr = moderatec(constr)
+    constr(m_lcon + 1:m) = moderatec(constr(m_lcon + 1:m))
 end if
 
 !====================!