subroutine lincob (n, npt, m, amat, b, x, rhobeg, rhoend, iprint, maxfun, xbase, xpt, &
fval, xsav, xopt, gopt, hq, pq, bmat, zmat, ndim, step, sp, xnew, &
iact, rescon, qfac, rfac, pqw, w, calfun)
implicit real (wp) (a-h, o-z)
dimension amat (n,*), b (*), x (*), xbase (*), xpt (npt,*), fval (*), xsav (*), &
xopt (*), gopt (*), hq (*), pq (*), bmat (ndim,*), zmat (npt,*), &
step (*), sp (*), xnew (*), iact (*), rescon (*), qfac (n,*), rfac (*), &
pqw (*), w (*)
procedure (func) :: calfun
!
! The arguments N, NPT, M, X, RHOBEG, RHOEND, IPRINT and MAXFUN are
! identical to the corresponding arguments in SUBROUTINE LINCOA.
! AMAT is a matrix whose columns are the constraint gradients, scaled
! so that they have unit length.
! B contains on entry the right hand sides of the constraints, scaled
! as above, but later B is modified for variables relative to XBASE.
! XBASE holds a shift of origin that should reduce the contributions
! from rounding errors to values of the model and Lagrange functions.
! XPT contains the interpolation point coordinates relative to XBASE.
! FVAL holds the values of F at the interpolation points.
! XSAV holds the best feasible vector of variables so far, without any
! shift of origin.
! XOPT is set to XSAV-XBASE, which is the displacement from XBASE of
! the feasible vector of variables that provides the least calculated
! F so far, this vector being the current trust region centre.
! GOPT holds the gradient of the quadratic model at XSAV = XBASE+XOPT.
! HQ holds the explicit second derivatives of the quadratic model.
! PQ contains the parameters of the implicit second derivatives of the
! quadratic model.
! BMAT holds the last N columns of the big inverse matrix H.
! ZMAT holds the factorization of the leading NPT by NPT submatrix
! of H, this factorization being ZMAT times Diag(DZ) times ZMAT^T,
! where the elements of DZ are plus or minus one, as specified by IDZ.
! NDIM is the first dimension of BMAT and has the value NPT+N.
! STEP is employed for trial steps from XOPT. It is also used for working
! space when XBASE is shifted and in PRELIM.
! SP is reserved for the scalar products XOPT^T XPT(K,.), K=1,2,...,NPT,
! followed by STEP^T XPT(K,.), K=1,2,...,NPT.
! XNEW is the displacement from XBASE of the vector of variables for
! the current calculation of F, except that SUBROUTINE TRSTEP uses it
! for working space.
! IACT is an integer array for the indices of the active constraints.
! RESCON holds useful information about the constraint residuals. Every
! nonnegative RESCON(J) is the residual of the J-th constraint at the
! current trust region centre. Otherwise, if RESCON(J) is negative, the
! J-th constraint holds as a strict inequality at the trust region
! centre, its residual being at least |RESCON(J)|; further, the value
! of |RESCON(J)| is at least the current trust region radius DELTA.
! QFAC is the orthogonal part of the QR factorization of the matrix of
! active constraint gradients, these gradients being ordered in
! accordance with IACT. When NACT is less than N, columns are added
! to QFAC to complete an N by N orthogonal matrix, which is important
! for keeping calculated steps sufficiently close to the boundaries
! of the active constraints.
! RFAC is the upper triangular part of this QR factorization, beginning
! with the first diagonal element, followed by the two elements in the
! upper triangular part of the second column and so on.
! PQW is used for working space, mainly for storing second derivative
! coefficients of quadratic functions. Its length is NPT+N.
! The array W is also used for working space. The required number of
! elements, namely MAX[M+3*N,2*M+N,2*NPT], is set in LINCOA.
!
! Set some constants.
!
real(wp),parameter :: half = 0.5_wp
real(wp),parameter :: one = 1.0_wp
real(wp),parameter :: tenth = 0.1_wp
real(wp),parameter :: zero = 0.0_wp
np = n + 1
nh = (n*np) / 2
nptm = npt - np
!
! Set the elements of XBASE, XPT, FVAL, XSAV, XOPT, GOPT, HQ, PQ, BMAT,
! ZMAT and SP for the first iteration. An important feature is that,
! if the interpolation point XPT(K,.) is not feasible, where K is any
! integer from [1,NPT], then a change is made to XPT(K,.) if necessary
! so that the constraint violation is at least 0.2*RHOBEG. Also KOPT
! is set so that XPT(KOPT,.) is the initial trust region centre.
!
call prelim (n, npt, m, amat, b, x, rhobeg, iprint, xbase, xpt, fval, xsav, xopt, &
gopt, kopt, hq, pq, bmat, zmat, idz, ndim, sp, rescon, step, pqw, w, &
calfun)
!
! Begin the iterative procedure.
!
nf = npt
fopt = fval (kopt)
rho = rhobeg
delta = rho
ifeas = 0
nact = 0
itest = 3
10 knew = 0
nvala = 0
nvalb = 0
!
! Shift XBASE if XOPT may be too far from XBASE. First make the changes
! to BMAT that do not depend on ZMAT.
!
20 fsave = fopt
xoptsq = zero
do i = 1, n
xoptsq = xoptsq + xopt (i) ** 2
end do
if (xoptsq >= 1.0e4_wp*delta*delta) then
qoptsq = 0.25_wp * xoptsq
do k = 1, npt
sum = zero
do i = 1, n
sum = sum + xpt (k, i) * xopt (i)
end do
sum = sum - half * xoptsq
w (npt+k) = sum
sp (k) = zero
do i = 1, n
xpt (k, i) = xpt (k, i) - half * xopt (i)
step (i) = bmat (k, i)
w (i) = sum * xpt (k, i) + qoptsq * xopt (i)
ip = npt + i
do j = 1, i
bmat (ip, j) = bmat (ip, j) + step (i) * w (j) + w (i) * step (j)
end do
end do
end do
!
! Then the revisions of BMAT that depend on ZMAT are calculated.
!
do k = 1, nptm
sumz = zero
do i = 1, npt
sumz = sumz + zmat (i, k)
w (i) = w (npt+i) * zmat (i, k)
end do
do j = 1, n
sum = qoptsq * sumz * xopt (j)
do i = 1, npt
sum = sum + w (i) * xpt (i, j)
end do
step (j) = sum
if (k < idz) sum = - sum
do i = 1, npt
bmat (i, j) = bmat (i, j) + sum * zmat (i, k)
end do
end do
do i = 1, n
ip = i + npt
temp = step (i)
if (k < idz) temp = - temp
do j = 1, i
bmat (ip, j) = bmat (ip, j) + temp * step (j)
end do
end do
end do
!
! Update the right hand sides of the constraints.
!
if (m > 0) then
do j = 1, m
temp = zero
do i = 1, n
temp = temp + amat (i, j) * xopt (i)
end do
b (j) = b (j) - temp
end do
end if
!
! The following instructions complete the shift of XBASE, including the
! changes to the parameters of the quadratic model.
!
ih = 0
do j = 1, n
w (j) = zero
do k = 1, npt
w (j) = w (j) + pq (k) * xpt (k, j)
xpt (k, j) = xpt (k, j) - half * xopt (j)
end do
do i = 1, j
ih = ih + 1
hq (ih) = hq (ih) + w (i) * xopt (j) + xopt (i) * w (j)
bmat (npt+i, j) = bmat (npt+j, i)
end do
end do
do j = 1, n
xbase (j) = xbase (j) + xopt (j)
xopt (j) = zero
xpt (kopt, j) = zero
end do
end if
!
! In the case KNEW=0, generate the next trust region step by calling
! TRSTEP, where SNORM is the current trust region radius initially.
! The final value of SNORM is the length of the calculated step,
! except that SNORM is zero on return if the projected gradient is
! unsuitable for starting the conjugate gradient iterations.
!
delsav = delta
ksave = knew
if (knew == 0) then
snorm = delta
do i = 1, n
xnew (i) = gopt (i)
end do
call trstep (n, npt, m, amat, b, xpt, hq, pq, nact, iact, rescon, qfac, rfac, &
snorm, step, xnew, w, w(m+1), pqw, pqw(np), w(m+np))
!
! A trust region step is applied whenever its length, namely SNORM, is at
! least HALF*DELTA. It is also applied if its length is at least 0.1999
! times DELTA and if a line search of TRSTEP has caused a change to the
! active set. Otherwise there is a branch below to label 530 or 560.
!
temp = half * delta
if (xnew(1) >= half) temp = 0.1999_wp * delta
if (snorm <= temp) then
delta = half * delta
if (delta <= 1.4_wp*rho) delta = rho
nvala = nvala + 1
nvalb = nvalb + 1
temp = snorm / rho
if (delsav > rho) temp = one
if (temp >= half) nvala = zero
if (temp >= tenth) nvalb = zero
if (delsav > rho) go to 530
if (nvala < 5 .and. nvalb < 3) go to 530
if (snorm > zero) ksave = - 1
go to 560
end if
nvala = zero
nvalb = zero
!
! Alternatively, KNEW is positive. Then the model step is calculated
! within a trust region of radius DEL, after setting the gradient at
! XBASE and the second derivative parameters of the KNEW-th Lagrange
! function in W(1) to W(N) and in PQW(1) to PQW(NPT), respectively.
!
else
del = max (tenth*delta, rho)
do i = 1, n
w (i) = bmat (knew, i)
end do
do k = 1, npt
pqw (k) = zero
end do
do j = 1, nptm
temp = zmat (knew, j)
if (j < idz) temp = - temp
do k = 1, npt
pqw (k) = pqw (k) + temp * zmat (k, j)
end do
end do
call qmstep (n, npt, m, amat, b, xpt, xopt, nact, iact, rescon, qfac, kopt, &
knew, del, step, w, pqw, w(np), w(np+m), ifeas)
end if
!
! Set VQUAD to the change to the quadratic model when the move STEP is
! made from XOPT. If STEP is a trust region step, then VQUAD should be
! negative. If it is nonnegative due to rounding errors in this case,
! there is a branch to label 530 to try to improve the model.
!
vquad = zero
ih = 0
do j = 1, n
vquad = vquad + step (j) * gopt (j)
do i = 1, j
ih = ih + 1
temp = step (i) * step (j)
if (i == j) temp = half * temp
vquad = vquad + temp * hq (ih)
end do
end do
do k = 1, npt
temp = zero
do j = 1, n
temp = temp + xpt (k, j) * step (j)
sp (npt+k) = temp
end do
vquad = vquad + half * pq (k) * temp * temp
end do
if (ksave == 0 .and. vquad >= zero) go to 530
!
! Calculate the next value of the objective function. The difference
! between the actual new value of F and the value predicted by the
! model is recorded in DIFF.
!
220 nf = nf + 1
if (nf > maxfun) then
nf = nf - 1
if (iprint > 0) print 230
230 format (/ 4 x, 'Return from LINCOA because CALFUN has been',&
' called MAXFUN times.')
go to 600
end if
xdiff = zero
do i = 1, n
xnew (i) = xopt (i) + step (i)
x (i) = xbase (i) + xnew (i)
xdiff = xdiff + (x(i)-xsav(i)) ** 2
end do
xdiff = sqrt (xdiff)
if (ksave ==-1) xdiff = rho
if (xdiff <= tenth*rho .or. xdiff >= delta+delta) then
ifeas = 0
if (iprint > 0) print 250
250 format (/ 4 x, 'Return from LINCOA because rounding errors',&
' prevent reasonable changes to X.')
go to 600
end if
if (ksave <= 0) ifeas = 1
f = real (ifeas, wp)
call calfun (n, x, f)
if (iprint == 3) then
print 260, nf, f, (x(i), i=1, n)
260 format (/ 4 x, 'Function number', i6, ' F =', 1 pd18.10,&
' The corresponding X is:' / (2 x, 5d15.6))
end if
if (ksave ==-1) go to 600
diff = f - fopt - vquad
!
! If X is feasible, then set DFFALT to the difference between the new
! value of F and the value predicted by the alternative model.
!
if (ifeas == 1 .and. itest < 3) then
do k = 1, npt
pqw (k) = zero
w (k) = fval (k) - fval (kopt)
end do
do j = 1, nptm
sum = zero
do i = 1, npt
sum = sum + w (i) * zmat (i, j)
end do
if (j < idz) sum = - sum
do k = 1, npt
pqw (k) = pqw (k) + sum * zmat (k, j)
end do
end do
vqalt = zero
do k = 1, npt
sum = zero
do j = 1, n
sum = sum + bmat (k, j) * step (j)
end do
vqalt = vqalt + sum * w (k)
vqalt = vqalt + pqw (k) * sp (npt+k) * (half*sp(npt+k)+sp(k))
end do
dffalt = f - fopt - vqalt
end if
if (itest == 3) then
dffalt = diff
itest = 0
end if
!
! Pick the next value of DELTA after a trust region step.
!
if (ksave == 0) then
ratio = (f-fopt) / vquad
if (ratio <= tenth) then
delta = half * delta
else if (ratio <= 0.7_wp) then
delta = max (half*delta, snorm)
else
temp = sqrt (2.0_wp) * delta
delta = max (half*delta, snorm+snorm)
delta = min (delta, temp)
end if
if (delta <= 1.4_wp*rho) delta = rho
end if
!
! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
! can be moved. If STEP is a trust region step, then KNEW is zero at
! present, but a positive value is picked by subroutine UPDATE.
!
call update (n, npt, xpt, bmat, zmat, idz, ndim, sp, step, kopt, knew, pqw, w)
if (knew == 0) then
if (iprint > 0) print 320
320 format (/ 4 x, &
'Return from LINCOA because the denominator of the updating formula is zero.')
go to 600
end if
!
! If ITEST is increased to 3, then the next quadratic model is the
! one whose second derivative matrix is least subject to the new
! interpolation conditions. Otherwise the new model is constructed
! by the symmetric Broyden method in the usual way.
!
if (ifeas == 1) then
itest = itest + 1
if (abs(dffalt) >= tenth*abs(diff)) itest = 0
end if
!
! Update the second derivatives of the model by the symmetric Broyden
! method, using PQW for the second derivative parameters of the new
! KNEW-th Lagrange function. The contribution from the old parameter
! PQ(KNEW) is included in the second derivative matrix HQ. W is used
! later for the gradient of the new KNEW-th Lagrange function.
!
if (itest < 3) then
do k = 1, npt
pqw (k) = zero
end do
do j = 1, nptm
temp = zmat (knew, j)
if (temp /= zero) then
if (j < idz) temp = - temp
do k = 1, npt
pqw (k) = pqw (k) + temp * zmat (k, j)
end do
end if
end do
ih = 0
do i = 1, n
w (i) = bmat (knew, i)
temp = pq (knew) * xpt (knew, i)
do j = 1, i
ih = ih + 1
hq (ih) = hq (ih) + temp * xpt (knew, j)
end do
end do
pq (knew) = zero
do k = 1, npt
pq (k) = pq (k) + diff * pqw (k)
end do
end if
!
! Include the new interpolation point with the corresponding updates of
! SP. Also make the changes of the symmetric Broyden method to GOPT at
! the old XOPT if ITEST is less than 3.
!
fval (knew) = f
sp (knew) = sp (kopt) + sp (npt+kopt)
ssq = zero
do i = 1, n
xpt (knew, i) = xnew (i)
ssq = ssq + step (i) ** 2
end do
sp (npt+knew) = sp (npt+kopt) + ssq
if (itest < 3) then
do k = 1, npt
temp = pqw (k) * sp (k)
do i = 1, n
w (i) = w (i) + temp * xpt (k, i)
end do
end do
do i = 1, n
gopt (i) = gopt (i) + diff * w (i)
end do
end if
!
! Update FOPT, XSAV, XOPT, KOPT, RESCON and SP if the new F is the
! least calculated value so far with a feasible vector of variables.
!
if (f < fopt .and. ifeas == 1) then
fopt = f
do j = 1, n
xsav (j) = x (j)
xopt (j) = xnew (j)
end do
kopt = knew
snorm = sqrt (ssq)
do j = 1, m
if (rescon(j) >= delta+snorm) then
rescon (j) = snorm - rescon (j)
else
rescon (j) = rescon (j) + snorm
if (rescon(j)+delta > zero) then
temp = b (j)
do i = 1, n
temp = temp - xopt (i) * amat (i, j)
end do
temp = max (temp, zero)
if (temp >= delta) temp = - temp
rescon (j) = temp
end if
end if
end do
do k = 1, npt
sp (k) = sp (k) + sp (npt+k)
end do
!
! Also revise GOPT when symmetric Broyden updating is applied.
!
if (itest < 3) then
ih = 0
do j = 1, n
do i = 1, j
ih = ih + 1
if (i < j) gopt (j) = gopt (j) + hq (ih) * step (i)
gopt (i) = gopt (i) + hq (ih) * step (j)
end do
end do
do k = 1, npt
temp = pq (k) * sp (npt+k)
do i = 1, n
gopt (i) = gopt (i) + temp * xpt (k, i)
end do
end do
end if
end if
!
! Replace the current model by the least Frobenius norm interpolant if
! this interpolant gives substantial reductions in the predictions
! of values of F at feasible points.
!
if (itest == 3) then
do k = 1, npt
pq (k) = zero
w (k) = fval (k) - fval (kopt)
end do
do j = 1, nptm
sum = zero
do i = 1, npt
sum = sum + w (i) * zmat (i, j)
end do
if (j < idz) sum = - sum
do k = 1, npt
pq (k) = pq (k) + sum * zmat (k, j)
end do
end do
do j = 1, n
gopt (j) = zero
do i = 1, npt
gopt (j) = gopt (j) + w (i) * bmat (i, j)
end do
end do
do k = 1, npt
temp = pq (k) * sp (k)
do i = 1, n
gopt (i) = gopt (i) + temp * xpt (k, i)
end do
end do
do ih = 1, nh
hq (ih) = zero
end do
end if
!
! If a trust region step has provided a sufficient decrease in F, then
! branch for another trust region calculation. Every iteration that
! takes a model step is followed by an attempt to take a trust region
! step.
!
knew = 0
if (ksave > 0) go to 20
if (ratio >= tenth) go to 20
!
! Alternatively, find out if the interpolation points are close enough
! to the best point so far.
!
530 distsq = max (delta*delta, 4.0_wp*rho*rho)
do k = 1, npt
sum = zero
do j = 1, n
sum = sum + (xpt(k, j)-xopt(j)) ** 2
end do
if (sum > distsq) then
knew = k
distsq = sum
end if
end do
!
! If KNEW is positive, then branch back for the next iteration, which
! will generate a "model step". Otherwise, if the current iteration
! has reduced F, or if DELTA was above its lower bound when the last
! trust region step was calculated, then try a "trust region" step
! instead.
!
if (knew > 0) go to 20
knew = 0
if (fopt < fsave) go to 20
if (delsav > rho) go to 20
!
! The calculations with the current value of RHO are complete.
! Pick the next value of RHO.
!
560 if (rho > rhoend) then
delta = half * rho
if (rho > 250.0_wp*rhoend) then
rho = tenth * rho
else if (rho <= 16.0_wp*rhoend) then
rho = rhoend
else
rho = sqrt (rho*rhoend)
end if
delta = max (delta, rho)
if (iprint >= 2) then
if (iprint >= 3) print 570
570 format (5 x)
print 580, rho, nf
580 format (/ 4 x, 'New RHO =', 1 pd11.4, 5 x, 'Number of',&
' function values =', i6)
print 590, fopt, (xbase(i)+xopt(i), i=1, n)
590 format (4 x, 'Least value of F =', 1 pd23.15, 9 x,&
'The corresponding X is:'/(2 x, 5d15.6))
end if
go to 10
end if
!
! Return from the calculation, after branching to label 220 for another
! Newton-Raphson step if it has not been tried before.
!
if (ksave ==-1) go to 220
600 if (fopt <= f .or. ifeas == 0) then
do i = 1, n
x (i) = xsav (i)
end do
f = fopt
end if
if (iprint >= 1) then
print 620, nf
620 format (/ 4 x, 'At the return from LINCOA', 5 x,&
'Number of function values =', i6)
print 590, f, (x(i), i=1, n)
end if
w (1) = f
w (2) = real (nf, wp) + half
end subroutine lincob