subroutine cobylb (n, m, mpp, x, rhobeg, rhoend, iprint, maxfun, con, sim, simi, &
datmat, a, vsig, veta, sigbar, dx, w, iact, calcfc)
implicit real (wp) (a-h, o-z)
dimension x (*), con (*), sim (n,*), simi (n,*), datmat (mpp,*), a (n,*), vsig(*),&
veta (*), sigbar (*), dx (*), w (*), iact (*)
procedure (func) :: calcfc
!
! Set the initial values of some parameters. The last column of SIM holds
! the optimal vertex of the current simplex, and the preceding N columns
! hold the displacements from the optimal vertex to the other vertices.
! Further, SIMI holds the inverse of the matrix that is contained in the
! first N columns of SIM.
!
iptem = min (n, 5)
iptemp = iptem + 1
np = n + 1
mp = m + 1
alpha = 0.25_wp
beta = 2.1_wp
gamma = 0.5_wp
delta = 1.1_wp
rho = rhobeg
parmu = 0.0_wp
if (iprint >= 2) print 10, rho
10 format (/ 3 x, 'The initial value of RHO is', 1 pe13.6, 2 x,&
'and PARMU is set to zero.')
nfvals = 0
temp = 1.0_wp / rho
do i = 1, n
sim (i, np) = x (i)
do j = 1, n
sim (i, j) = 0.0_wp
simi (i, j) = 0.0_wp
end do
sim (i, i) = rho
simi (i, i) = temp
end do
jdrop = np
ibrnch = 0
!
! Make the next call of the user-supplied subroutine CALCFC. These
! instructions are also used for calling CALCFC during the iterations of
! the algorithm.
!
40 if (nfvals >= maxfun .and. nfvals > 0) then
if (iprint >= 1) print 50
50 format (/ 3 x, 'Return from subroutine COBYLA because the ',&
'MAXFUN limit has been reached.')
go to 600
end if
nfvals = nfvals + 1
call calcfc (n, m, x, f, con)
resmax = 0.0_wp
if (m > 0) then
do k = 1, m
resmax = max (resmax,-con(k))
end do
end if
if (nfvals == iprint-1 .or. iprint == 3) then
print 70, nfvals, f, resmax, (x(i), i=1, iptem)
70 format (/ 3 x, 'NFVALS =', i5, 3 x, 'F =', 1 pe13.6, 4 x, 'MAXCV =', 1 pe13.6 &
& / 3 x, 'X =', 1 pe13.6, 1 p4e15.6)
if (iptem < n) print 80, (x(i), i=iptemp, n)
80 format (1 pe19.6, 1 p4e15.6)
end if
con (mp) = f
con (mpp) = resmax
if (ibrnch == 1) go to 440
!
! Set the recently calculated function values in a column of DATMAT. This
! array has a column for each vertex of the current simplex, the entries of
! each column being the values of the constraint functions (if any)
! followed by the objective function and the greatest constraint violation
! at the vertex.
!
do k = 1, mpp
datmat (k, jdrop) = con (k)
end do
if (nfvals > np) go to 130
!
! Exchange the new vertex of the initial simplex with the optimal vertex if
! necessary. Then, if the initial simplex is not complete, pick its next
! vertex and calculate the function values there.
!
if (jdrop <= n) then
if (datmat(mp, np) <= f) then
x (jdrop) = sim (jdrop, np)
else
sim (jdrop, np) = x (jdrop)
do k = 1, mpp
datmat (k, jdrop) = datmat (k, np)
datmat (k, np) = con (k)
end do
do k = 1, jdrop
sim (jdrop, k) = - rho
temp = 0.0_wp
do i = k, jdrop
temp = temp - simi (i, k)
end do
simi (jdrop, k) = temp
end do
end if
end if
if (nfvals <= n) then
jdrop = nfvals
x (jdrop) = x (jdrop) + rho
go to 40
end if
130 ibrnch = 1
!
! Identify the optimal vertex of the current simplex.
!
140 phimin = datmat (mp, np) + parmu * datmat (mpp, np)
nbest = np
do j = 1, n
temp = datmat (mp, j) + parmu * datmat (mpp, j)
if (temp < phimin) then
nbest = j
phimin = temp
else if (temp == phimin .and. parmu == 0.0_wp) then
if (datmat(mpp, j) < datmat(mpp, nbest)) nbest = j
end if
end do
!
! Switch the best vertex into pole position if it is not there already,
! and also update SIM, SIMI and DATMAT.
!
if (nbest <= n) then
do i = 1, mpp
temp = datmat (i, np)
datmat (i, np) = datmat (i, nbest)
datmat (i, nbest) = temp
end do
do i = 1, n
temp = sim (i, nbest)
sim (i, nbest) = 0.0_wp
sim (i, np) = sim (i, np) + temp
tempa = 0.0_wp
do k = 1, n
sim (i, k) = sim (i, k) - temp
tempa = tempa - simi (k, i)
end do
simi (nbest, i) = tempa
end do
end if
!
! Make an error return if SIGI is a poor approximation to the inverse of
! the leading N by N submatrix of SIG.
!
error = 0.0_wp
do i = 1, n
do j = 1, n
temp = 0.0_wp
if (i == j) temp = temp - 1.0_wp
do k = 1, n
temp = temp + simi (i, k) * sim (k, j)
end do
error = max (error, abs(temp))
end do
end do
if (error > 0.1_wp) then
if (iprint >= 1) print 210
210 format (/ 3 x, 'Return from subroutine COBYLA because ',&
'rounding errors are becoming damaging.')
go to 600
end if
!
! Calculate the coefficients of the linear approximations to the objective
! and constraint functions, placing minus the objective function gradient
! after the constraint gradients in the array A. The vector W is used for
! working space.
!
do k = 1, mp
con (k) = - datmat (k, np)
do j = 1, n
w (j) = datmat (k, j) + con (k)
end do
do i = 1, n
temp = 0.0_wp
do j = 1, n
temp = temp + w (j) * simi (j, i)
end do
if (k == mp) temp = - temp
a (i, k) = temp
end do
end do
!
! Calculate the values of sigma and eta, and set IFLAG=0 if the current
! simplex is not acceptable.
!
iflag = 1
parsig = alpha * rho
pareta = beta * rho
do j = 1, n
wsig = 0.0_wp
weta = 0.0_wp
do i = 1, n
wsig = wsig + simi (j, i) ** 2
weta = weta + sim (i, j) ** 2
end do
vsig (j) = 1.0_wp / sqrt (wsig)
veta (j) = sqrt (weta)
if (vsig(j) < parsig .or. veta(j) > pareta) iflag = 0
end do
!
! If a new vertex is needed to improve acceptability, then decide which
! vertex to drop from the simplex.
!
if (ibrnch == 1 .or. iflag == 1) go to 370
jdrop = 0
temp = pareta
do j = 1, n
if (veta(j) > temp) then
jdrop = j
temp = veta (j)
end if
end do
if (jdrop == 0) then
do j = 1, n
if (vsig(j) < temp) then
jdrop = j
temp = vsig (j)
end if
end do
end if
!
! Calculate the step to the new vertex and its sign.
!
temp = gamma * rho * vsig (jdrop)
do i = 1, n
dx (i) = temp * simi (jdrop, i)
end do
cvmaxp = 0.0_wp
cvmaxm = 0.0_wp
do k = 1, mp
sum = 0.0_wp
do i = 1, n
sum = sum + a (i, k) * dx (i)
end do
if (k < mp) then
temp = datmat (k, np)
cvmaxp = max (cvmaxp,-sum-temp)
cvmaxm = max (cvmaxm, sum-temp)
end if
end do
dxsign = 1.0_wp
if (parmu*(cvmaxp-cvmaxm) > sum+sum) dxsign = - 1.0_wp
!
! Update the elements of SIM and SIMI, and set the next X.
!
temp = 0.0_wp
do i = 1, n
dx (i) = dxsign * dx (i)
sim (i, jdrop) = dx (i)
temp = temp + simi (jdrop, i) * dx (i)
end do
do i = 1, n
simi (jdrop, i) = simi (jdrop, i) / temp
end do
do j = 1, n
if (j /= jdrop) then
temp = 0.0_wp
do i = 1, n
temp = temp + simi (j, i) * dx (i)
end do
do i = 1, n
simi (j, i) = simi (j, i) - temp * simi (jdrop, i)
end do
end if
x (j) = sim (j, np) + dx (j)
end do
go to 40
!
! Calculate DX=x(*)-x(0). Branch if the length of DX is less than 0.5*RHO.
!
370 iz = 1
izdota = iz + n * n
ivmc = izdota + n
isdirn = ivmc + mp
idxnew = isdirn + n
ivmd = idxnew + n
call trstlp (n, m, a, con, rho, dx, ifull, iact, w(iz), w(izdota), w(ivmc), &
& w(isdirn), w(idxnew), w(ivmd))
if (ifull == 0) then
temp = 0.0_wp
do i = 1, n
temp = temp + dx (i) ** 2
end do
if (temp < 0.25_wp*rho*rho) then
ibrnch = 1
go to 550
end if
end if
!
! Predict the change to F and the new maximum constraint violation if the
! variables are altered from x(0) to x(0)+DX.
!
resnew = 0.0_wp
con (mp) = 0.0_wp
do k = 1, mp
sum = con (k)
do i = 1, n
sum = sum - a (i, k) * dx (i)
end do
if (k < mp) resnew = max (resnew, sum)
end do
!
! Increase PARMU if necessary and branch back if this change alters the
! optimal vertex. Otherwise PREREM and PREREC will be set to the predicted
! reductions in the merit function and the maximum constraint violation
! respectively.
!
barmu = 0.0_wp
prerec = datmat (mpp, np) - resnew
if (prerec > 0.0_wp) barmu = sum / prerec
if (parmu < 1.5_wp*barmu) then
parmu = 2.0_wp * barmu
if (iprint >= 2) print 410, parmu
410 format (/ 3 x, 'Increase in PARMU to', 1 pe13.6)
phi = datmat (mp, np) + parmu * datmat (mpp, np)
do j = 1, n
temp = datmat (mp, j) + parmu * datmat (mpp, j)
if (temp < phi) go to 140
if (temp == phi .and. parmu == 0.0_wp) then
if (datmat(mpp, j) < datmat(mpp, np)) go to 140
end if
end do
end if
prerem = parmu * prerec - sum
!
! Calculate the constraint and objective functions at x(*). Then find the
! actual reduction in the merit function.
!
do i = 1, n
x (i) = sim (i, np) + dx (i)
end do
ibrnch = 1
go to 40
440 vmold = datmat (mp, np) + parmu * datmat (mpp, np)
vmnew = f + parmu * resmax
trured = vmold - vmnew
if (parmu == 0.0_wp .and. f == datmat(mp, np)) then
prerem = prerec
trured = datmat (mpp, np) - resmax
end if
!
! Begin the operations that decide whether x(*) should replace one of the
! vertices of the current simplex, the change being mandatory if TRURED is
! positive. Firstly, JDROP is set to the index of the vertex that is to be
! replaced.
!
ratio = 0.0_wp
if (trured <= 0.0_wp) ratio = 1.0_wp
jdrop = 0
do j = 1, n
temp = 0.0_wp
do i = 1, n
temp = temp + simi (j, i) * dx (i)
end do
temp = abs (temp)
if (temp > ratio) then
jdrop = j
ratio = temp
end if
sigbar (j) = temp * vsig (j)
end do
!
! Calculate the value of ell.
!
edgmax = delta * rho
l = 0
do j = 1, n
if (sigbar(j) >= parsig .or. sigbar(j) >= vsig(j)) then
temp = veta (j)
if (trured > 0.0_wp) then
temp = 0.0_wp
do i = 1, n
temp = temp + (dx(i)-sim(i, j)) ** 2
end do
temp = sqrt (temp)
end if
if (temp > edgmax) then
l = j
edgmax = temp
end if
end if
end do
if (l > 0) jdrop = l
if (jdrop == 0) go to 550
!
! Revise the simplex by updating the elements of SIM, SIMI and DATMAT.
!
temp = 0.0_wp
do i = 1, n
sim (i, jdrop) = dx (i)
temp = temp + simi (jdrop, i) * dx (i)
end do
do i = 1, n
simi (jdrop, i) = simi (jdrop, i) / temp
end do
do j = 1, n
if (j /= jdrop) then
temp = 0.0_wp
do i = 1, n
temp = temp + simi (j, i) * dx (i)
end do
do i = 1, n
simi (j, i) = simi (j, i) - temp * simi (jdrop, i)
end do
end if
end do
do k = 1, mpp
datmat (k, jdrop) = con (k)
end do
!
! Branch back for further iterations with the current RHO.
!
if (trured > 0.0_wp .and. trured >= 0.1_wp*prerem) go to 140
550 if (iflag == 0) then
ibrnch = 0
go to 140
end if
!
! Otherwise reduce RHO if it is not at its least value and reset PARMU.
!
if (rho > rhoend) then
rho = 0.5_wp * rho
if (rho <= 1.5_wp*rhoend) rho = rhoend
if (parmu > 0.0_wp) then
denom = 0.0_wp
do k = 1, mp
cmin = datmat (k, np)
cmax = cmin
do i = 1, n
cmin = min (cmin, datmat(k, i))
cmax = max (cmax, datmat(k, i))
end do
if (k <= m .and. cmin < 0.5_wp*cmax) then
temp = max (cmax, 0.0_wp) - cmin
if (denom <= 0.0_wp) then
denom = temp
else
denom = min (denom, temp)
end if
end if
end do
if (denom == 0.0_wp) then
parmu = 0.0_wp
else if (cmax-cmin < parmu*denom) then
parmu = (cmax-cmin) / denom
end if
end if
if (iprint >= 2) print 580, rho, parmu
580 format (/ 3 x, 'Reduction in RHO to', 1 pe13.6, ' and PARMU =', 1 pe13.6)
if (iprint == 2) then
print 70, nfvals, datmat (mp, np), datmat (mpp, np), (sim(i, np), i=1, &
& iptem)
if (iptem < n) print 80, (x(i), i=iptemp, n)
end if
go to 140
end if
!
! Return the best calculated values of the variables.
!
if (iprint >= 1) print 590
590 format (/ 3 x, 'Normal return from subroutine COBYLA')
if (ifull == 1) go to 620
600 do i = 1, n
x (i) = sim (i, np)
end do
f = datmat (mp, np)
resmax = datmat (mpp, np)
620 if (iprint >= 1) then
print 70, nfvals, f, resmax, (x(i), i=1, iptem)
if (iptem < n) print 80, (x(i), i=iptemp, n)
end if
maxfun = nfvals
end subroutine cobylb