diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index 18e24edace..a01410745b 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -17,7 +17,7 @@ module cobylb_mod ! ! Started: July 2021 ! -! Last Modified: Tuesday, March 19, 2024 PM08:44:31 +! Last Modified: Tuesday, March 19, 2024 PM10:40:31 !--------------------------------------------------------------------------------------------------! implicit none @@ -163,8 +163,8 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! 1. FACTOR_ALPHA < FACTOR_GAMMA < 1 < FACTOR_BETA. ! 2. FACTOR_GAMMA has nothing to do with GAMMA1 and GAMMA2, which are the contracting/expanding ! factors for updating the trust-region radius DELTA. -! 3. Powell's used one more factor FACTOR_DELTA = 1.1 (in general, 1 < FACTOR_DELTA <= FACTOR_BETA). -! It has nothing to do with DELTA, which is the trust-region radius. It was used when defining +! 3. Powell used one more factor FACTOR_DELTA = 1.1 (in general, 1 < FACTOR_DELTA <= FACTOR_BETA). +! It had nothing to do with DELTA, which is the trust-region radius. It was used when defining ! JDROP_TR. We use a completely different scheme (see SETDROP_TR), which does not need FACTOR_DELTA. real(RP), parameter :: factor_alpha = QUART ! The factor alpha in the COBYLA paper real(RP), parameter :: factor_beta = 2.1_RP ! The factor beta in the COBYLA paper @@ -579,9 +579,9 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! objective and constraints at X, assuming them to have the values at the closest point. ! N.B.: ! 1. If this happens, do NOT include X into the filter, as F and CONSTR are inaccurate. - ! 2. In precise arithmetic, X - SIM(:, N+1) = D has a length of FACTOR_ALPHA * DELTA > 0. - ! Due to rounding, however, X may be quite close to SIM(:, N+1). In experiments with single - ! precision on 20240317, X = SIM(:, N+1) did happen. + ! 2. In precise arithmetic, the geometry improving step ensures that the distance between X + ! and any interpolation point is at least FACTOR_GAMMA*DELTA, yet X may be close to them due + ! to rounding. In an experiment with single precision on 20240317, X = SIM(:, N+1) occurred. distsq(n + 1) = sum((x - sim(:, n + 1))**2) distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)] ! Implied do-loop !!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1) % Implicit expansion @@ -681,7 +681,6 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! Print a return message according to IPRINT. call retmsg(solver, info, iprint, nf, f, x, cstrv, constr) - !====================! ! Calculation ends ! !====================! diff --git a/fortran/common/history.f90 b/fortran/common/history.f90 index 54c282582a..c050e1d2e3 100644 --- a/fortran/common/history.f90 +++ b/fortran/common/history.f90 @@ -7,7 +7,7 @@ module history_mod ! ! Started: July 2020 ! -! Last Modified: Tuesday, March 19, 2024 PM06:52:07 +! Last Modified: Tuesday, March 19, 2024 PM10:31:28 !--------------------------------------------------------------------------------------------------! implicit none @@ -287,19 +287,21 @@ subroutine savehist(nf, x, xhist, f, fhist, cstrv, chist, constr, conhist) ! The following code checks that XHIST does not contain a segment that repeats. If such a segment ! is found, we believe that the solver has encountered an infinite cycle, which would be a bug. - ! Note that we check this only if NF > (N+1)*(N+2)/2, so that the initialization has (surely) - ! finished. This is because XHIST may contain repeating segments during the initialization if - ! X0 + RHOBEG == X0, which can happen if all entries of X0 are excessively large compared with - ! RHOBEG. It is possible to revise the initialization subroutine to avoid repetition, but we - ! choose not to, a motivation being to keep the initialization parallelizable. + ! N.B.: + ! 1. We check this only if NF > (N+1)*(N+2)/2, when the initialization has surely finished. This + ! is because XHIST may contain repeating segments during the initialization if X0 + RHOBEG = X0, + ! which can happen if all entries of X0 are excessively large compared with RHOBEG. It is + ! possible to revise the initialization subroutine to avoid repetition, but we choose not to, + ! a motivation being to keep the initialization parallelizable. + ! 2. We skip the test if N = 1, as false positive may occur (also possible when N > 1, but rare). nhist = min(nf, maxxhist) n = int(size(x), kind(n)) - do i = 1, 10 - if (nhist >= 2 * i .and. nf > (n + 1) * (n + 2) / 2) then + if (n > 1 .and. nf > (n + 1) * (n + 2) / 2) then + do i = 1, min(10_IK, nhist / 2_IK) call wassert(.not. all(abs(xhist(:, nhist - i + 1:nhist) - xhist(:, nhist - 2 * i + 1:nhist - i)) <= 0), & & 'XHIST does not contain a repeating segment of length '//num2str(i), srname) - end if - end do + end do + end if end if end subroutine savehist