diff --git a/fortran/cobyla/cobyla.f90 b/fortran/cobyla/cobyla.f90 index bc8ee2d163..9a94dd81f7 100644 --- a/fortran/cobyla/cobyla.f90 +++ b/fortran/cobyla/cobyla.f90 @@ -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 diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index 6a40027f51..d84756a496 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -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 diff --git a/fortran/cobyla/initialize.f90 b/fortran/cobyla/initialize.f90 index 624a331e79..2618814922 100644 --- a/fortran/cobyla/initialize.f90 +++ b/fortran/cobyla/initialize.f90 @@ -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]) diff --git a/fortran/common/evaluate.f90 b/fortran/common/evaluate.f90 index 0819e1258f..fbc183a4cc 100644 --- a/fortran/common/evaluate.f90 +++ b/fortran/common/evaluate.f90 @@ -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 !====================!