240319.224110.HKT fortran: revise the checking of reapeating segement…

…s in `xhist`

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  !
 !====================!

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