libprima · pawlenartowicz · Sep 2, 2026
@@ -290,15 +290,15 @@ subroutine bobyqa(calfun, x, &
         xl_loc = xl
     end if
 end if
-xl_loc(trueloc(is_nan(xl_loc) .or. xl_loc < -BOUNDMAX)) = -BOUNDMAX
+where (is_nan(xl_loc) .or. xl_loc < -BOUNDMAX) xl_loc = -BOUNDMAX
 
 xu_loc = BOUNDMAX
 if (present(xu)) then
     if (size(xu) > 0) then
         xu_loc = xu
     end if
 end if
-xu_loc(trueloc(is_nan(xu_loc) .or. xu_loc > BOUNDMAX)) = BOUNDMAX
+where (is_nan(xu_loc) .or. xu_loc > BOUNDMAX) xu_loc = BOUNDMAX
 
 ! The solver requires that MINVAL(XU-XL) >= 2*RHOBEG, and we return if MINVAL(XU-XL) < 2*EPS.
 ! It would be better to fix the variables at (XU+XL)/2 wherever XU and XL almost equal, as is done

@@ -708,7 +708,7 @@ function errbd(crvmin, d, gopt, hq, moderr_rec, pq, rho, sl, su, xopt, xpt) resu
 use, non_intrinsic :: consts_mod, only : RP, IK, HALF, DEBUGGING
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_finite
-use, non_intrinsic :: linalg_mod, only : matprod, diag, issymmetric, trueloc
+use, non_intrinsic :: linalg_mod, only : matprod, diag, issymmetric
 use, non_intrinsic :: powalg_mod, only : hess_mul
 
 implicit none
@@ -766,8 +766,8 @@ function errbd(crvmin, d, gopt, hq, moderr_rec, pq, rho, sl, su, xopt, xpt) resu
 xnew = xopt + d
 gnew = gopt + hess_mul(d, xpt, pq, hq)
 bfirst = maxval(abs(moderr_rec))
-bfirst(trueloc(xnew <= sl)) = gnew(trueloc(xnew <= sl)) * rho
-bfirst(trueloc(xnew >= su)) = -gnew(trueloc(xnew >= su)) * rho
+where (xnew <= sl) bfirst = gnew * rho
+where (xnew >= su) bfirst = -gnew * rho
 bsecond = HALF * (diag(hq) + matprod(xpt**2, pq)) * rho**2
 ebound = minval(max(bfirst, bfirst + bsecond))
 if (crvmin > 0) then

@@ -41,7 +41,7 @@ function setdrop_tr(kopt, ximproved, bmat, d, delta, rho, xpt, zmat) result(knew
 use, non_intrinsic :: consts_mod, only : RP, IK, ONE, DEBUGGING
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan, is_finite
-use, non_intrinsic :: linalg_mod, only : issymmetric, trueloc
+use, non_intrinsic :: linalg_mod, only : issymmetric
 use, non_intrinsic :: powalg_mod, only : calden
 
 implicit none
@@ -140,7 +140,7 @@ function setdrop_tr(kopt, ximproved, bmat, d, delta, rho, xpt, zmat) result(knew
 end if
 
 ! SCORE(K) = NaN implies DEN(K) = NaN. We exclude such K as we want DEN to be big.
-score(trueloc(is_nan(score))) = -ONE
+where (is_nan(score)) score = -ONE
 
 knew = 0
 ! The following IF works slightly better than `IF (ANY(SCORE > 0))` from Powell's BOBYQA/LINCOA code.
@@ -198,7 +198,7 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d)
 use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, ONE, TWO, HALF, TEN, EPS, DEBUGGING
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan, is_finite
-use, non_intrinsic :: linalg_mod, only : matprod, inprod, trueloc, norm, issymmetric
+use, non_intrinsic :: linalg_mod, only : matprod, inprod, norm, issymmetric
 use, non_intrinsic :: powalg_mod, only : hess_mul, calden
 
 implicit none
@@ -374,8 +374,8 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d)
 
     ! First, revise SLBD. Note that SLBD_TEST <= 0 unless the input violates XOPT >= SL.
     slbd_test = slbd
-    slbd_test(trueloc(xdiff > 0)) = lfrac(trueloc(xdiff > 0))
-    slbd_test(trueloc(xdiff < 0)) = ufrac(trueloc(xdiff < 0))
+    where (xdiff > 0) slbd_test = lfrac
+    where (xdiff < 0) slbd_test = ufrac
     if (any(slbd_test > slbd)) then
         ilbd = int(maxloc(slbd_test, mask=(.not. is_nan(slbd_test)), dim=1), kind(ilbd))
         slbd = slbd_test(ilbd)
@@ -387,8 +387,8 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d)
 
     ! Second, revise SUBD. Note that SUBD_TEST >= 0 unless the input violates XOPT <= SU.
     subd_test = subd
-    subd_test(trueloc(xdiff > 0)) = ufrac(trueloc(xdiff > 0))
-    subd_test(trueloc(xdiff < 0)) = lfrac(trueloc(xdiff < 0))
+    where (xdiff > 0) subd_test = ufrac
+    where (xdiff < 0) subd_test = lfrac
     if (any(subd_test < subd)) then
         iubd = int(minloc(subd_test, mask=(.not. is_nan(subd_test)), dim=1), kind(iubd))
         subd = max(sumin, subd_test(iubd))
@@ -512,8 +512,8 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d)
     end if
     s = ZERO
     mask_free = (min(xopt - sl, glag) > 0 .or. max(xopt - su, glag) < 0)
-    s(trueloc(mask_free)) = bigstp
-    ggfree = sum(glag(trueloc(mask_free))**2)
+    where (mask_free) s = bigstp
+    ggfree = sum(glag**2, mask=mask_free)
     ! In Powell's code, the subroutine returns immediately if GGFREE is 0. However, GGFREE depends
     ! on GLAG, which in turn depends on UPHILL. It can happen that GGFREE is 0 when UPHILL = 0 but
     ! not so when UPHILL= 1. Thus we skip the iteration for the current UPHILL but do not return.
@@ -538,22 +538,22 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d)
         mask_fixl = (s >= bigstp .and. xtemp <= sl)  ! S == BIGSTP & XTEMP == SL
         mask_fixu = (s >= bigstp .and. xtemp >= su)  ! S == BIGSTP & XTEMP == SU
         mask_free = (s >= bigstp .and. .not. (mask_fixl .or. mask_fixu))
-        s(trueloc(mask_fixl)) = sl(trueloc(mask_fixl)) - xopt(trueloc(mask_fixl))
-        s(trueloc(mask_fixu)) = su(trueloc(mask_fixu)) - xopt(trueloc(mask_fixu))
-        sfixsq = sfixsq + sum(s(trueloc(mask_fixl .or. mask_fixu))**2)
-        ggfree = sum(glag(trueloc(mask_free))**2)
+        where (mask_fixl) s = sl - xopt
+        where (mask_fixu) s = su - xopt
+        sfixsq = sfixsq + sum(s**2, mask=(mask_fixl .or. mask_fixu))
+        ggfree = sum(glag**2, mask=mask_free)
         if (.not. (sfixsq > ssqsav .and. ggfree > 0)) then
             exit
         end if
     end do
 
     ! Set the remaining free components of S and all components of XCAUCHY. S may be scaled later.
-    x(trueloc(glag > 0)) = sl(trueloc(glag > 0))
-    x(trueloc(glag <= 0)) = su(trueloc(glag <= 0))
-    x(trueloc(abs(s) <= 0)) = xopt(trueloc(abs(s) <= 0))
+    where (glag > 0) x = sl
+    where (glag <= 0) x = su
+    where (abs(s) <= 0) x = xopt
     xtemp = max(sl, min(su, xopt - grdstp * glag))
-    x(trueloc(s >= bigstp)) = xtemp(trueloc(s >= bigstp))  ! S == BIGSTP
-    s(trueloc(s >= bigstp)) = -grdstp * glag(trueloc(s >= bigstp))  ! S == BIGSTP
+    where (s >= bigstp) x = xtemp  ! S == BIGSTP
+    where (s >= bigstp) s = -grdstp * glag  ! S == BIGSTP
     gs = inprod(glag, s)
 
     ! Set CURV to the curvature of the KNEW-th Lagrange function along S. Scale S by a factor less

@@ -98,7 +98,7 @@ subroutine rescue(calfun, solver, iprint, maxfun, delta, ftarget, xl, xu, kopt,
 use, non_intrinsic :: history_mod, only : savehist
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf, is_finite
 use, non_intrinsic :: infos_mod, only : MAXFUN_REACHED, INFO_DFT
-use, non_intrinsic :: linalg_mod, only : issymmetric, matprod, inprod, r1update, r2update, trueloc
+use, non_intrinsic :: linalg_mod, only : issymmetric, matprod, inprod, r1update, r2update
 use, non_intrinsic :: message_mod, only : fmsg
 use, non_intrinsic :: pintrf_mod, only : OBJ
 use, non_intrinsic :: powalg_mod, only : hess_mul, setij
@@ -255,9 +255,12 @@ subroutine rescue(calfun, solver, iprint, maxfun, delta, ftarget, xl, xu, kopt,
 ptsaux(1, :) = min(delta, su)
 ptsaux(2, :) = max(-delta, sl)
 mask = (ptsaux(1, :) + ptsaux(2, :) < 0)
-ptsaux([1, 2], trueloc(mask)) = ptsaux([2, 1], trueloc(mask))
+where (mask)
+    ptsaux(1, :) = max(-delta, sl)
+    ptsaux(2, :) = min(delta, su)
+end where
 mask = (abs(ptsaux(2, :)) < HALF * abs(ptsaux(1, :)))
-ptsaux(2, trueloc(mask)) = HALF * ptsaux(1, trueloc(mask))
+where (mask) ptsaux(2, :) = HALF * ptsaux(1, :)
 
 ! Set the identifiers of the artificial interpolation points that are along a coordinate direction
 ! from XOPT, and set the corresponding nonzero elements of BMAT and ZMAT.
@@ -394,8 +397,8 @@ subroutine rescue(calfun, solver, iprint, maxfun, delta, ftarget, xl, xu, kopt,
     ! For all K with PTSID(K) > 0, calculate the denominator DEN(K) = SIGMA in the updating formula
     ! of H for XPT(:, KORIG) to replace XPT_PROV(:, K).
     den = ZERO
-    hdiag(trueloc(ptsid > 0)) = sum(zmat(trueloc(ptsid > 0), :)**2, dim=2)
-    den(trueloc(ptsid > 0)) = hdiag(trueloc(ptsid > 0)) * beta + vlag(trueloc(ptsid > 0))**2
+    where (ptsid > 0) hdiag = sum(zmat**2, dim=2)
+    where (ptsid > 0) den = hdiag * beta + vlag(1:npt)**2
 
     ! Attempt setting KPROV to the index of the provisional point to be replaced with the KORIG-th
     ! original interpolation point. We choose KPROV by maximizing DEN(KPROV), which will be the
@@ -551,7 +554,7 @@ subroutine rescue(calfun, solver, iprint, maxfun, delta, ftarget, xl, xu, kopt,
         moderr = f - vquad
         gopt = gopt + moderr * bmat(:, kpt)
         pqinc = moderr * matprod(zmat, zmat(kpt, :))
-        pq(trueloc(ptsid <= 0)) = pq(trueloc(ptsid <= 0)) + pqinc(trueloc(ptsid <= 0))
+        where (ptsid <= 0) pq = pq + pqinc
         do k = 1, npt
             if (ptsid(k) <= 0) then
                 cycle

@@ -57,7 +57,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     & DEBUGGING
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan, is_finite
-use, non_intrinsic :: linalg_mod, only : inprod, issymmetric, trueloc, norm
+use, non_intrinsic :: linalg_mod, only : issymmetric, norm
 use, non_intrinsic :: powalg_mod, only : hess_mul
 use, non_intrinsic :: univar_mod, only : interval_max
 
@@ -189,8 +189,8 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
 ! value at XOPT. Thus XBDI(I) shows whether or not to fix the I-th variable at one of its bounds
 ! initially, with NACT being set to the number of fixed variables.
 xbdi = 0
-xbdi(trueloc(xopt >= su .and. gopt <= 0)) = 1
-xbdi(trueloc(xopt <= sl .and. gopt >= 0)) = -1
+where (xopt >= su .and. gopt <= 0) xbdi = 1
+where (xopt <= sl .and. gopt >= 0) xbdi = -1
 nact = int(count(xbdi /= 0), kind(nact))
 
 ! Initialized D and CRVMIN.
@@ -199,7 +199,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
 
 ! GNEW is the gradient at the current iterate.
 gnew = gopt
-gredsq = sum(gnew(trueloc(xbdi == 0))**2)
+gredsq = sum(gnew**2, mask=(xbdi == 0))
 ! DELSQ is the upper bound on the sum of squares of the free variables.
 delsq = delta * delta
 ! QRED is the reduction in Q so far.
@@ -221,7 +221,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
 ! approximate solution to the subproblem (1.8), even if there are hundreds of variables."
 maxiter = int(min(10**min(4, range(0_IK)), int(n - nact)**2), IK)
 do iter = 1, maxiter
-    resid = delsq - sum(d(trueloc(xbdi == 0))**2)
+    resid = delsq - sum(d**2, mask=(xbdi == 0))
     if (resid <= 0) then
         twod_search = .true.
         exit
@@ -237,9 +237,9 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     else
         s = beta * s - gnew
     end if
-    s(trueloc(xbdi /= 0)) = ZERO
+    where (xbdi /= 0) s = ZERO
     stepsq = sum(s**2)
-    ds = inprod(d(trueloc(xbdi == 0)), s(trueloc(xbdi == 0)))
+    ds = sum(d * s, mask=(xbdi == 0))
 
     if (.not. (stepsq > EPS * delsq .and. gredsq * delsq > (tol * qred)**2 .and. .not. is_nan(ds))) then
         exit
@@ -270,7 +270,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     end if
 
     hs = hess_mul(s, xpt, pq, hq)
-    shs = inprod(s(trueloc(xbdi == 0)), hs(trueloc(xbdi == 0)))
+    shs = sum(s * hs, mask=(xbdi == 0))
     stplen = bstep
     if (shs > 0) then
         stplen = min(bstep, gredsq / shs)
@@ -316,7 +316,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     !where (s > 0) sbound = min(stplen * s, su - xnew) / s
     !where (s < 0) sbound = max(stplen * s, sl - xnew) / s
     !----------------------------------------------------------------------------------------------!
-    sbound(trueloc(is_nan(sbound))) = stplen  ! Needed? No if we are sure that D and S are finite.
+    where (is_nan(sbound)) sbound = stplen  ! Needed? No if we are sure that D and S are finite.
     iact = 0
     if (any(sbound < stplen)) then
         iact = int(minloc(sbound, dim=1), kind(iact))
@@ -348,7 +348,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
         end if
         ggsav = gredsq
         gnew = gnew + stplen * hs
-        gredsq = sum(gnew(trueloc(xbdi == 0))**2)
+        gredsq = sum(gnew**2, mask=(xbdi == 0))
         dold = d
         d = d + stplen * s
 
@@ -380,7 +380,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
         end if
         beta = ZERO
         itercg = 0
-        gredsq = sum(gnew(trueloc(xbdi == 0))**2)
+        gredsq = sum(gnew**2, mask=(xbdi == 0))
     elseif (stplen < bstep) then
         ! Either apply another conjugate gradient iteration or exit.
         ! N.B. ITERCG > N - NACT is impossible.
@@ -426,20 +426,20 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     xnew = xopt + d
 
     ! Update XBDI. It indicates whether the lower (-1) or upper bound (+1) is reached or not (0).
-    xbdi(trueloc(xbdi == 0 .and. (xnew >= su))) = 1
-    xbdi(trueloc(xbdi == 0 .and. (xnew <= sl))) = -1
+    where (xbdi == 0 .and. (xnew >= su)) xbdi = 1
+    where (xbdi == 0 .and. (xnew <= sl)) xbdi = -1
     nact = int(count(xbdi /= 0), kind(nact))
     if (nact >= n - 1) then
         exit
     end if
 
     ! Update GREDSQ, DREDG, DREDSQ.
-    gredsq = sum(gnew(trueloc(xbdi == 0))**2)
-    dredg = inprod(d(trueloc(xbdi == 0)), gnew(trueloc(xbdi == 0)))
+    gredsq = sum(gnew**2, mask=(xbdi == 0))
+    dredg = sum(d * gnew, mask=(xbdi == 0))
     if (iter == 1 .or. nact > nactsav) then
-        dredsq = sum(d(trueloc(xbdi == 0))**2) ! In theory, DREDSQ changes only when NACT increases.
+        dredsq = sum(d**2, mask=(xbdi == 0)) ! In theory, DREDSQ changes only when NACT increases.
         dred = d
-        dred(trueloc(xbdi /= 0)) = ZERO
+        where (xbdi /= 0) dred = ZERO
         hdred = hess_mul(dred, xpt, pq, hq)
         nactsav = nact
     end if
@@ -452,7 +452,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     end if
     temp = sqrt(temp)
     s = (dredg * d - dredsq * gnew) / temp
-    s(trueloc(xbdi /= 0)) = ZERO
+    where (xbdi /= 0) s = ZERO
     sredg = -temp
 
     ! By considering the simple bounds on the free variables, calculate an upper bound on the
@@ -479,7 +479,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     sqdscr = -REALMAX
     where (xbdi == 0 .and. su - xopt < sqrt(ssq)) sqdscr = sqrt(max(ZERO, ssq - (su - xopt)**2))
     where (sqdscr + s > 0) tanbd = min(tanbd, (su - xnew) / (sqdscr + s))
-    tanbd(trueloc(is_nan(tanbd))) = ZERO
+    where (is_nan(tanbd)) tanbd = ZERO
     !----------------------------------------------------------------------------------------------!
     !!MATLAB code for defining TANBD:
     !!xfree = (xbdi == 0);
@@ -509,9 +509,9 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
 
     ! Calculate HS and some curvatures for the alternative iteration.
     hs = hess_mul(s, xpt, pq, hq)
-    shs = inprod(s(trueloc(xbdi == 0)), hs(trueloc(xbdi == 0)))
-    dhs = inprod(d(trueloc(xbdi == 0)), hs(trueloc(xbdi == 0)))
-    dhd = inprod(d(trueloc(xbdi == 0)), hdred(trueloc(xbdi == 0)))
+    shs = sum(s * hs, mask=(xbdi == 0))
+    dhs = sum(d * hs, mask=(xbdi == 0))
+    dhd = sum(d * hdred, mask=(xbdi == 0))
 
     ! Seek the greatest reduction in Q for a range of equally spaced values of HANGT in [0, ANGBD],
     ! with HANGT being the TANGENT of HALF the angle of the alternative iteration.
@@ -538,7 +538,7 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
     sth = min((hangt + hangt) / (ONE + hangt**2), hangt + hangt)
     gnew = gnew + (cth - ONE) * hdred + sth * hs
     dold = d
-    d(trueloc(xbdi == 0)) = cth * d(trueloc(xbdi == 0)) + sth * s(trueloc(xbdi == 0))
+    where (xbdi == 0) d = cth * d + sth * s
 
     ! Exit in case of Inf/NaN in D.
     if (.not. is_finite(sum(abs(d)))) then
@@ -558,8 +558,8 @@ subroutine trsbox(delta, gopt_in, hq_in, pq_in, sl, su, tol, xopt, xpt, crvmin,
 
 ! Set D, giving careful attention to the bounds.
 xnew = max(sl, min(su, xopt + d))
-xnew(trueloc(xbdi == -1)) = sl(trueloc(xbdi == -1))
-xnew(trueloc(xbdi == 1)) = su(trueloc(xbdi == 1))
+where (xbdi == -1) xnew = sl
+where (xbdi == 1) xnew = su
 d = xnew - xopt
 
 ! Set CRVMIN to ZERO if it has never been set or becomes NaN due to ill conditioning.

@@ -390,7 +390,7 @@ subroutine tryqalt(bmat, fval, ratio, sl, su, xopt, xpt, zmat, itest, gopt, hq,
 use, non_intrinsic :: consts_mod, only : RP, IK, ZERO, TEN, TENTH, DEBUGGING
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan, is_posinf
-use, non_intrinsic :: linalg_mod, only : matprod, inprod, issymmetric, trueloc
+use, non_intrinsic :: linalg_mod, only : matprod, inprod, issymmetric
 use, non_intrinsic :: powalg_mod, only : hess_mul
 
 implicit none
@@ -456,17 +456,17 @@ subroutine tryqalt(bmat, fval, ratio, sl, su, xopt, xpt, zmat, itest, gopt, hq,
 
 ! Calculate the norm square of the projected gradient.
 pgopt = gopt
-pgopt(trueloc(xopt >= su)) = max(ZERO, gopt(trueloc(xopt >= su)))
-pgopt(trueloc(xopt <= sl)) = min(ZERO, gopt(trueloc(xopt <= sl)))
+where (xopt >= su) pgopt = max(ZERO, gopt)
+where (xopt <= sl) pgopt = min(ZERO, gopt)
 
 ! Calculate the parameters of the least Frobenius norm interpolant to the current data.
 pqalt = matprod(zmat, matprod(fval, zmat))
 galt = matprod(bmat(:, 1:npt), fval) + hess_mul(xopt, xpt, pqalt)
 
 ! Calculate the norm square of the projected alternative gradient.
 pgalt = galt
-pgalt(trueloc(xopt >= su)) = max(ZERO, galt(trueloc(xopt >= su)))
-pgalt(trueloc(xopt <= sl)) = min(ZERO, galt(trueloc(xopt <= sl)))
+where (xopt >= su) pgalt = max(ZERO, galt)
+where (xopt <= sl) pgalt = min(ZERO, galt)
 
 ! Test whether to replace the new quadratic model by the least Frobenius norm interpolant,
 ! making the replacement if the test is satisfied.