subroutine trstlp (n, m, a, b, rho, dx, ifull, iact, z, zdota, vmultc, sdirn, dxnew, &
vmultd)
implicit real (wp) (a-h, o-z)
dimension a (n,*), b (*), dx (*), iact (*), z (n,*), zdota (*), vmultc (*), &
sdirn (*), dxnew (*), vmultd (*)
!
! This subroutine calculates an N-component vector DX by applying the
! following two stages. In the first stage, DX is set to the shortest
! vector that minimizes the greatest violation of the constraints
! A(1,K)*DX(1)+A(2,K)*DX(2)+...+A(N,K)*DX(N) .GE. B(K), K=2,3,...,M,
! subject to the Euclidean length of DX being at most RHO. If its length is
! strictly less than RHO, then we use the resultant freedom in DX to
! minimize the objective function
! -A(1,M+1)*DX(1)-A(2,M+1)*DX(2)-...-A(N,M+1)*DX(N)
! subject to no increase in any greatest constraint violation. This
! notation allows the gradient of the objective function to be regarded as
! the gradient of a constraint. Therefore the two stages are distinguished
! by MCON .EQ. M and MCON .GT. M respectively. It is possible that a
! degeneracy may prevent DX from attaining the target length RHO. Then the
! value IFULL=0 would be set, but usually IFULL=1 on return.
!
! In general NACT is the number of constraints in the active set and
! IACT(1),...,IACT(NACT) are their indices, while the remainder of IACT
! contains a permutation of the remaining constraint indices. Further, Z is
! an orthogonal matrix whose first NACT columns can be regarded as the
! result of Gram-Schmidt applied to the active constraint gradients. For
! J=1,2,...,NACT, the number ZDOTA(J) is the scalar product of the J-th
! column of Z with the gradient of the J-th active constraint. DX is the
! current vector of variables and here the residuals of the active
! constraints should be zero. Further, the active constraints have
! nonnegative Lagrange multipliers that are held at the beginning of
! VMULTC. The remainder of this vector holds the residuals of the inactive
! constraints at DX, the ordering of the components of VMULTC being in
! agreement with the permutation of the indices of the constraints that is
! in IACT. All these residuals are nonnegative, which is achieved by the
! shift RESMAX that makes the least residual zero.
!
! Initialize Z and some other variables. The value of RESMAX will be
! appropriate to DX=0, while ICON will be the index of a most violated
! constraint if RESMAX is positive. Usually during the first stage the
! vector SDIRN gives a search direction that reduces all the active
! constraint violations by one simultaneously.
!
ifull = 1
mcon = m
nact = 0
resmax = 0.0_wp
do i = 1, n
do j = 1, n
z (i, j) = 0.0_wp
end do
z (i, i) = 1.0_wp
dx (i) = 0.0_wp
end do
if (m >= 1) then
do k = 1, m
if (b(k) > resmax) then
resmax = b (k)
icon = k
end if
end do
do k = 1, m
iact (k) = k
vmultc (k) = resmax - b (k)
end do
end if
if (resmax == 0.0_wp) go to 480
do i = 1, n
sdirn (i) = 0.0_wp
end do
!
! End the current stage of the calculation if 3 consecutive iterations
! have either failed to reduce the best calculated value of the objective
! function or to increase the number of active constraints since the best
! value was calculated. This strategy prevents cycling, but there is a
! remote possibility that it will cause premature termination.
!
60 optold = 0.0_wp
icount = 0
70 if (mcon == m) then
optnew = resmax
else
optnew = 0.0_wp
do i = 1, n
optnew = optnew - dx (i) * a (i, mcon)
end do
end if
if (icount == 0 .or. optnew < optold) then
optold = optnew
nactx = nact
icount = 3
else if (nact > nactx) then
nactx = nact
icount = 3
else
icount = icount - 1
if (icount == 0) go to 490
end if
!
! If ICON exceeds NACT, then we add the constraint with index IACT(ICON) to
! the active set. Apply Givens rotations so that the last N-NACT-1 columns
! of Z are orthogonal to the gradient of the new constraint, a scalar
! product being set to zero if its nonzero value could be due to computer
! rounding errors. The array DXNEW is used for working space.
!
if (icon <= nact) go to 260
kk = iact (icon)
do i = 1, n
dxnew (i) = a (i, kk)
end do
tot = 0.0_wp
k = n
100 if (k > nact) then
sp = 0.0_wp
spabs = 0.0_wp
do i = 1, n
temp = z (i, k) * dxnew (i)
sp = sp + temp
spabs = spabs + abs (temp)
end do
acca = spabs + 0.1_wp * abs (sp)
accb = spabs + 0.2_wp * abs (sp)
if (spabs >= acca .or. acca >= accb) sp = 0.0_wp
if (tot == 0.0_wp) then
tot = sp
else
kp = k + 1
temp = sqrt (sp*sp+tot*tot)
alpha = sp / temp
beta = tot / temp
tot = temp
do i = 1, n
temp = alpha * z (i, k) + beta * z (i, kp)
z (i, kp) = alpha * z (i, kp) - beta * z (i, k)
z (i, k) = temp
end do
end if
k = k - 1
go to 100
end if
!
! Add the new constraint if this can be done without a deletion from the
! active set.
!
if (tot /= 0.0_wp) then
nact = nact + 1
zdota (nact) = tot
vmultc (icon) = vmultc (nact)
vmultc (nact) = 0.0_wp
go to 210
end if
!
! The next instruction is reached if a deletion has to be made from the
! active set in order to make room for the new active constraint, because
! the new constraint gradient is a linear combination of the gradients of
! the old active constraints. Set the elements of VMULTD to the multipliers
! of the linear combination. Further, set IOUT to the index of the
! constraint to be deleted, but branch if no suitable index can be found.
!
ratio = - 1.0_wp
k = nact
130 zdotv = 0.0_wp
zdvabs = 0.0_wp
do i = 1, n
temp = z (i, k) * dxnew (i)
zdotv = zdotv + temp
zdvabs = zdvabs + abs (temp)
end do
acca = zdvabs + 0.1_wp * abs (zdotv)
accb = zdvabs + 0.2_wp * abs (zdotv)
if (zdvabs < acca .and. acca < accb) then
temp = zdotv / zdota (k)
if (temp > 0.0_wp .and. iact(k) <= m) then
tempa = vmultc (k) / temp
if (ratio < 0.0_wp .or. tempa < ratio) then
ratio = tempa
iout = k
end if
end if
if (k >= 2) then
kw = iact (k)
do i = 1, n
dxnew (i) = dxnew (i) - temp * a (i, kw)
end do
end if
vmultd (k) = temp
else
vmultd (k) = 0.0_wp
end if
k = k - 1
if (k > 0) go to 130
if (ratio < 0.0_wp) go to 490
!
! Revise the Lagrange multipliers and reorder the active constraints so
! that the one to be replaced is at the end of the list. Also calculate the
! new value of ZDOTA(NACT) and branch if it is not acceptable.
!
do k = 1, nact
vmultc (k) = max (0.0_wp, vmultc(k)-ratio*vmultd(k))
end do
if (icon < nact) then
isave = iact (icon)
vsave = vmultc (icon)
k = icon
170 kp = k + 1
kw = iact (kp)
sp = 0.0_wp
do i = 1, n
sp = sp + z (i, k) * a (i, kw)
end do
temp = sqrt (sp*sp+zdota(kp)**2)
alpha = zdota (kp) / temp
beta = sp / temp
zdota (kp) = alpha * zdota (k)
zdota (k) = temp
do i = 1, n
temp = alpha * z (i, kp) + beta * z (i, k)
z (i, kp) = alpha * z (i, k) - beta * z (i, kp)
z (i, k) = temp
end do
iact (k) = kw
vmultc (k) = vmultc (kp)
k = kp
if (k < nact) go to 170
iact (k) = isave
vmultc (k) = vsave
end if
temp = 0.0_wp
do i = 1, n
temp = temp + z (i, nact) * a (i, kk)
end do
if (temp == 0.0_wp) go to 490
zdota (nact) = temp
vmultc (icon) = 0.0_wp
vmultc (nact) = ratio
!
! Update IACT and ensure that the objective function continues to be
! treated as the last active constraint when MCON>M.
!
210 iact (icon) = iact (nact)
iact (nact) = kk
if (mcon > m .and. kk /= mcon) then
k = nact - 1
sp = 0.0_wp
do i = 1, n
sp = sp + z (i, k) * a (i, kk)
end do
temp = sqrt (sp*sp+zdota(nact)**2)
alpha = zdota (nact) / temp
beta = sp / temp
zdota (nact) = alpha * zdota (k)
zdota (k) = temp
do i = 1, n
temp = alpha * z (i, nact) + beta * z (i, k)
z (i, nact) = alpha * z (i, k) - beta * z (i, nact)
z (i, k) = temp
end do
iact (nact) = iact (k)
iact (k) = kk
temp = vmultc (k)
vmultc (k) = vmultc (nact)
vmultc (nact) = temp
end if
!
! If stage one is in progress, then set SDIRN to the direction of the next
! change to the current vector of variables.
!
if (mcon > m) go to 320
kk = iact (nact)
temp = 0.0_wp
do i = 1, n
temp = temp + sdirn (i) * a (i, kk)
end do
temp = temp - 1.0_wp
temp = temp / zdota (nact)
do i = 1, n
sdirn (i) = sdirn (i) - temp * z (i, nact)
end do
go to 340
!
! Delete the constraint that has the index IACT(ICON) from the active set.
!
260 if (icon < nact) then
isave = iact (icon)
vsave = vmultc (icon)
k = icon
270 kp = k + 1
kk = iact (kp)
sp = 0.0_wp
do i = 1, n
sp = sp + z (i, k) * a (i, kk)
end do
temp = sqrt (sp*sp+zdota(kp)**2)
alpha = zdota (kp) / temp
beta = sp / temp
zdota (kp) = alpha * zdota (k)
zdota (k) = temp
do i = 1, n
temp = alpha * z (i, kp) + beta * z (i, k)
z (i, kp) = alpha * z (i, k) - beta * z (i, kp)
z (i, k) = temp
end do
iact (k) = kk
vmultc (k) = vmultc (kp)
k = kp
if (k < nact) go to 270
iact (k) = isave
vmultc (k) = vsave
end if
nact = nact - 1
!
! If stage one is in progress, then set SDIRN to the direction of the next
! change to the current vector of variables.
!
if (mcon > m) go to 320
temp = 0.0_wp
do i = 1, n
temp = temp + sdirn (i) * z (i, nact+1)
end do
do i = 1, n
sdirn (i) = sdirn (i) - temp * z (i, nact+1)
end do
go to 340
!
! Pick the next search direction of stage two.
!
320 temp = 1.0_wp / zdota (nact)
do i = 1, n
sdirn (i) = temp * z (i, nact)
end do
!
! Calculate the step to the boundary of the trust region or take the step
! that reduces RESMAX to zero. The two statements below that include the
! factor 1.0E-6 prevent some harmless underflows that occurred in a test
! calculation. Further, we skip the step if it could be zero within a
! reasonable tolerance for computer rounding errors.
!
340 dd = rho * rho
sd = 0.0_wp
ss = 0.0_wp
do i = 1, n
if (abs(dx(i)) >= 1.0e-6_wp*rho) dd = dd - dx (i) ** 2
sd = sd + dx (i) * sdirn (i)
ss = ss + sdirn (i) ** 2
end do
if (dd <= 0.0_wp) go to 490
temp = sqrt (ss*dd)
if (abs(sd) >= 1.0e-6_wp*temp) temp = sqrt (ss*dd+sd*sd)
stpful = dd / (temp+sd)
step = stpful
if (mcon == m) then
acca = step + 0.1_wp * resmax
accb = step + 0.2_wp * resmax
if (step >= acca .or. acca >= accb) go to 480
step = min (step, resmax)
end if
!
! Set DXNEW to the new variables if STEP is the steplength, and reduce
! RESMAX to the corresponding maximum residual if stage one is being done.
! Because DXNEW will be changed during the calculation of some Lagrange
! multipliers, it will be restored to the following value later.
!
do i = 1, n
dxnew (i) = dx (i) + step * sdirn (i)
end do
if (mcon == m) then
resold = resmax
resmax = 0.0_wp
do k = 1, nact
kk = iact (k)
temp = b (kk)
do i = 1, n
temp = temp - a (i, kk) * dxnew (i)
end do
resmax = max (resmax, temp)
end do
end if
!
! Set VMULTD to the VMULTC vector that would occur if DX became DXNEW. A
! device is included to force VMULTD(K)=0.0 if deviations from this value
! can be attributed to computer rounding errors. First calculate the new
! Lagrange multipliers.
!
k = nact
390 zdotw = 0.0_wp
zdwabs = 0.0_wp
do i = 1, n
temp = z (i, k) * dxnew (i)
zdotw = zdotw + temp
zdwabs = zdwabs + abs (temp)
end do
acca = zdwabs + 0.1_wp * abs (zdotw)
accb = zdwabs + 0.2_wp * abs (zdotw)
if (zdwabs >= acca .or. acca >= accb) zdotw = 0.0_wp
vmultd (k) = zdotw / zdota (k)
if (k >= 2) then
kk = iact (k)
do i = 1, n
dxnew (i) = dxnew (i) - vmultd (k) * a (i, kk)
end do
k = k - 1
go to 390
end if
if (mcon > m) vmultd (nact) = max (0.0_wp, vmultd(nact))
!
! Complete VMULTC by finding the new constraint residuals.
!
do i = 1, n
dxnew (i) = dx (i) + step * sdirn (i)
end do
if (mcon > nact) then
kl = nact + 1
do k = kl, mcon
kk = iact (k)
sum = resmax - b (kk)
sumabs = resmax + abs (b(kk))
do i = 1, n
temp = a (i, kk) * dxnew (i)
sum = sum + temp
sumabs = sumabs + abs (temp)
end do
acca = sumabs + 0.1_wp * abs (sum)
accb = sumabs + 0.2_wp * abs (sum)
if (sumabs >= acca .or. acca >= accb) sum = 0.0_wp
vmultd (k) = sum
end do
end if
!
! Calculate the fraction of the step from DX to DXNEW that will be taken.
!
ratio = 1.0_wp
icon = 0
do k = 1, mcon
if (vmultd(k) < 0.0_wp) then
temp = vmultc (k) / (vmultc(k)-vmultd(k))
if (temp < ratio) then
ratio = temp
icon = k
end if
end if
end do
!
! Update DX, VMULTC and RESMAX.
!
temp = 1.0_wp - ratio
do i = 1, n
dx (i) = temp * dx (i) + ratio * dxnew (i)
end do
do k = 1, mcon
vmultc (k) = max (0.0_wp, temp*vmultc(k)+ratio*vmultd(k))
end do
if (mcon == m) resmax = resold + ratio * (resmax-resold)
!
! If the full step is not acceptable then begin another iteration.
! Otherwise switch to stage two or end the calculation.
!
if (icon > 0) go to 70
if (step == stpful) return
480 mcon = m + 1
icon = mcon
iact (mcon) = mcon
vmultc (mcon) = 0.0_wp
go to 60
!
! We employ any freedom that may be available to reduce the objective
! function before returning a DX whose length is less than RHO.
!
490 if (mcon == m) go to 480
ifull = 0
end subroutine trstlp