From 8d044d159a1eeac9f6d89405713f7efa78a46484 Mon Sep 17 00:00:00 2001 From: Nickolai Belakovski Date: Sun, 25 Aug 2024 11:34:13 -0400 Subject: [PATCH] Refactor calcfc_internal to remove nested function This is in support of efforts to integrate PRIMA with SciPy. They provide wheels for musl, and in musl you cannot have nested functions in Fortran (at least, not without gcc 14, which is cutting edge, and therefore it's not practical for scipy to switch to it). --- fortran/cobyla/cobyla.f90 | 5 +++-- fortran/cobyla/cobylb.f90 | 27 ++++----------------------- fortran/cobyla/initialize.f90 | 6 ++++-- fortran/common/evaluate.f90 | 22 +++++++++++++++++----- 4 files changed, 28 insertions(+), 32 deletions(-) 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 !====================!