| 1 |
************************************************************************ |
| 2 |
* |
| 3 |
* File mi80ncon fortran. |
| 4 |
* |
| 5 |
* m8ajac m8augl m8aug1 m8chkj m8prtj m8sclj |
| 6 |
* m8setj m8viol |
| 7 |
* |
| 8 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 9 |
|
| 10 |
subroutine m8ajac( gotjac, nncon, nnjac, njac, |
| 11 |
$ ne, nka, a, ha, ka, |
| 12 |
$ fcon, fcon2, gcon, gcon2, xn, y, z, nwcore ) |
| 13 |
|
| 14 |
implicit double precision (a-h,o-z) |
| 15 |
logical gotjac |
| 16 |
integer*4 ha(ne) |
| 17 |
integer ka(nka) |
| 18 |
double precision a(ne) |
| 19 |
double precision fcon(nncon), fcon2(nncon), |
| 20 |
$ gcon(njac ), gcon2(njac ), |
| 21 |
$ xn (nnjac), y(nncon), z(nwcore) |
| 22 |
|
| 23 |
* ------------------------------------------------------------------ |
| 24 |
* m8ajac stores the Jacobian into A. |
| 25 |
* If gotjac is true, the Jacobian is already in gcon. |
| 26 |
* ------------------------------------------------------------------ |
| 27 |
|
| 28 |
common /m5log1/ idebug,ierr,lprint |
| 29 |
common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 30 |
|
| 31 |
if ( .not. gotjac ) then |
| 32 |
if (lderiv .lt. 2) |
| 33 |
$ call m6dmmy( njac, gcon ) |
| 34 |
|
| 35 |
call m6fcon( 2, nncon, nnjac, njac, fcon, gcon, |
| 36 |
$ ne, nka, ha, ka, |
| 37 |
$ xn, z, nwcore ) |
| 38 |
|
| 39 |
if (lderiv .lt. 2) |
| 40 |
$ call m6dcon( nncon, nnjac, njac, |
| 41 |
$ ne, nka, ha, ka, |
| 42 |
$ fcon, fcon2, gcon, gcon2, xn, y, z, nwcore ) |
| 43 |
if (ierr .ne. 0) return |
| 44 |
end if |
| 45 |
|
| 46 |
l = 0 |
| 47 |
do 400 j = 1, nnjac |
| 48 |
k1 = ka(j) |
| 49 |
k2 = ka(j+1) - 1 |
| 50 |
do 300 k = k1, k2 |
| 51 |
ir = ha(k) |
| 52 |
if (ir .gt. nncon) go to 400 |
| 53 |
l = l + 1 |
| 54 |
a(k) = gcon(l) |
| 55 |
300 continue |
| 56 |
400 continue |
| 57 |
|
| 58 |
* end of m8ajac |
| 59 |
end |
| 60 |
|
| 61 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 62 |
|
| 63 |
subroutine m8augl( mode, m, n, nb, ns, inform, |
| 64 |
$ ne, nka, a, ha, ka, |
| 65 |
$ hs, bl, bu, xn, z, nwcore ) |
| 66 |
|
| 67 |
implicit double precision (a-h,o-z) |
| 68 |
integer*4 ha(ne), hs(nb) |
| 69 |
integer ka(nka) |
| 70 |
double precision a(ne), bl(nb), bu(nb), xn(nb), z(nwcore) |
| 71 |
|
| 72 |
* ------------------------------------------------------------------ |
| 73 |
* m8augl is a front-end for m8aug1. It is called by misolv and |
| 74 |
* m5solv at various stages of the augmented Lagrangian algorithm. |
| 75 |
* 11 Oct 1991: a, ha, ka etc added as parameters. |
| 76 |
* 05 Mar 1992: mode 3 now makes nonlinear rows free for first major. |
| 77 |
* 11 Sep 1992: hs added as parameter. |
| 78 |
* ------------------------------------------------------------------ |
| 79 |
|
| 80 |
common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu , |
| 81 |
$ lhrtyp,lhs ,lkb |
| 82 |
common /m5loc / lpi ,lpi2 ,lw ,lw2 , |
| 83 |
$ lx ,lx2 ,ly ,ly2 , |
| 84 |
$ lgsub ,lgsub2,lgrd ,lgrd2 , |
| 85 |
$ lr ,lrg ,lrg2 ,lxn |
| 86 |
common /m8len / njac ,nncon ,nncon0,nnjac |
| 87 |
common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold , |
| 88 |
$ lblslk,lbuslk,lxlam ,lrhs , |
| 89 |
$ lgcon ,lgcon2,lxdif ,lxold |
| 90 |
|
| 91 |
ms = m + ns |
| 92 |
call m8aug1( mode, ms, nncon, nnjac, njac, n, nb, inform, |
| 93 |
$ ne, nka, a, ha, ka, |
| 94 |
$ hs, z(lkb), bl, bu, z(lbbl), z(lbbu), |
| 95 |
$ z(lblslk), z(lbuslk), |
| 96 |
$ z(lgcon ), z(lgcon2), z(lxlam), xn ) |
| 97 |
|
| 98 |
* end of m8augl |
| 99 |
end |
| 100 |
|
| 101 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 102 |
|
| 103 |
subroutine m8aug1( mode, ms, nncon, nnjac, njac, n, nb, inform, |
| 104 |
$ ne, nka, a, ha, ka, |
| 105 |
$ hs, kb, bl, bu, bbl, bbu, blslk, buslk, |
| 106 |
$ gcon, gcon2, xlam, xn ) |
| 107 |
|
| 108 |
implicit double precision (a-h,o-z) |
| 109 |
integer*4 ha(ne), hs(nb) |
| 110 |
integer ka(nka), kb(ms) |
| 111 |
double precision a(ne), bl(nb), bu(nb), bbl(ms), bbu(ms), |
| 112 |
$ blslk(nncon), buslk(nncon), |
| 113 |
$ gcon (njac ), gcon2(njac ), |
| 114 |
$ xlam (nncon), xn(nb) |
| 115 |
|
| 116 |
* ------------------------------------------------------------------ |
| 117 |
* m8aug1 does the work for m8augl. |
| 118 |
* 11 Sep 1992: hs added as parameter for mode 5, to allow |
| 119 |
* nonbasic slacks to be moved onto the relaxed bounds. |
| 120 |
* ------------------------------------------------------------------ |
| 121 |
|
| 122 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 123 |
common /m1file/ iread,iprint,isumm |
| 124 |
common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj |
| 125 |
common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm |
| 126 |
common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor |
| 127 |
common /m8al2 / radius,rhsmod,modpen,modrhs |
| 128 |
|
| 129 |
intrinsic max, min |
| 130 |
parameter ( zero = 0.0d+0 ) |
| 131 |
|
| 132 |
inform = 0 |
| 133 |
bplus = 0.9*plinfy |
| 134 |
bminus = - bplus |
| 135 |
|
| 136 |
if (mode .eq. 1) then |
| 137 |
* --------------------------------------------------------------- |
| 138 |
* Save the Jacobian from the initial a(*) in gcon, gcon2. |
| 139 |
* --------------------------------------------------------------- |
| 140 |
l = 0 |
| 141 |
do 140 j = 1, nnjac |
| 142 |
k1 = ka(j) |
| 143 |
k2 = ka(j+1) - 1 |
| 144 |
do 120 k = k1, k2 |
| 145 |
ir = ha(k) |
| 146 |
if (ir .gt. nncon) go to 140 |
| 147 |
l = l + 1 |
| 148 |
gcon (l) = a(k) |
| 149 |
gcon2(l) = a(k) |
| 150 |
120 continue |
| 151 |
140 continue |
| 152 |
|
| 153 |
else if (mode .eq. 2) then |
| 154 |
* --------------------------------------------------------------- |
| 155 |
* Make sure the initial xn values are feasible. |
| 156 |
* Do the slacks too, in case the RHS has been perturbed. |
| 157 |
* --------------------------------------------------------------- |
| 158 |
do 220 j = 1, nb |
| 159 |
xn(j) = max( xn(j), bl(j) ) |
| 160 |
xn(j) = min( xn(j), bu(j) ) |
| 161 |
220 continue |
| 162 |
|
| 163 |
else if (mode .eq. 3) then |
| 164 |
* --------------------------------------------------------------- |
| 165 |
* Start of the first major iteration. |
| 166 |
* 1. Initialize modpen, rhsmod. |
| 167 |
* 2. Save the bounds on the nonlinear rows, |
| 168 |
* 3. Make the nonlinear rows free (relax the bounds to infinity). |
| 169 |
* --------------------------------------------------------------- |
| 170 |
modpen = 0 |
| 171 |
rhsmod = zero |
| 172 |
call dcopy ( nncon, bl(n+1), 1, blslk, 1 ) |
| 173 |
call dcopy ( nncon, bu(n+1), 1, buslk, 1 ) |
| 174 |
|
| 175 |
do 320 i = 1, nncon |
| 176 |
bl(n+i) = - plinfy |
| 177 |
bu(n+i) = plinfy |
| 178 |
320 continue |
| 179 |
|
| 180 |
else if (mode .eq. 4) then |
| 181 |
* --------------------------------------------------------------- |
| 182 |
* Unbounded subproblem. |
| 183 |
* Increase the penalty parameter. |
| 184 |
* --------------------------------------------------------------- |
| 185 |
modpen = modpen + 1 |
| 186 |
if (modpen .gt. 5 .or. nlag .eq. 0) then |
| 187 |
inform = 1 |
| 188 |
else |
| 189 |
t = nncon |
| 190 |
if (penpar .le. zero) penpar = t / 100.0 |
| 191 |
penpar = 10.0 * penpar |
| 192 |
if (iprint .gt. 0) write(iprint, 1400) penpar |
| 193 |
if (isumm .gt. 0) write(isumm , 1400) penpar |
| 194 |
end if |
| 195 |
|
| 196 |
else if (mode .eq. 5) then |
| 197 |
* --------------------------------------------------------------- |
| 198 |
* Infeasible subproblem. |
| 199 |
* Relax the bounds on the linearized constraints. |
| 200 |
* Also, move nonbasic slacks onto the new bounds |
| 201 |
* so there won't be large numbers of them floating in between. |
| 202 |
* --------------------------------------------------------------- |
| 203 |
modrhs = modrhs + 1 |
| 204 |
if (modrhs .gt. 5) then |
| 205 |
inform = 1 |
| 206 |
else |
| 207 |
if (modrhs .eq. 1) rhsmod = tolx |
| 208 |
t = sinf / nncon |
| 209 |
rhsmod = 10.0 * max( rhsmod, t ) |
| 210 |
if (iprint .gt. 0) write(iprint, 1500) rhsmod |
| 211 |
if (isumm .gt. 0) write(isumm , 1500) rhsmod |
| 212 |
|
| 213 |
do 520 i = 1, nncon |
| 214 |
j = n + i |
| 215 |
if (bl(j) .gt. bminus) then |
| 216 |
bl(j) = blslk(i) - rhsmod |
| 217 |
if (hs(j) .eq. 0) xn(j) = bl(j) |
| 218 |
end if |
| 219 |
if (bu(j) .lt. bplus ) then |
| 220 |
bu(j) = buslk(i) + rhsmod |
| 221 |
if (hs(j) .eq. 1) xn(j) = bu(j) |
| 222 |
end if |
| 223 |
520 continue |
| 224 |
|
| 225 |
do 540 k = 1, ms |
| 226 |
j = kb(k) |
| 227 |
bbl(k) = bl(j) |
| 228 |
bbu(k) = bu(j) |
| 229 |
540 continue |
| 230 |
end if |
| 231 |
|
| 232 |
else if (mode .eq. 6) then |
| 233 |
* --------------------------------------------------------------- |
| 234 |
* Restore the bounds on the nonlinear constraints. |
| 235 |
* --------------------------------------------------------------- |
| 236 |
do 620 i = 1, nncon |
| 237 |
jslack = n + i |
| 238 |
bl(jslack) = blslk(i) |
| 239 |
bu(jslack) = buslk(i) |
| 240 |
620 continue |
| 241 |
end if |
| 242 |
return |
| 243 |
|
| 244 |
1400 format(' Penalty parameter increased to', 1p, e12.2) |
| 245 |
1500 format(' XXX Infeasible subproblem. ', |
| 246 |
$ ' RHS perturbation increased to', 1p, e12.2) |
| 247 |
|
| 248 |
* end of m8aug1 |
| 249 |
end |
| 250 |
|
| 251 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 252 |
|
| 253 |
subroutine m8chkj( m, n, njac, nx, |
| 254 |
$ ne, nka, ha, ka, |
| 255 |
$ bl, bu, f, f2, g, g2, |
| 256 |
$ x, y, y2, z, nwcore ) |
| 257 |
|
| 258 |
implicit double precision (a-h,o-z) |
| 259 |
integer*4 ha(ne) |
| 260 |
integer ka(nka) |
| 261 |
double precision bl(n), bu(n), f(m), f2(m), g(njac), g2(njac), |
| 262 |
$ x(n), y(n), y2(nx), z(nwcore) |
| 263 |
|
| 264 |
* ------------------------------------------------------------------ |
| 265 |
* m8chkj verifies the Jacobian gcon using |
| 266 |
* finite differences on the constraint functions f. |
| 267 |
* The various verify levels are described in m7chkg. |
| 268 |
* ------------------------------------------------------------------ |
| 269 |
|
| 270 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 271 |
common /m1file/ iread,iprint,isumm |
| 272 |
common /m3scal/ sclobj,scltol,lscale |
| 273 |
common /m5log1/ idebug,ierr,lprint |
| 274 |
common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 275 |
common /m8veri/ jverif(4),lverif(2) |
| 276 |
|
| 277 |
intrinsic abs, max, min |
| 278 |
parameter ( zero = 0.0d+0, one = 1.0d+0, ok = 0.1d+0 ) |
| 279 |
|
| 280 |
logical cheap, first |
| 281 |
character*4 key |
| 282 |
character*4 lbad , lgood |
| 283 |
data lbad /'bad?'/, lgood /'ok '/ |
| 284 |
|
| 285 |
lvl = lverif(1) |
| 286 |
if (lvl .lt. 0) return |
| 287 |
j1 = max( jverif(3), 1 ) |
| 288 |
j2 = min( jverif(4), n ) |
| 289 |
cheap = lvl .le. 1 .or. j1 .gt. j2 |
| 290 |
lssave = lscale |
| 291 |
lscale = 0 |
| 292 |
|
| 293 |
* Evaluate f and g at the base point x. |
| 294 |
|
| 295 |
if (lderiv .le. 1) call m6dmmy( njac, g ) |
| 296 |
call m6fcon( 2, m, n, njac, f, g, |
| 297 |
$ ne, nka, ha, ka, |
| 298 |
$ x, z, nwcore ) |
| 299 |
if (ierr .ne. 0) go to 900 |
| 300 |
if (knowng(2) .eq. 0) go to 900 |
| 301 |
|
| 302 |
if (iprint .gt. 0) then |
| 303 |
if (cheap) then |
| 304 |
write(iprint, 1100) |
| 305 |
else |
| 306 |
write(iprint, 1000) |
| 307 |
end if |
| 308 |
end if |
| 309 |
|
| 310 |
* -------------------------- |
| 311 |
* Cheap test. |
| 312 |
* -------------------------- |
| 313 |
|
| 314 |
* Generate a direction in which to perturb x. |
| 315 |
|
| 316 |
yj = one/n |
| 317 |
do 10 j = 1, n |
| 318 |
y(j) = yj |
| 319 |
y2(j) = yj |
| 320 |
yj = - yj * 0.99999 |
| 321 |
10 continue |
| 322 |
|
| 323 |
* Define a difference interval. |
| 324 |
* If needed, alter y to ensure that it will be a feasible direction. |
| 325 |
* If this gives zero, go back to original y and forget feasibility. |
| 326 |
|
| 327 |
dx = difint(1) * (one + dnorm1( n, x, 1 )) |
| 328 |
call m7chkd( n, bl, bu, x, dx, y, nfeas ) |
| 329 |
if (nfeas .eq. 0) call dcopy ( n, y2, 1, y, 1 ) |
| 330 |
|
| 331 |
* ------------------------------------------------------------------ |
| 332 |
* We must not perturb x(j) if the j-th column of the Jacobian |
| 333 |
* contains any unknown elements. |
| 334 |
* ------------------------------------------------------------------ |
| 335 |
nder = 0 |
| 336 |
l = 0 |
| 337 |
do 30 j = 1, n |
| 338 |
k1 = ka(j) |
| 339 |
k2 = ka(j+1) - 1 |
| 340 |
|
| 341 |
do 20 k = k1, k2 |
| 342 |
ir = ha(k) |
| 343 |
if (ir .gt. m) go to 25 |
| 344 |
l = l + 1 |
| 345 |
if (g(l) .eq. gdummy) y(j) = zero |
| 346 |
20 continue |
| 347 |
|
| 348 |
25 if (y(j) .ne. zero) nder = nder + 1 |
| 349 |
30 continue |
| 350 |
|
| 351 |
if (nder .eq. 0) go to 900 |
| 352 |
|
| 353 |
* Set f2 = constraint values at a short step along y. |
| 354 |
|
| 355 |
do 40 j = 1, n |
| 356 |
y2(j) = x(j) + dx*y(j) |
| 357 |
40 continue |
| 358 |
|
| 359 |
call m6fcon( 0, m, n, njac, f2, g2, |
| 360 |
$ ne, nka, ha, ka, |
| 361 |
$ y2, z, nwcore ) |
| 362 |
if (ierr .ne. 0) go to 900 |
| 363 |
|
| 364 |
* Set y2 = (f2 - f)/dx - Jacobian*y. This should be small. |
| 365 |
* At the same time, find the first Jacobian element in column j1. |
| 366 |
|
| 367 |
do 60 i = 1, m |
| 368 |
y2(i) = (f2(i) - f(i)) / dx |
| 369 |
60 continue |
| 370 |
|
| 371 |
l = 0 |
| 372 |
lsave = 0 |
| 373 |
do 100 j = 1, n |
| 374 |
yj = y(j) |
| 375 |
k1 = ka(j) |
| 376 |
k2 = ka(j+1) - 1 |
| 377 |
do 80 k = k1, k2 |
| 378 |
ir = ha(k) |
| 379 |
if (ir .gt. m) go to 100 |
| 380 |
l = l + 1 |
| 381 |
if (j .lt. j1) lsave = l |
| 382 |
y2(ir) = y2(ir) - g(l)*yj |
| 383 |
80 continue |
| 384 |
100 continue |
| 385 |
|
| 386 |
imax = idamax( m,y2,1 ) |
| 387 |
gmax = (f2(imax) - f(imax)) / dx |
| 388 |
emax = abs( y2(imax) ) / (one + abs( gmax )) |
| 389 |
if (emax .le. ok) then |
| 390 |
if (iprint .gt. 0) write(iprint, 1400) |
| 391 |
else |
| 392 |
if (iprint .gt. 0) write(iprint, 1500) |
| 393 |
if (isumm .gt. 0) write(isumm , 1500) |
| 394 |
end if |
| 395 |
if (iprint .gt. 0) write(iprint, 1600) emax, imax |
| 396 |
if (cheap) go to 900 |
| 397 |
|
| 398 |
* ---------------------------------------------------- |
| 399 |
* Proceed with the verification of columns j1 thru j2. |
| 400 |
* ---------------------------------------------------- |
| 401 |
if (iprint .gt. 0) write(iprint, 2000) |
| 402 |
l = lsave |
| 403 |
nwrong = 0 |
| 404 |
ngood = 0 |
| 405 |
emax = - one |
| 406 |
jmax = 0 |
| 407 |
|
| 408 |
do 200 j = j1, j2 |
| 409 |
|
| 410 |
* See if there are any known gradients in this column. |
| 411 |
|
| 412 |
k1 = ka(j) |
| 413 |
k2 = ka(j+1) - 1 |
| 414 |
lsave = l |
| 415 |
|
| 416 |
do 120 k = k1, k2 |
| 417 |
ir = ha(k) |
| 418 |
if (ir .gt. m) go to 200 |
| 419 |
l = l + 1 |
| 420 |
if (g(l) .ne. gdummy) go to 140 |
| 421 |
120 continue |
| 422 |
go to 200 |
| 423 |
|
| 424 |
* Found one. |
| 425 |
|
| 426 |
140 xj = x(j) |
| 427 |
dx = difint(1) * (one + abs( xj )) |
| 428 |
if (bl(j) .lt. bu(j) .and. xj .ge. bu(j)) dx = -dx |
| 429 |
x(j) = xj + dx |
| 430 |
call m6fcon( 0, m, n, njac, f2, g2, |
| 431 |
$ ne, nka, ha, ka, |
| 432 |
$ x, z, nwcore ) |
| 433 |
if (ierr .ne. 0) go to 900 |
| 434 |
|
| 435 |
* Check nonzeros in the jth column of the Jacobian. |
| 436 |
* Don't bother printing a line if it looks like an exact zero. |
| 437 |
|
| 438 |
l = lsave |
| 439 |
first = .true. |
| 440 |
|
| 441 |
do 160 k = k1, k2 |
| 442 |
ir = ha(k) |
| 443 |
if (ir .gt. m) go to 180 |
| 444 |
l = l + 1 |
| 445 |
gi = g(l) |
| 446 |
if (gi .eq. gdummy) go to 160 |
| 447 |
gdiff = (f2(ir) - f(ir))/dx |
| 448 |
err = abs( gdiff - gi ) / (one + abs( gi )) |
| 449 |
|
| 450 |
if (emax .lt. err) then |
| 451 |
emax = err |
| 452 |
imax = ir |
| 453 |
jmax = j |
| 454 |
end if |
| 455 |
|
| 456 |
key = lgood |
| 457 |
if (err .gt. eps5) key = lbad |
| 458 |
if (key .eq. lbad) nwrong = nwrong + 1 |
| 459 |
if (key .eq.lgood) ngood = ngood + 1 |
| 460 |
if (abs( gi ) + err .le. eps0) go to 160 |
| 461 |
if (iprint .gt. 0) then |
| 462 |
if (first) then |
| 463 |
write(iprint, 2100) j,xj,dx,l,ir,gi,gdiff,key |
| 464 |
else |
| 465 |
write(iprint, 2200) l,ir,gi,gdiff,key |
| 466 |
end if |
| 467 |
end if |
| 468 |
first = .false. |
| 469 |
160 continue |
| 470 |
|
| 471 |
180 x(j) = xj |
| 472 |
200 continue |
| 473 |
|
| 474 |
if (iprint .gt. 0) then |
| 475 |
if (nwrong .eq. 0) then |
| 476 |
write(iprint, 2500) ngood ,j1,j2 |
| 477 |
else |
| 478 |
write(iprint, 2600) nwrong,j1,j2 |
| 479 |
end if |
| 480 |
write(iprint, 2700) emax,imax,jmax |
| 481 |
end if |
| 482 |
|
| 483 |
if (emax .ge. one) then |
| 484 |
|
| 485 |
* Bad gradients in funcon. |
| 486 |
|
| 487 |
ierr = 8 |
| 488 |
call m1envt( 1 ) |
| 489 |
if (iprint .gt. 0) write(iprint, 3800) |
| 490 |
if (isumm .gt. 0) write(isumm , 3800) |
| 491 |
end if |
| 492 |
|
| 493 |
* Exit. |
| 494 |
|
| 495 |
900 lscale = lssave |
| 496 |
return |
| 497 |
|
| 498 |
1000 format(/// ' Verification of constraint gradients', |
| 499 |
$ ' returned by subroutine funcon.') |
| 500 |
1100 format(/ ' Cheap test on funcon...') |
| 501 |
1400 format( ' The Jacobian seems to be OK.') |
| 502 |
1500 format( ' XXX The Jacobian seems to be incorrect.') |
| 503 |
1600 format( ' The largest discrepancy was', 1p, e12.2, |
| 504 |
$ ' in constraint', i6) |
| 505 |
2000 format(// ' Column x(j) dx(j)', 3x, |
| 506 |
$ ' Element no. Row Jacobian value Difference approxn') |
| 507 |
2100 format(/ i7, 1p, e16.8, e10.2, 2i10, 2e18.8, 2x, a4) |
| 508 |
2200 format( 33x, 2i10, 1pe18.8, e18.8, 2x, a4) |
| 509 |
2500 format(/ i7, ' Jacobian elements in cols ', i6, ' thru', i6, |
| 510 |
$ ' seem to be OK.') |
| 511 |
2600 format(/ ' XXX There seem to be', i6, |
| 512 |
$ ' incorrect Jacobian elements in cols', i6, ' thru', i6) |
| 513 |
2700 format(/ ' XXX The largest relative error was', 1p, e12.2, |
| 514 |
$ ' in row', i6, ', column', i6 /) |
| 515 |
3800 format(// ' EXIT -- subroutine funcon appears to be', |
| 516 |
$ ' giving incorrect gradients') |
| 517 |
|
| 518 |
* end of m8chkj |
| 519 |
end |
| 520 |
|
| 521 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 522 |
|
| 523 |
subroutine m8prtj( n, nb, nncon, nnjac, lprint, majits, nlag,nscl, |
| 524 |
$ ne, nka, a, ha, ka, hs, |
| 525 |
$ ascale, bl, bu, fcon, xlam, xn ) |
| 526 |
|
| 527 |
implicit double precision (a-h,o-z) |
| 528 |
integer*4 ha(ne), hs(nnjac) |
| 529 |
integer ka(nka) |
| 530 |
double precision a(ne), ascale(nscl), bl(nb), bu(nb), |
| 531 |
$ fcon(nncon), xlam(nncon), xn(nb) |
| 532 |
|
| 533 |
* ------------------------------------------------------------------ |
| 534 |
* m8prtj prints x, lambda, f(x) and/or the Jacobian |
| 535 |
* at the start of each major iteration. |
| 536 |
* |
| 537 |
* 08 Apr 1992: Internal values of hs(*) allow more values for key. |
| 538 |
* ------------------------------------------------------------------ |
| 539 |
|
| 540 |
common /m1file/ iread,iprint,isumm |
| 541 |
common /m3scal/ sclobj,scltol,lscale |
| 542 |
|
| 543 |
intrinsic min, mod |
| 544 |
|
| 545 |
logical prtx, prtl, prtf, prtj, scaled |
| 546 |
character*2 key(-1:4) |
| 547 |
data key/'FR', 'LO', 'UP', 'SB', 'BS', 'FX'/ |
| 548 |
|
| 549 |
if (iprint .le. 0) return |
| 550 |
|
| 551 |
* Unscale everything if necessary. |
| 552 |
|
| 553 |
scaled = lscale .ge. 2 |
| 554 |
if (scaled) call m2scla( 2, nncon, n, nb, ne, nka, |
| 555 |
$ ha, ka, a, ascale, bl, bu, xlam, xn ) |
| 556 |
if (scaled) call ddscl ( nncon, ascale(n+1), 1, fcon, 1 ) |
| 557 |
|
| 558 |
nxn = min( nnjac, 5 ) |
| 559 |
nlam = min( nncon, 5 ) |
| 560 |
l = lprint/10 |
| 561 |
prtx = mod( l,10 ) .gt. 0 |
| 562 |
l = l/10 |
| 563 |
prtl = mod( l,10 ) .gt. 0 .and. majits .gt. 1 |
| 564 |
l = l/10 |
| 565 |
prtf = mod( l,10 ) .gt. 0 |
| 566 |
l = l/10 |
| 567 |
prtj = mod( l,10 ) .gt. 0 |
| 568 |
|
| 569 |
if ( prtx ) write(iprint, 1100) (xn(j), j=1,nnjac) |
| 570 |
if ( prtl ) write(iprint, 1200) (xlam(i), i=1,nncon) |
| 571 |
if ( prtf ) write(iprint, 1300) (fcon(i), i=1,nncon) |
| 572 |
if ( prtj ) then |
| 573 |
write(iprint, 1400) |
| 574 |
do 100 j = 1, nnjac |
| 575 |
k1 = ka(j) |
| 576 |
k2 = ka(j+1) - 1 |
| 577 |
do 50 k = k1, k2 |
| 578 |
ir = ha(k) |
| 579 |
if (ir .gt. nncon) go to 60 |
| 580 |
50 continue |
| 581 |
k = k2 + 1 |
| 582 |
60 k2 = k - 1 |
| 583 |
l = hs(j) |
| 584 |
write(iprint, 1410) j,xn(j),key(l), (ha(k),a(k), k=k1,k2) |
| 585 |
100 continue |
| 586 |
end if |
| 587 |
|
| 588 |
* Rescale if necessary. |
| 589 |
|
| 590 |
if (scaled) call m2scla( 1, nncon, n, nb, ne, nka, |
| 591 |
$ ha, ka, a, ascale, bl, bu, xlam, xn ) |
| 592 |
if (scaled) call dddiv ( nncon, ascale(n+1), 1, fcon, 1 ) |
| 593 |
return |
| 594 |
|
| 595 |
1100 format(/ ' Jacobian variables' |
| 596 |
$ / ' ------------------' / 1p, (5e16.7)) |
| 597 |
1200 format(/ ' Multiplier estimates' |
| 598 |
$ / ' --------------------' / 1p, (5e16.7)) |
| 599 |
1300 format(/ ' Constraint functions' |
| 600 |
$ / ' --------------------' / 1p, (5e16.7)) |
| 601 |
1400 format(/ ' x and Jacobian' / ' ----------------') |
| 602 |
1410 format(i6, 1p, e13.5, 1x, a2, 1x, 4(i9, e13.5) |
| 603 |
$ / (22x, 4(i9, e13.5))) |
| 604 |
|
| 605 |
* end of m8prtj |
| 606 |
end |
| 607 |
|
| 608 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 609 |
|
| 610 |
subroutine m8sclj( mode, nncon, nn, ng, n, nb, ne, nka, |
| 611 |
$ ascale, ha, ka, g ) |
| 612 |
|
| 613 |
implicit double precision (a-h,o-z) |
| 614 |
integer*4 ha(ne) |
| 615 |
integer ka(nka) |
| 616 |
double precision ascale(nb), g(ng) |
| 617 |
|
| 618 |
* ------------------------------------------------------------------ |
| 619 |
* If mode = 1, m8sclj scales the objective gradient. |
| 620 |
* If mode = 2, m8sclj scales the Jacobian. |
| 621 |
* nn, ng, g are nnobj, nnobj, gobj or nnjac, njac, gcon resp. |
| 622 |
* |
| 623 |
* m8sclj is called by m6fobj and m6fcon only if modefg = 2. |
| 624 |
* Hence, it is used to scale known gradient elements (if any), |
| 625 |
* but is not called when missing gradients are being estimated |
| 626 |
* by m6dobj or m6dcon. |
| 627 |
* ------------------------------------------------------------------ |
| 628 |
|
| 629 |
common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 630 |
|
| 631 |
if (knowng(mode) .eq. 0) return |
| 632 |
|
| 633 |
if (mode .eq. 1) then |
| 634 |
* --------------------------------------------------------------- |
| 635 |
* Scale known objective gradients. |
| 636 |
* --------------------------------------------------------------- |
| 637 |
do 100 j = 1, nn |
| 638 |
grad = g(j) |
| 639 |
if (grad .ne. gdummy) g(j) = grad * ascale(j) |
| 640 |
100 continue |
| 641 |
|
| 642 |
else |
| 643 |
* --------------------------------------------------------------- |
| 644 |
* Scale known Jacobian elements. |
| 645 |
* --------------------------------------------------------------- |
| 646 |
l = 0 |
| 647 |
|
| 648 |
do 300 j = 1, nn |
| 649 |
k1 = ka(j) |
| 650 |
k2 = ka(j+1) - 1 |
| 651 |
cscale = ascale(j) |
| 652 |
|
| 653 |
do 250 k = k1, k2 |
| 654 |
ir = ha(k) |
| 655 |
if (ir .gt. nncon) go to 300 |
| 656 |
l = l + 1 |
| 657 |
grad = g(l) |
| 658 |
if (grad .ne. gdummy) |
| 659 |
$ g(l) = grad * cscale / ascale(n + ir) |
| 660 |
250 continue |
| 661 |
300 continue |
| 662 |
end if |
| 663 |
|
| 664 |
* end of m8sclj |
| 665 |
end |
| 666 |
|
| 667 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 668 |
|
| 669 |
subroutine m8setj( lcstat, m, n, nb, nn, nncon, nnjac, njac, nscl, |
| 670 |
$ lcrash, ns, |
| 671 |
$ nconvg, nomove, nreduc, optsub, convgd, |
| 672 |
$ ne, nka, a, ha, ka, hs, |
| 673 |
$ ascale, bl, bu, |
| 674 |
$ blslk, buslk, fcon, fcon2, |
| 675 |
$ fdif , fold , gcon, gcon2, xlam, |
| 676 |
$ rhs , xdif , xold, pi, xn, y, z, nwcore ) |
| 677 |
|
| 678 |
implicit double precision (a-h,o-z) |
| 679 |
logical optsub, convgd |
| 680 |
integer*4 ha(ne), hs(nb) |
| 681 |
integer ka(nka) |
| 682 |
double precision a(ne), ascale(nscl), bl(nb), bu(nb) |
| 683 |
double precision blslk(nncon), buslk(nncon), |
| 684 |
$ fcon (nncon), fcon2(nncon), |
| 685 |
$ fdif (nncon), fold (nncon), |
| 686 |
$ gcon (njac ), gcon2(njac ), |
| 687 |
$ xlam (nncon), rhs (nncon), xdif(nn), xold(nn), |
| 688 |
$ pi(m), xn(nb), y(m), z(nwcore) |
| 689 |
|
| 690 |
* ------------------------------------------------------------------ |
| 691 |
* m8setj sets up the linearly constrained subproblem for the |
| 692 |
* next major iteration of the augmented Lagrangian algorithm. |
| 693 |
* |
| 694 |
* lcrash is an input parameter. If lcrash = 5, the nonlinear |
| 695 |
* constraints have been relaxed so far. On major itn 2, |
| 696 |
* Crash is called to find a basis including them. |
| 697 |
* |
| 698 |
* nconvg counts no. of times optimality test has been satisfied. |
| 699 |
* nreduc counts no. of times penpar has been reduced. |
| 700 |
* optsub (input ) says if the previous subproblem was optimal. |
| 701 |
* convgd (output) says if it is time to stop. |
| 702 |
* |
| 703 |
* 06 Mar 1992: The first major iteration now relaxes nonlinear rows |
| 704 |
* and terminates as soon as the linear constraints are |
| 705 |
* satisfied. xold is not defined until Major 2. |
| 706 |
* 08 Apr 1992: Internal values of hs(*) simplify setting xlam2. |
| 707 |
* 22 May 1992: penpar reduced for more majors before setting to zero. |
| 708 |
* Stops Steinke2 (and maybe OTIS?) from suddenly |
| 709 |
* getting a huge x and lambda. |
| 710 |
* 04 Jun 1992: lcrash added. |
| 711 |
* 03 Sep 1992: Penalty parameter is now a multiple of the |
| 712 |
* default value 100/nncon. |
| 713 |
* ------------------------------------------------------------------ |
| 714 |
|
| 715 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 716 |
common /m1file/ iread,iprint,isumm |
| 717 |
common /m2lu3 / lenl,lenu,ncp,lrow,lcol |
| 718 |
common /m2parm/ dparm(30),iparm(30) |
| 719 |
common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu , |
| 720 |
$ lhrtyp,lhs ,lkb |
| 721 |
common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj |
| 722 |
common /m5log1/ idebug,ierr,lprint |
| 723 |
common /m5log2/ jq1,jq2,jr1,jr2,lines1,lines2 |
| 724 |
logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed |
| 725 |
common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed |
| 726 |
common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi |
| 727 |
common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg |
| 728 |
common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor |
| 729 |
common /m8al2 / radius,rhsmod,modpen,modrhs |
| 730 |
*- common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 731 |
common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2 |
| 732 |
common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave |
| 733 |
common /cycle2/ objtru,suminf,numinf |
| 734 |
|
| 735 |
intrinsic abs, max, min |
| 736 |
|
| 737 |
parameter ( zero = 0.0d+0, one = 1.0d+0, |
| 738 |
$ biglam = 1.0d+5, dmaxlm = 1.0d+8 ) |
| 739 |
|
| 740 |
logical gotjac, feasbl, major1, major2, phead , zerlam |
| 741 |
double precision objold |
| 742 |
save objold |
| 743 |
|
| 744 |
character*8 line |
| 745 |
character*1 lstate(4), label |
| 746 |
data line /'--------'/ |
| 747 |
data lstate /' ','I','U','T'/ |
| 748 |
|
| 749 |
rhsmod = zero |
| 750 |
modrhs = 0 |
| 751 |
bplus = 0.9*plinfy |
| 752 |
bminus = - bplus |
| 753 |
dlmax = zero |
| 754 |
dlrel = zero |
| 755 |
dxrel = zero |
| 756 |
step = one |
| 757 |
feasbl = ninf .eq. 0 |
| 758 |
zerlam = .true. |
| 759 |
major1 = majits .eq. 1 |
| 760 |
major2 = majits .eq. 2 |
| 761 |
|
| 762 |
* Output heading for terse log. |
| 763 |
|
| 764 |
phead = newhed .or. |
| 765 |
$ majits .gt. 2 .and. mod( majits, 20 ) .eq. 2 |
| 766 |
if (prnt0 .and. phead) write(iprint, 1100) |
| 767 |
if ( phead) newhed = .false. |
| 768 |
phead = major1 .or. |
| 769 |
$ majits .gt. 2 .and. mod( majits, 10 ) .eq. 2 |
| 770 |
if (summ0 .and. phead) write(isumm , 1110) |
| 771 |
|
| 772 |
* Output heading for detailed log. |
| 773 |
|
| 774 |
if (prnt1 .and. .not. major1) call m1page( 0 ) |
| 775 |
if (prnt1) write(iprint, 2010) (line, j=1,16), majits, penpar |
| 776 |
if (summ1) write(isumm , 2020) (line, j=1, 8), majits, penpar |
| 777 |
|
| 778 |
* ------------------------------------------------------------------ |
| 779 |
* Beginning of first major iteration. |
| 780 |
* misolv has already called m8augl( 2,... ) to make sure that |
| 781 |
* all nonlinear variables are inside their bounds. |
| 782 |
* ------------------------------------------------------------------ |
| 783 |
if (major1) then |
| 784 |
|
| 785 |
* For interest, get the initial constraint violation. |
| 786 |
|
| 787 |
call m8viol( n, nb, nncon, |
| 788 |
$ ne, nka, a, ha, ka, |
| 789 |
$ bl, bu, fcon, xn, y, z, nwcore ) |
| 790 |
|
| 791 |
* Save the nonlinear constraint bounds and relax them. |
| 792 |
* (The first major will just satisfy the linear constraints.) |
| 793 |
* Then set internal values for hs. |
| 794 |
|
| 795 |
call m8augl( 3, m, n, nb, ns, inform, |
| 796 |
$ ne, nka, a, ha, ka, |
| 797 |
$ hs, bl, bu, xn, z, nwcore ) |
| 798 |
call m5hs ( 'Internal', nb, bl, bu, hs, xn ) |
| 799 |
|
| 800 |
* xlam has been initialized from pi in m4getb. |
| 801 |
* Initialize xold and fold, so they are defined when m5frmc |
| 802 |
* evaluates the augmented Lagrangian at the first feasible point. |
| 803 |
|
| 804 |
call dcopy ( nn , xn , 1, xold, 1 ) |
| 805 |
call dcopy ( nncon, fcon, 1, fold, 1 ) |
| 806 |
go to 500 |
| 807 |
end if |
| 808 |
|
| 809 |
* ------------------------------------------------------------------ |
| 810 |
* Beginning of later major iterations. |
| 811 |
* ------------------------------------------------------------------ |
| 812 |
|
| 813 |
* Restore the bounds on the nonlinear constraints. |
| 814 |
|
| 815 |
call m8augl( 6, m, n, nb, ns, inform, |
| 816 |
$ ne, nka, a, ha, ka, |
| 817 |
$ hs, bl, bu, xn, z, nwcore ) |
| 818 |
|
| 819 |
* Make sure all variables are inside their bounds. |
| 820 |
* Then set internal values for hs. |
| 821 |
|
| 822 |
call m8augl( 2, m, n, nb, ns, inform, |
| 823 |
$ ne, nka, a, ha, ka, |
| 824 |
$ hs, bl, bu, xn, z, nwcore ) |
| 825 |
call m5hs ( 'Internal', nb, bl, bu, hs, xn ) |
| 826 |
|
| 827 |
if (major2) then |
| 828 |
|
| 829 |
* Initialize xold. |
| 830 |
* xlam still contains pi from m4getb. Check it is OK. |
| 831 |
* fold is set later. |
| 832 |
|
| 833 |
call dcopy ( nn, xn, 1, xold, 1 ) |
| 834 |
step = one |
| 835 |
|
| 836 |
do 120 i = 1, nncon |
| 837 |
j = n + i |
| 838 |
js = hs(j) |
| 839 |
py = xlam(i) |
| 840 |
xlam2 = py |
| 841 |
if (js .eq. 0 .and. py .gt. zero) xlam2 = zero |
| 842 |
if (js .eq. 1 .and. py .lt. zero) xlam2 = zero |
| 843 |
xlam2 = min( xlam2, dmaxlm ) |
| 844 |
xlam2 = max( xlam2, (- dmaxlm) ) |
| 845 |
xlam(i) = xlam2 |
| 846 |
120 continue |
| 847 |
|
| 848 |
go to 400 |
| 849 |
end if |
| 850 |
|
| 851 |
* ------------------------------------------------------------------ |
| 852 |
* Major 3 and later. |
| 853 |
* If xn is feasible, reset the multiplier estimates. |
| 854 |
* Since the previous subproblem could have been stopped early, |
| 855 |
* the pi's for inequality constraints must be checked. |
| 856 |
* If constraint i is active and pi(i) has the wrong sign, |
| 857 |
* we just use zero. |
| 858 |
* 22 May 1992: If the subproblem is optimal and a slack is |
| 859 |
* superbasic, the pi(i) should again be zero. |
| 860 |
* ------------------------------------------------------------------ |
| 861 |
if ( feasbl ) then |
| 862 |
do 150 i = 1, nncon |
| 863 |
j = n + i |
| 864 |
js = hs(j) |
| 865 |
py = pi(i) |
| 866 |
xlam2 = py |
| 867 |
if (js .eq. 0 .and. py .gt. zero) xlam2 = zero |
| 868 |
if (js .eq. 1 .and. py .lt. zero) xlam2 = zero |
| 869 |
if (js .eq. 2 .and. optsub ) xlam2 = zero |
| 870 |
xlam2 = min( xlam2, dmaxlm ) |
| 871 |
xlam2 = max( xlam2, (- dmaxlm) ) |
| 872 |
pi(i) = xlam2 |
| 873 |
y(i) = xlam2 - xlam(i) |
| 874 |
150 continue |
| 875 |
|
| 876 |
* Find the relative change in xlam. |
| 877 |
* If xlam is very small, it is probably still zero (never set). |
| 878 |
* We treat it specially to avoid taking a small step below. |
| 879 |
|
| 880 |
imax = idamax( nncon, y , 1 ) |
| 881 |
dlmax = abs ( y(imax) ) |
| 882 |
dlnorm = dasum ( nncon, y , 1 ) |
| 883 |
xlnorm = dasum ( nncon, xlam, 1 ) |
| 884 |
zerlam = xlnorm .le. eps0 |
| 885 |
dlrel = dlnorm / (one + xlnorm) |
| 886 |
end if |
| 887 |
|
| 888 |
* ------------------------------------------------------------------ |
| 889 |
* Find the relative change in the Jacobian variables. |
| 890 |
* Note that xdif and xold allow for all nonlinear variables, |
| 891 |
* in case "step" below turns out not to be one. |
| 892 |
* ------------------------------------------------------------------ |
| 893 |
do 220 j = 1, nn |
| 894 |
xdif(j) = xn(j) - xold(j) |
| 895 |
220 continue |
| 896 |
|
| 897 |
imax = idamax( nnjac, xdif, 1 ) |
| 898 |
dxmax = abs ( xdif(imax) ) |
| 899 |
dxrel = dasum ( nnjac, xdif, 1 ) / (one + dasum ( nnjac,xold,1 )) |
| 900 |
if (prnt1) write(iprint, 2400) dxmax, dxrel, dlmax, dlrel |
| 901 |
if (summ1) write(isumm , 2410) dxmax, dlmax |
| 902 |
|
| 903 |
* ------------------------------------------------------------------ |
| 904 |
* Determine a step to be used as follows: |
| 905 |
* xlam = xlam + step*y, xn = xold + step*xdif. |
| 906 |
* |
| 907 |
* damp = 2.0 (for example) prevents x or lambda from changing |
| 908 |
* by more than 200 per cent. This is an exceedingly crude attempt |
| 909 |
* to avoid trouble when the new subproblem solution xn is wildly |
| 910 |
* different from the previous solution xold. |
| 911 |
* |
| 912 |
* Note that the current active set is retained even if step ends |
| 913 |
* up being less than 1.0. Some nonbasic variables may lie strictly |
| 914 |
* between their bounds (which is ok now that xn retains all values). |
| 915 |
* ------------------------------------------------------------------ |
| 916 |
damp = dparm(4) |
| 917 |
step = one |
| 918 |
if (.not. zerlam) then |
| 919 |
step = min( step, damp/(dlrel + eps) ) |
| 920 |
end if |
| 921 |
step = min( step, damp/(dxrel + eps) ) |
| 922 |
|
| 923 |
if (step .ge. 0.99) then |
| 924 |
* ===================== |
| 925 |
* A unit step looks ok. |
| 926 |
* ===================== |
| 927 |
step = one |
| 928 |
if (feasbl) call dcopy ( nncon, pi, 1, xlam, 1 ) |
| 929 |
call dcopy ( nn, xn, 1, xold, 1 ) |
| 930 |
|
| 931 |
* see if xn is the same as xold. |
| 932 |
|
| 933 |
if (dxrel .le. eps0) then |
| 934 |
nomove = nomove + 1 |
| 935 |
if (prnt1) write(iprint, 2100) |
| 936 |
if (summ1) write(isumm , 2100) |
| 937 |
go to 450 |
| 938 |
end if |
| 939 |
else |
| 940 |
* ==================== |
| 941 |
* Take a shorter step. |
| 942 |
* ==================== |
| 943 |
if (feasbl) call daxpy ( nncon, step, y, 1, xlam, 1 ) |
| 944 |
call daxpy ( nn, step, xdif, 1, xold, 1 ) |
| 945 |
call dcopy ( nn, xold, 1, xn, 1 ) |
| 946 |
|
| 947 |
if (prnt1) write(iprint, 2500) step |
| 948 |
if (summ1) write(isumm , 2500) step |
| 949 |
end if |
| 950 |
|
| 951 |
* ------------------------------------------------------------------ |
| 952 |
* Evaluate the nonlinear constraints and the Jacobian |
| 953 |
* and copy the Jacobian into A to form the next linearization. |
| 954 |
* ------------------------------------------------------------------ |
| 955 |
400 nomove = 0 |
| 956 |
gotjac = feasbl .and. nlag .eq. 1 .and. step .eq. one |
| 957 |
|
| 958 |
call m8ajac( gotjac, nncon, nnjac, njac, |
| 959 |
$ ne, nka, a, ha, ka, |
| 960 |
$ fcon, fcon2, gcon, gcon2, xn, y, z, nwcore ) |
| 961 |
if (ierr .ne. 0) go to 900 |
| 962 |
|
| 963 |
* Save the function values in fold. |
| 964 |
|
| 965 |
call dcopy ( nncon, fcon, 1, fold, 1 ) |
| 966 |
|
| 967 |
* Find how well the new xn satisfies the nonlinear constraints. |
| 968 |
|
| 969 |
450 call m8viol( n, nb, nncon, |
| 970 |
$ ne, nka, a, ha, ka, |
| 971 |
$ bl, bu, fcon, xn, y, z, nwcore ) |
| 972 |
|
| 973 |
* ------------------------------------------------------------------ |
| 974 |
* Print a line. |
| 975 |
* ------------------------------------------------------------------ |
| 976 |
500 label = lstate(lcstat + 1) |
| 977 |
lu = lenl + lenu |
| 978 |
if (prnt0) write(iprint, 1200) majits - 1, minits, label, itn, |
| 979 |
$ ninf, objtru, vimax, rgnorm, nssave, |
| 980 |
$ nfcon(1), dxrel, dlrel, step, lu, penpar |
| 981 |
C write(isumm , 1210) majits - 1, minits, label, |
| 982 |
if (summ0) write(isumm , 1210) majits - 1, minits, label, |
| 983 |
$ objtru, vimax, rgnorm, nssave, |
| 984 |
$ nfcon(1), dxrel, dlrel, step |
| 985 |
|
| 986 |
if (prnt1 .and. (major1 .or. major2)) write(iprint, 1000) |
| 987 |
if (prnt1) write(iprint, 2600) vimax, virel |
| 988 |
if (summ1) write(isumm , 2700) vimax |
| 989 |
|
| 990 |
* ------------------------------------------------------------------ |
| 991 |
* Compute rhs = fcon - Jacobian*xn. |
| 992 |
* This is minus the required rhs for the new subproblem. |
| 993 |
* ------------------------------------------------------------------ |
| 994 |
if (lprint .gt. 1) |
| 995 |
$call m8prtj( n, nb, nncon, nnjac, lprint, majits, nlag, nscl, |
| 996 |
$ ne, nka, a, ha, ka, hs, |
| 997 |
$ ascale, bl, bu, fcon, xlam, xn ) |
| 998 |
|
| 999 |
call dcopy ( nncon, fcon, 1, rhs, 1 ) |
| 1000 |
call m2aprd( 7, xn, nnjac, rhs, nncon, |
| 1001 |
$ ne, nka, a, ha, ka, |
| 1002 |
$ z, nwcore ) |
| 1003 |
call dscal ( nncon, (- one), rhs, 1 ) |
| 1004 |
|
| 1005 |
* ------------------------------------------------------------------ |
| 1006 |
* Do Crash on nonlinear rows at start of Major 2 |
| 1007 |
* unless a basis is known already. |
| 1008 |
* ------------------------------------------------------------------ |
| 1009 |
if (major2 .and. lcrash .eq. 5) then |
| 1010 |
|
| 1011 |
* m8viol has set y to be the exact slacks on nonlinear rows. |
| 1012 |
* Copy these into xn(n+i) and make sure they are feasible. |
| 1013 |
* m2crsh uses them to decide which slacks to grab for the basis. |
| 1014 |
|
| 1015 |
do 550 i = 1, nncon |
| 1016 |
j = n + i |
| 1017 |
xn(j) = max( y(i), bl(j) ) |
| 1018 |
xn(j) = min( xn(j), bu(j) ) |
| 1019 |
550 continue |
| 1020 |
|
| 1021 |
* Set row types in hrtype. |
| 1022 |
* m2crsh uses kb as workspace. It may alter xn(n+i) for |
| 1023 |
* nonlinear slacks. |
| 1024 |
* Set internal values for hs again to match xn. |
| 1025 |
|
| 1026 |
call m2amat( 2, m, n, nb, |
| 1027 |
$ ne, nka, a, ha, ka, |
| 1028 |
$ bl, bu, z(lhrtyp) ) |
| 1029 |
call m2crsh( lcrash, m , n , nb , nn, |
| 1030 |
$ ne , nka , a , ha , ka, |
| 1031 |
$ z(lkb), hs , z(lhrtyp), bl , bu, |
| 1032 |
$ xn , z , nwcore ) |
| 1033 |
lcrash = 0 |
| 1034 |
call m5hs ( 'Internal', nb, bl, bu, hs, xn ) |
| 1035 |
end if |
| 1036 |
|
| 1037 |
* ------------------------------------------------------------------ |
| 1038 |
* Test for optimality. |
| 1039 |
* For safety, we require the convergence test to be satisfied |
| 1040 |
* for 2 subproblems in a row. |
| 1041 |
* ------------------------------------------------------------------ |
| 1042 |
if (major1 .or. major2) then |
| 1043 |
objold = objtru |
| 1044 |
convgd = .false. |
| 1045 |
nconvg = 0 |
| 1046 |
else |
| 1047 |
dfrel = abs( objtru - objold ) / (one + objold) |
| 1048 |
objold = objtru |
| 1049 |
|
| 1050 |
convgd = dxrel .le. 0.1 .and. dlrel .le. 0.1 .and. |
| 1051 |
$ virel .le. rowtol .and. optsub .and. |
| 1052 |
$ dfrel .le. 0.001 |
| 1053 |
nconvg = nconvg + 1 |
| 1054 |
if (.not. convgd) nconvg = 0 |
| 1055 |
convgd = nconvg .ge. 2 |
| 1056 |
|
| 1057 |
* Change from Partial Completion to Full Completion if |
| 1058 |
* the iterations seem to be converging. |
| 1059 |
* (This may be begging the question! More theory needed.) |
| 1060 |
|
| 1061 |
if (ncom .ne. 1) then |
| 1062 |
if (dxrel .le. 0.1 .and. dlrel .le. 0.1 .and. |
| 1063 |
$ virel .le. 0.1 .and. optsub ) then |
| 1064 |
ncom = 1 |
| 1065 |
if (summ1) write(isumm , 2800) |
| 1066 |
if (prnt1) write(iprint, 2800) |
| 1067 |
end if |
| 1068 |
end if |
| 1069 |
|
| 1070 |
* ============================= |
| 1071 |
* Adjust the penalty parameter. |
| 1072 |
* ============================= |
| 1073 |
if (feasbl .and. optsub) then |
| 1074 |
|
| 1075 |
if (virel .le. radius .and. dlrel .le. radius) then |
| 1076 |
if (penpar .gt. zero) then |
| 1077 |
nreduc = nreduc + 1 |
| 1078 |
penpar = penpar / 10.0 |
| 1079 |
|
| 1080 |
* The next line caused trouble on Steinke2 |
| 1081 |
* (suddenly a huge search direction). |
| 1082 |
* For safety we should always reduce penpar gradually. |
| 1083 |
*------> if (nconvg .ge. 1 .or. nreduc .ge. 5) penpar = zero |
| 1084 |
|
| 1085 |
if (penpar .le. 1.0d-6 .or. nreduc .ge. 10) then |
| 1086 |
penpar = zero |
| 1087 |
end if |
| 1088 |
end if |
| 1089 |
|
| 1090 |
else if (majits .ge. 5 .and. dlrel .ge. 10.0) then |
| 1091 |
nreduc = 0 |
| 1092 |
oldpen = penpar |
| 1093 |
if (penpar .le. zero) penpar = 1.0d+0 |
| 1094 |
*-04 Sep 1992: Raising the penalty parameter (below) seems to cause |
| 1095 |
*- more trouble than it fixes. Give it up for now. |
| 1096 |
*- fac = 1.5 |
| 1097 |
*- if (dlrel .ge. biglam) fac = 2.0 |
| 1098 |
*- penpar = penpar * fac |
| 1099 |
*- penmax = 1.0d+3 |
| 1100 |
*- penpar = min( penpar, penmax ) |
| 1101 |
if (iprint .gt. 0 .and. penpar .gt. oldpen) then |
| 1102 |
write(iprint, 2900) penpar |
| 1103 |
end if |
| 1104 |
end if |
| 1105 |
end if |
| 1106 |
end if |
| 1107 |
|
| 1108 |
* ------------------------------------------------------------------ |
| 1109 |
* Exit. |
| 1110 |
* ------------------------------------------------------------------ |
| 1111 |
900 return |
| 1112 |
|
| 1113 |
1000 format(' ') |
| 1114 |
1100 format(/ ' Major minor total ninf true objective viol', |
| 1115 |
$ ' rg nsb ncon xchange lchange step', |
| 1116 |
$ ' LU penalty') |
| 1117 |
1110 format(/ ' Major minor true objective viol', |
| 1118 |
$ ' rg nsb ncon xchange lchange step') |
| 1119 |
C1210 FORMAT(2I6, A1, I5, 1PE16.8, 2E8.1, I5, I6, 2E8.1) |
| 1120 |
1210 format(1p, 2i6, a1, e16.8, 2e8.1, i5, i7, 3e8.1) |
| 1121 |
|
| 1122 |
1200 format(1p, 2i6, a1, i7, i5, e16.8, 2e8.1, i5, i7, 3e8.1, i7, e8.1) |
| 1123 |
2010 format( 1x, 16a8 / ' Start of major itn', i4, |
| 1124 |
$ 10x, ' Penalty parameter =', 1p, e12.2) |
| 1125 |
2020 format(/ 1x, 8a8 / ' Start of major itn', i4, |
| 1126 |
$ 10x, ' Penalty parameter =', 1p, e12.2) |
| 1127 |
2100 format(/ ' Jacobian variables have not changed --', |
| 1128 |
$ ' old linearization kept') |
| 1129 |
2400 format(/ ' Maximum change in Jacobian vars =', 1p, e12.4, |
| 1130 |
$ 4x, ' ( =', e11.4, ' normalized)' |
| 1131 |
$ / ' Maximum change in multipliers =', e12.4, |
| 1132 |
$ 4x, ' ( =', e11.4, ' normalized)' ) |
| 1133 |
2410 format(' Change in Jacobn vars =', 1p, e12.4 |
| 1134 |
$ / ' Change in multipliers =', e12.4) |
| 1135 |
2500 format( ' Step size for x and lamda =', f12.5) |
| 1136 |
2600 format( ' Maximum constraint violation =', 1p, e12.4, |
| 1137 |
$ 4x, ' ( =', e11.4, ' normalized)' ) |
| 1138 |
2700 format(' Constraint violation =', 1p, e12.4) |
| 1139 |
2800 format(/ ' Completion Full requested as from now') |
| 1140 |
2900 format(/ ' Penalty parameter increased to', 1p, e12.2 /) |
| 1141 |
|
| 1142 |
* end of m8setj |
| 1143 |
end |
| 1144 |
|
| 1145 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1146 |
|
| 1147 |
subroutine m8viol( n, nb, nncon, |
| 1148 |
$ ne, nka, a, ha, ka, |
| 1149 |
$ bl, bu, fcon, xn, y, z, nwcore ) |
| 1150 |
|
| 1151 |
implicit double precision (a-h,o-z) |
| 1152 |
integer*4 ha(ne) |
| 1153 |
integer ka(nka) |
| 1154 |
double precision a(ne), bl(nb), bu(nb), fcon(nncon), |
| 1155 |
$ xn(nb), y(nncon), z(nwcore) |
| 1156 |
|
| 1157 |
* ------------------------------------------------------------------ |
| 1158 |
* m8viol finds how much xn violates the nonlinear constraints. |
| 1159 |
* y(*) is output as the vector of violations. |
| 1160 |
* maxvi points to the biggest violation. |
| 1161 |
* vimax is the biggest violation. |
| 1162 |
* virel is the biggest violation normalized by xnorm. |
| 1163 |
* ------------------------------------------------------------------ |
| 1164 |
|
| 1165 |
common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm |
| 1166 |
common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave |
| 1167 |
|
| 1168 |
intrinsic max |
| 1169 |
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 1170 |
|
| 1171 |
* Set y = - fcon - (linear A)*xn, excluding slacks. |
| 1172 |
|
| 1173 |
do 520 i = 1, nncon |
| 1174 |
y(i) = - fcon(i) |
| 1175 |
520 continue |
| 1176 |
call m2aprd( 8, xn, n, y, nncon, |
| 1177 |
$ ne, nka, a, ha, ka, |
| 1178 |
$ z, nwcore ) |
| 1179 |
|
| 1180 |
* See how much y violates the bounds on the nonlinear slacks. |
| 1181 |
|
| 1182 |
vimax = zero |
| 1183 |
maxvi = 1 |
| 1184 |
|
| 1185 |
do 550 i = 1, nncon |
| 1186 |
j = n + i |
| 1187 |
slack = y(i) |
| 1188 |
viol = max( zero, bl(j) - slack, slack - bu(j) ) |
| 1189 |
if (vimax .lt. viol) then |
| 1190 |
vimax = viol |
| 1191 |
maxvi = i |
| 1192 |
end if |
| 1193 |
550 continue |
| 1194 |
|
| 1195 |
xnorm = dnorm1( nb, xn, 1 ) |
| 1196 |
virel = vimax / (one + xnorm) |
| 1197 |
|
| 1198 |
* end of m8viol |
| 1199 |
end |
Parent Directory
| 
Revision Log