| 1 |
************************************************************************ |
| 2 |
* |
| 3 |
* File mi70nobj fortran. |
| 4 |
* |
| 5 |
* m7bsg m7chkd m7chkg m7chzq m7fixb |
| 6 |
* m7rg m7rgit m7sdir m7sscv |
| 7 |
* |
| 8 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 9 |
|
| 10 |
subroutine m7bsg ( ms, nn, kb, gsub, grd ) |
| 11 |
|
| 12 |
implicit double precision (a-h,o-z) |
| 13 |
integer kb(ms) |
| 14 |
double precision gsub(nn), grd(ms) |
| 15 |
|
| 16 |
* ------------------------------------------------------------------ |
| 17 |
* m7bsg sets grd = basic and superbasic components of gsub. |
| 18 |
* ------------------------------------------------------------------ |
| 19 |
|
| 20 |
common /m3scal/ sclobj,scltol,lscale |
| 21 |
common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj |
| 22 |
|
| 23 |
parameter ( zero = 0.0d+0 ) |
| 24 |
|
| 25 |
do 20 k = 1, ms |
| 26 |
j = kb(k) |
| 27 |
if (j .le. nn) then |
| 28 |
grd(k) = gsub(j) |
| 29 |
else |
| 30 |
grd(k) = zero |
| 31 |
end if |
| 32 |
20 continue |
| 33 |
|
| 34 |
if (iobj .ne. 0) grd(iobj) = - minimz * sclobj |
| 35 |
|
| 36 |
* end of m7bsg |
| 37 |
end |
| 38 |
|
| 39 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 40 |
|
| 41 |
subroutine m7chkd( n, bl, bu, x, dx, d, nfeas ) |
| 42 |
|
| 43 |
implicit double precision (a-h,o-z) |
| 44 |
double precision bl(n), bu(n), x(n), dx, d(n) |
| 45 |
|
| 46 |
* ------------------------------------------------------------------ |
| 47 |
* m7chkd checks that x + dx*d is feasible. |
| 48 |
* It is used by m7chkg and m8chkj for the cheap gradient checks. |
| 49 |
* Original: Just looked at the sign of d for variables on a bound. |
| 50 |
* 13 Mar 1992: dx added as a parameter to make certain that |
| 51 |
* x + dx*d does not lie outside the bounds. |
| 52 |
* d may be altered to achieve this. |
| 53 |
* ------------------------------------------------------------------ |
| 54 |
|
| 55 |
parameter ( zero = 0.0d+0 ) |
| 56 |
|
| 57 |
nfeas = 0 |
| 58 |
do 500 j = 1, n |
| 59 |
xj = x(j) |
| 60 |
b1 = bl(j) |
| 61 |
b2 = bu(j) |
| 62 |
if (b1 .eq. b2 ) d(j) = zero |
| 63 |
|
| 64 |
if (d(j) .ne. zero) then |
| 65 |
|
| 66 |
* x(j) is not fixed, so there is room to move. |
| 67 |
* If xj + dx*dj is beyond one bound, reverse dj |
| 68 |
* and make sure it is not beyond the other. |
| 69 |
* Give up and use set dj = zero if both bounds are too close. |
| 70 |
|
| 71 |
dj = d(j) |
| 72 |
xnew = xj + dx*dj |
| 73 |
|
| 74 |
if (dj .gt. zero) then |
| 75 |
if (xnew .gt. b2) then |
| 76 |
dj = - dj |
| 77 |
xnew = xj + dx*dj |
| 78 |
if (xnew .lt. b1) dj = zero |
| 79 |
end if |
| 80 |
else |
| 81 |
if (xnew .lt. b1) then |
| 82 |
dj = - dj |
| 83 |
xnew = xj + dx*dj |
| 84 |
if (xnew .gt. b2) dj = zero |
| 85 |
end if |
| 86 |
end if |
| 87 |
|
| 88 |
d(j) = dj |
| 89 |
if (dj .ne. zero) nfeas = nfeas + 1 |
| 90 |
end if |
| 91 |
500 continue |
| 92 |
|
| 93 |
* end of m7chkd |
| 94 |
end |
| 95 |
|
| 96 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 97 |
|
| 98 |
subroutine m7chkg( n, bl, bu, g, g2, |
| 99 |
$ x, da, db, z, nwcore ) |
| 100 |
|
| 101 |
implicit double precision (a-h,o-z) |
| 102 |
double precision bl(n), bu(n), da(n), db(n), |
| 103 |
$ g(n), g2(n), x(n), z(nwcore) |
| 104 |
|
| 105 |
* ------------------------------------------------------------------ |
| 106 |
* This routine checks that the gradient of an n-dimensional |
| 107 |
* function has been defined and programmed correctly. |
| 108 |
* |
| 109 |
* First, a cheap heuristic test is performed, as in |
| 110 |
* subroutine chkgrd by the following authors: |
| 111 |
* Philip E. Gill, Walter Murray, Susan M. Picken and Hazel M. Barber |
| 112 |
* D.N.A.C., National Physical Laboratory, England (circa 1975). |
| 113 |
* |
| 114 |
* Next, a more reliable test is performed on each component of the |
| 115 |
* gradient, for indices in the range jverif(1) thru jverif(2). |
| 116 |
* |
| 117 |
* lverif(1) is the verify level, which has the following meaning: |
| 118 |
* |
| 119 |
* -1 do not perform any check. |
| 120 |
* 0 do the cheap test only. |
| 121 |
* 1 or 3 do both cheap and full test on objective gradients. |
| 122 |
* 2 or 3 do both cheap and full test on the jacobian. |
| 123 |
* ------------------------------------------------------------------ |
| 124 |
|
| 125 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 126 |
common /m1file/ iread,iprint,isumm |
| 127 |
common /m3scal/ sclobj,scltol,lscale |
| 128 |
common /m5log1/ idebug,ierr,lprint |
| 129 |
common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 130 |
common /m8veri/ jverif(4),lverif(2) |
| 131 |
|
| 132 |
intrinsic abs, max, min |
| 133 |
|
| 134 |
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 135 |
|
| 136 |
logical cheap, nodd |
| 137 |
character*4 key |
| 138 |
character*4 lbad , lgood |
| 139 |
data lbad/'bad?'/, lgood/'ok '/ |
| 140 |
|
| 141 |
lvl = lverif(1) |
| 142 |
if (lvl .lt. 0) return |
| 143 |
|
| 144 |
j1 = max( jverif(1), 1 ) |
| 145 |
j2 = min( jverif(2), n ) |
| 146 |
cheap = lvl .eq. 0 .or. lvl .eq. 2 .or. j1 .gt. j2 |
| 147 |
lssave = lscale |
| 148 |
lscale = 0 |
| 149 |
|
| 150 |
* Evaluate the function f and gradient g at the base point x. |
| 151 |
|
| 152 |
call m6dmmy( n, g ) |
| 153 |
call m6fobj( 2, n, f, g, x, z, nwcore ) |
| 154 |
if (ierr .ne. 0) go to 900 |
| 155 |
if (knowng(1) .eq. 0) go to 900 |
| 156 |
|
| 157 |
if (iprint .gt. 0) then |
| 158 |
if ( cheap ) then |
| 159 |
write(iprint, 1100) |
| 160 |
else |
| 161 |
write(iprint, 1000) |
| 162 |
end if |
| 163 |
end if |
| 164 |
|
| 165 |
* -------------------------- |
| 166 |
* Cheap test. |
| 167 |
* -------------------------- |
| 168 |
|
| 169 |
* If n is odd, u is set to 1/(n - 1). Otherwise, u = 1/n. |
| 170 |
|
| 171 |
nodd = 2*(n/2) .ne. n |
| 172 |
rn = n |
| 173 |
r = one/rn |
| 174 |
u = r |
| 175 |
if (nodd .and. n .ne. 1) u = one/(rn - one) |
| 176 |
|
| 177 |
* Set arrays da and db to be (almost) orthogonal. |
| 178 |
* We must not perturb x(j) if g(j) is unknown. |
| 179 |
|
| 180 |
call dload ( n, zero, da, 1 ) |
| 181 |
call dload ( n, zero, db, 1 ) |
| 182 |
do 20 j = 1, n |
| 183 |
if (g(j) .eq. gdummy) go to 20 |
| 184 |
db(j) = r |
| 185 |
da(j) = u |
| 186 |
20 continue |
| 187 |
|
| 188 |
do 30 j = 1, n, 2 |
| 189 |
if (da(j) .ne. zero) da(j) = - da(j) - j / rn |
| 190 |
30 continue |
| 191 |
if (nodd) da(n) = zero |
| 192 |
|
| 193 |
* Define a difference interval. |
| 194 |
* Make sure da and db are feasible directions. |
| 195 |
|
| 196 |
dx = difint(1) * (one + dnorm1( n, x, 1 )) |
| 197 |
call m7chkd( n, bl, bu, x, dx, da, nfeas1 ) |
| 198 |
call m7chkd( n, bl, bu, x, dx, db, nfeas2 ) |
| 199 |
|
| 200 |
if (nfeas1 + nfeas2 .eq. 0) then |
| 201 |
if (iprint .gt. 0) write(iprint, 1200) |
| 202 |
else |
| 203 |
|
| 204 |
* Set w1 = g(t)*da and w2 = g(t)*db. |
| 205 |
|
| 206 |
w1 = ddot ( n, da, 1, g, 1 ) |
| 207 |
w2 = ddot ( n, db, 1, g, 1 ) |
| 208 |
|
| 209 |
* Make a forward-difference approximation to the gradient |
| 210 |
* along da and db. |
| 211 |
|
| 212 |
do 40 i = 1, n |
| 213 |
da(i) = x(i) + dx*da(i) |
| 214 |
40 continue |
| 215 |
f1 = f |
| 216 |
if (n .ne. 1) call m6fobj( 0, n, f1, g2, da, z, nwcore ) |
| 217 |
if (ierr .ne. 0) go to 900 |
| 218 |
v1 = (f1 - f)/dx |
| 219 |
|
| 220 |
do 50 i = 1, n |
| 221 |
db(i) = x(i) + dx*db(i) |
| 222 |
50 continue |
| 223 |
call m6fobj( 0, n, f2, g2, db, z, nwcore ) |
| 224 |
if (ierr .ne. 0) go to 900 |
| 225 |
v2 = (f2 - f)/dx |
| 226 |
|
| 227 |
* c1 and c2 are the differences between approximated and |
| 228 |
* programmed gradient projected along da and db respectively. |
| 229 |
|
| 230 |
c1 = v1 - w1 |
| 231 |
c2 = v2 - w2 |
| 232 |
|
| 233 |
* Set an error indicator if c1 or c2 is too large. |
| 234 |
|
| 235 |
ifail = 0 |
| 236 |
if (c1*c1 .ge. dx*(w1*w1 + one) .or. |
| 237 |
$ c2*c2 .ge. dx*(w2*w2 + one)) ifail = 2 |
| 238 |
|
| 239 |
if (ifail .eq. 0) then |
| 240 |
if (iprint .gt. 0) write(iprint, 1400) |
| 241 |
else |
| 242 |
if (iprint .gt. 0) write(iprint, 1500) |
| 243 |
if (isumm .gt. 0) write(isumm , 1500) |
| 244 |
end if |
| 245 |
if (iprint .gt. 0) write(iprint, 1600) w1,w2,v1,v2 |
| 246 |
end if |
| 247 |
|
| 248 |
* ------------------------------------------------------------------ |
| 249 |
* Check each component by differencing along |
| 250 |
* the coordinate directions. |
| 251 |
* Don't bother printing a line if it looks like an exact zero. |
| 252 |
* ------------------------------------------------------------------ |
| 253 |
|
| 254 |
if (cheap) go to 900 |
| 255 |
if (iprint .gt. 0) write(iprint, 2000) |
| 256 |
top = (one + abs( f ))*eps2 |
| 257 |
emax = - one |
| 258 |
jmax = 0 |
| 259 |
nwrong = 0 |
| 260 |
ngood = 0 |
| 261 |
|
| 262 |
do 200 j = j1, j2 |
| 263 |
xj = x(j) |
| 264 |
gj = g(j) |
| 265 |
if (gj .eq. gdummy) go to 200 |
| 266 |
gabs = one + abs( gj ) |
| 267 |
dx = top / gabs |
| 268 |
x(j) = xj + dx |
| 269 |
call m6fobj( 0, n, fforwd, g2, x, z, nwcore ) |
| 270 |
if (ierr .ne. 0) go to 900 |
| 271 |
|
| 272 |
gdiff = (fforwd - f) / dx |
| 273 |
err = abs( gdiff - gj ) / gabs |
| 274 |
|
| 275 |
if (emax .lt. err) then |
| 276 |
emax = err |
| 277 |
jmax = j |
| 278 |
end if |
| 279 |
|
| 280 |
key = lgood |
| 281 |
if (err .gt. eps5) key = lbad |
| 282 |
if (key .eq. lbad) nwrong = nwrong + 1 |
| 283 |
if (key .eq.lgood) ngood = ngood + 1 |
| 284 |
if (abs( gj ) + err .gt. eps0) then |
| 285 |
if (iprint .gt. 0) write(iprint, 2100) j,xj,dx,gj,gdiff,key |
| 286 |
end if |
| 287 |
x(j) = xj |
| 288 |
200 continue |
| 289 |
|
| 290 |
if (iprint .gt. 0) then |
| 291 |
if (nwrong .eq. 0) then |
| 292 |
write(iprint, 2500) ngood ,j1,j2 |
| 293 |
else |
| 294 |
write(iprint, 2600) nwrong,j1,j2 |
| 295 |
end if |
| 296 |
write(iprint, 2700) emax,jmax |
| 297 |
end if |
| 298 |
if (emax .lt. one) go to 900 |
| 299 |
|
| 300 |
* Bad gradients in funobj. |
| 301 |
|
| 302 |
ierr = 7 |
| 303 |
call m1envt( 1 ) |
| 304 |
if (iprint .gt. 0) write(iprint, 3700) |
| 305 |
if (isumm .gt. 0) write(isumm , 3700) |
| 306 |
|
| 307 |
* Exit. |
| 308 |
|
| 309 |
900 lscale = lssave |
| 310 |
return |
| 311 |
|
| 312 |
1000 format(/// ' Verification of objective gradients', |
| 313 |
$ ' returned by subroutine funobj.') |
| 314 |
1100 format(/ ' Cheap test on funobj...') |
| 315 |
1200 format(/ ' XXX Can''t find a feasible step -', |
| 316 |
$ ' objective gradients not verified.') |
| 317 |
1400 format( ' The objective gradients seem to be OK.') |
| 318 |
1500 format( ' XXX The objective gradients seem to be incorrect.') |
| 319 |
1600 format( ' Gradient projected in two directions', 1p, 2e20.11, |
| 320 |
$ / ' Difference approximations ', 2e20.11) |
| 321 |
2000 format(// 6x, 'j', 7x, 'x(j)', 8x, 'dx(j)', |
| 322 |
$ 11x, 'g(j)', 9x, 'Difference approxn' /) |
| 323 |
2100 format(i7, 1p, e16.8, e10.2, 2e18.8, 2x, a4) |
| 324 |
2500 format(/ i7, ' objective gradients out of', i6, ' thru', i6, |
| 325 |
$ ' seem to be OK.') |
| 326 |
2600 format(/ ' XXX There seem to be', i6, |
| 327 |
$ ' incorrect objective gradients in cols', i6, ' thru', i6) |
| 328 |
2700 format(/ ' XXX The largest relative error was', 1p, e12.2, |
| 329 |
$ ' in column', i6 /) |
| 330 |
3700 format(// ' EXIT -- subroutine funobj appears to be', |
| 331 |
$ ' giving incorrect gradients') |
| 332 |
|
| 333 |
* end of m7chkg |
| 334 |
end |
| 335 |
|
| 336 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 337 |
|
| 338 |
subroutine m7chzq( m, nb, ms, ns, Rset, jq, pivot, |
| 339 |
$ ne, nka, a, ha, ka, |
| 340 |
$ kb, bl, bu, x, y, z, nwcore ) |
| 341 |
|
| 342 |
implicit double precision (a-h,o-z) |
| 343 |
logical Rset |
| 344 |
integer*4 ha(ne) |
| 345 |
integer ka(ne), kb(ms) |
| 346 |
double precision a(ne), bl(nb), bu(nb), x(ms), y(ms), z(nwcore) |
| 347 |
|
| 348 |
* ------------------------------------------------------------------ |
| 349 |
* m7chzq selects a superbasic to replace the jp-th basic variable. |
| 350 |
* On entry, y(1:m) contains the jp-th row of B(inverse). |
| 351 |
* On exit, if Rset is true, y(S) = y(m1:ms) holds the vector v |
| 352 |
* needed by m6bswp to update R following a basis change. |
| 353 |
* |
| 354 |
* 29 Nov 1991: qnewtn changed to Rset (should fix a bug). |
| 355 |
* ------------------------------------------------------------------ |
| 356 |
|
| 357 |
common /m1file/ iread,iprint,isumm |
| 358 |
common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm |
| 359 |
|
| 360 |
intrinsic abs, min |
| 361 |
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 362 |
|
| 363 |
* Set y(S) = 0 - S(t)*y. Beware of the minus sign when using y(S). |
| 364 |
|
| 365 |
m1 = m + 1 |
| 366 |
call dload ( ns, zero, y(m1), 1 ) |
| 367 |
call m2aprd( 4, y, m, y(m1), ns, |
| 368 |
$ ne, nka, a, ha, ka, |
| 369 |
$ z, nwcore ) |
| 370 |
jq = m + idamax( ns, y(m1), 1 ) |
| 371 |
pivot = abs( y(jq) ) |
| 372 |
|
| 373 |
* Exit if the pivot is too small. |
| 374 |
|
| 375 |
if (pivot .lt. tolpiv) then |
| 376 |
if (iprint .gt. 0) then |
| 377 |
write(iprint, '(/ a, 1p, e11.1)') |
| 378 |
$ ' XXX m7chzq. Max pivot is too small:', pivot |
| 379 |
end if |
| 380 |
jq = - ms |
| 381 |
else |
| 382 |
|
| 383 |
* Choose one away from its bounds if possible. |
| 384 |
|
| 385 |
tol = 0.1*pivot |
| 386 |
dmax = - one |
| 387 |
|
| 388 |
do 200 k = m1, ms |
| 389 |
if (abs( y(k) ) .ge. tol) then |
| 390 |
j = kb(k) |
| 391 |
xj = x(k) |
| 392 |
d1 = xj - bl(j) |
| 393 |
d2 = bu(j) - xj |
| 394 |
d1 = min( abs( d1 ), abs( d2 ) ) |
| 395 |
if (dmax .le. d1) then |
| 396 |
dmax = d1 |
| 397 |
jq = k |
| 398 |
end if |
| 399 |
end if |
| 400 |
200 continue |
| 401 |
|
| 402 |
pivot = - y(jq) |
| 403 |
|
| 404 |
* Finish computing v, the vector needed to modify R. |
| 405 |
|
| 406 |
if ( Rset ) then |
| 407 |
y(jq) = - (one + pivot) |
| 408 |
call dscal ( ns, (one/pivot), y(m1), 1 ) |
| 409 |
end if |
| 410 |
end if |
| 411 |
|
| 412 |
* end of m7chzq |
| 413 |
end |
| 414 |
|
| 415 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 416 |
|
| 417 |
subroutine m7fixb( m, maxr, ms, n, nb, nr, ns, nx, inform, |
| 418 |
$ ne, nka, a, ha, ka, |
| 419 |
$ hs, kb, bl, bu, bbl, bbu, |
| 420 |
$ r, x, y, y2, z, nwcore ) |
| 421 |
|
| 422 |
implicit double precision (a-h,o-z) |
| 423 |
integer*4 ha(ne), hs(nb) |
| 424 |
integer ka(ne), kb(ms) |
| 425 |
double precision a(ne), bl(nb), bu(nb), bbl(ms), bbu(ms) |
| 426 |
double precision r(nr), x(ms), y(ms), y2(nx), z(nwcore) |
| 427 |
|
| 428 |
* ------------------------------------------------------------------ |
| 429 |
* m7fixb looks to see if the current basis B is ill-conditioned. |
| 430 |
* We assume that the factorization B = L*U has just been computed. |
| 431 |
* If the smallest diagonal of U seems too small, a superbasic is |
| 432 |
* chosen to be swapped with the corresponding column of B. |
| 433 |
* The reduced Hessian R is updated accordingly. |
| 434 |
* |
| 435 |
* On exit, |
| 436 |
* inform = 0 means no error condition. |
| 437 |
* inform = 5 means the LU update failed. |
| 438 |
* |
| 439 |
* 15 Mar 1992: First version, derived from bits of m7chzq, m7rgit. |
| 440 |
* 18 Mar 1992: If we call m7fixb only after m2bfac, we don't need |
| 441 |
* to make it update grd, pi and rg. |
| 442 |
* ------------------------------------------------------------------ |
| 443 |
|
| 444 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 445 |
common /m1file/ iread,iprint,isumm |
| 446 |
common /m2lu4 / parmlu(30),luparm(30) |
| 447 |
common /m2parm/ dparm(30),iparm(30) |
| 448 |
common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2 |
| 449 |
common /m5lp2 / invrq,invitn,invmod |
| 450 |
|
| 451 |
logical Rset, swap |
| 452 |
parameter ( zero = 0.0d+0, one = 1.0d+0 ) |
| 453 |
|
| 454 |
* LUSOL gives umax = the largest element of U, |
| 455 |
* dumin = the smallest diagonal of U, |
| 456 |
* and jumin = the corresponding column of B. |
| 457 |
|
| 458 |
inform = 0 |
| 459 |
umax = parmlu(12) |
| 460 |
dumin = parmlu(14) |
| 461 |
jumin = luparm(19) |
| 462 |
|
| 463 |
* Since modifying B is an expensive operation, |
| 464 |
* we don't do anything if dumin seems reasonably large. |
| 465 |
* The SPECS file gives swaptl = Swap tolerance = eps4 by default. |
| 466 |
* NOTE: If B has rank 0, dumin = umax = 0. |
| 467 |
* If B has rank 1, dumin = umax. |
| 468 |
* Hence we need both relative and absolute tests. |
| 469 |
|
| 470 |
swaptl = dparm(8) |
| 471 |
if (dumin .gt. swaptl * (umax + one)) return |
| 472 |
|
| 473 |
* U looks ill-conditioned. We want to replace the jumin-th column. |
| 474 |
* Use jp = jumin for brevity, and for the final call to m2bsol. |
| 475 |
|
| 476 |
inform = 1 |
| 477 |
jp = jumin |
| 478 |
Rset = r(1) .ne. zero |
| 479 |
m1 = m + 1 |
| 480 |
nz1 = min( ns, maxr ) |
| 481 |
|
| 482 |
* ------------------------------------------------------------------ |
| 483 |
* Solve B(t)*y = e(jp). |
| 484 |
* ------------------------------------------------------------------ |
| 485 |
call dload ( m, zero, y2, 1 ) |
| 486 |
y2(jp) = one |
| 487 |
call m2bsol( 3, m, y2, y, z, nwcore ) |
| 488 |
|
| 489 |
* ------------------------------------------------------------------ |
| 490 |
* Use y to "price" the superbasics. |
| 491 |
* For each column of S, y(S) returns (minus) the pivot element |
| 492 |
* that would arise if that column replaced the jumin-th basic. |
| 493 |
* ------------------------------------------------------------------ |
| 494 |
|
| 495 |
* Set y(S) = 0 - S(t)*y. Beware of the minus sign when using y(S). |
| 496 |
|
| 497 |
call dload ( ns, zero, y(m1), 1 ) |
| 498 |
call m2aprd( 4, y, m, y(m1), ns, |
| 499 |
$ ne, nka, a, ha, ka, |
| 500 |
$ z, nwcore ) |
| 501 |
kq = idamax( ns, y(m1), 1 ) |
| 502 |
jq = m + kq |
| 503 |
pivot = - y(jq) |
| 504 |
|
| 505 |
* Continue if pivot is noticeably larger than dumin. |
| 506 |
* It is not clear that this the best test, but it will do for now. |
| 507 |
|
| 508 |
cond0 = cond |
| 509 |
dtol = 2.0d+0 |
| 510 |
swap = abs( pivot ) .ge. dtol*dumin |
| 511 |
|
| 512 |
if ( swap ) then |
| 513 |
if (Rset .and. kq .le. maxr) then |
| 514 |
|
| 515 |
* Finish computing v, the vector needed to modify R. |
| 516 |
* v is stored in y(S) for use in m6bswp. |
| 517 |
|
| 518 |
y(jq) = - (one + pivot) |
| 519 |
call dscal ( ns, (one/pivot), y(m1), 1 ) |
| 520 |
|
| 521 |
* Modify R to account for the change in basis. |
| 522 |
* y2 is workspace. |
| 523 |
* Get cond(R) before and after to see how we're doing. |
| 524 |
|
| 525 |
call m6rcnd( maxr, nr, ns, r, dmax, dmin, cond0 ) |
| 526 |
call m6bswp( nz1, nr, r, y2, y(m1), kq, |
| 527 |
$ eps0, eps2, modr2 ) |
| 528 |
call m6rcnd( maxr, nr, ns, r, dmax, dmin, cond ) |
| 529 |
|
| 530 |
* We may want to cancel the swap if the new R |
| 531 |
* seems a lot worse than before. |
| 532 |
* At present we don't know enough about it, so let it go. |
| 533 |
|
| 534 |
worse = 1.0d+3 |
| 535 |
if (cond .gt. worse*cond0) then |
| 536 |
*--- swap = .false. |
| 537 |
end if |
| 538 |
|
| 539 |
* If cond is now pretty bad, set R = I. |
| 540 |
* (Alternatively we could set swap = .false. |
| 541 |
* and restore the old R by calling m6bswp.) |
| 542 |
|
| 543 |
if (cond .ge. one/eps1) then |
| 544 |
r(1) = zero |
| 545 |
call m6rset( maxr, nr, ns, r, cond ) |
| 546 |
end if |
| 547 |
end if |
| 548 |
|
| 549 |
* Swap columns jp and jq of (B S). |
| 550 |
|
| 551 |
if ( swap ) then |
| 552 |
jr1 = kb(jp) |
| 553 |
jq1 = kb(jq) |
| 554 |
kb(jp) = jq1 |
| 555 |
kb(jq) = jr1 |
| 556 |
hs(jq1) = 3 |
| 557 |
hs(jr1) = 2 |
| 558 |
|
| 559 |
if (iprint .gt. 0) |
| 560 |
$ write(iprint, 1000) jr1, jq1, dumin, pivot, cond0, cond |
| 561 |
if (isumm .gt. 0) |
| 562 |
$ write(isumm , 1000) jr1, jq1, dumin, pivot, cond0, cond |
| 563 |
|
| 564 |
t = bbl(jp) |
| 565 |
bbl(jp) = bbl(jq) |
| 566 |
bbl(jq) = t |
| 567 |
|
| 568 |
t = bbu(jp) |
| 569 |
bbu(jp) = bbu(jq) |
| 570 |
bbu(jq) = t |
| 571 |
|
| 572 |
t = x(jp) |
| 573 |
x(jp) = x(jq) |
| 574 |
x(jq) = t |
| 575 |
|
| 576 |
* Set y2 for modifying L and U. L*y2 = a(jq1). |
| 577 |
* Then modify L and U, using jp in /m5log3/. |
| 578 |
|
| 579 |
call m2unpk( jq1, m, n, ne, nka, a, ha, ka, y2 ) |
| 580 |
call m2bsol( 1, m, y2, y, z, nwcore ) |
| 581 |
call m2bsol( 4, m, y2, y, z, nwcore ) |
| 582 |
if (invrq .ne. 0) inform = 5 |
| 583 |
end if |
| 584 |
end if |
| 585 |
|
| 586 |
return |
| 587 |
|
| 588 |
1000 format(1p, ' Swap', 2i6, ' Diag, pivot =', 2e9.1, |
| 589 |
$ ' cond(R) =', 2e9.1) |
| 590 |
|
| 591 |
* end of m7fixb |
| 592 |
end |
| 593 |
|
| 594 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 595 |
|
| 596 |
subroutine m7rg ( m, ms, ns, grd, pi, rg, rgnorm, |
| 597 |
$ ne, nka, a, ha, ka, |
| 598 |
$ z, nwcore ) |
| 599 |
|
| 600 |
implicit double precision (a-h,o-z) |
| 601 |
integer*4 ha(ne) |
| 602 |
integer ka(nka) |
| 603 |
double precision a(ne), grd(ms), pi(m), rg(ns), z(nwcore) |
| 604 |
|
| 605 |
* ------------------------------------------------------------------ |
| 606 |
* m7rg calculates the reduced gradient rg = g(S) - S(t)*pi. |
| 607 |
* ------------------------------------------------------------------ |
| 608 |
|
| 609 |
intrinsic abs, min |
| 610 |
external idamax |
| 611 |
|
| 612 |
call dcopy ( ns, grd(m+1), 1, rg, 1 ) |
| 613 |
call m2aprd( 4, pi, m, rg, ns, |
| 614 |
$ ne, nka, a, ha, ka, |
| 615 |
$ z, nwcore ) |
| 616 |
kmax = idamax( ns, rg, 1 ) |
| 617 |
rgnorm = abs( rg(kmax) ) |
| 618 |
|
| 619 |
* end of m7rg |
| 620 |
end |
| 621 |
|
| 622 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 623 |
|
| 624 |
subroutine m7rgit( m, maxr, maxs, mbs, n, nb, incres, |
| 625 |
$ nn, nn0, nr, ns, nx, inform, nxtphs, |
| 626 |
$ ne, nka, a, ha, ka, |
| 627 |
$ hrtype, hs, kb, bl, bu, bbl, bbu, |
| 628 |
$ fsub, gsub, grd, grd2, |
| 629 |
$ pi, r, rg, rg2, x, xn, y, y2, z, nwcore ) |
| 630 |
|
| 631 |
implicit double precision (a-h,o-z) |
| 632 |
logical incres |
| 633 |
integer*4 ha(ne), hrtype(mbs), hs(nb) |
| 634 |
integer ka(nka), kb(mbs) |
| 635 |
double precision a(ne), bl(nb), bu(nb), bbl(mbs), bbu(mbs), |
| 636 |
$ gsub(nn0), grd(mbs), grd2(mbs), |
| 637 |
$ pi(m), r(nr), rg(maxs), rg2(maxs), |
| 638 |
$ x(mbs), xn(nb), y(nx), y2(nx), z(nwcore) |
| 639 |
|
| 640 |
* ------------------------------------------------------------------ |
| 641 |
* m7rgit performs an iteration of the reduced-gradient algorithm. |
| 642 |
* |
| 643 |
* incres in Phase 3 says if the new variable should increase or not. |
| 644 |
* It is used only in linear mode when m5pric is moving free |
| 645 |
* nonbasics toward their closest bound (and djq = zero). |
| 646 |
* |
| 647 |
* Rset says if a useful quasi-newton r exists. |
| 648 |
* It will be false if a feasible point has not yet been |
| 649 |
* found. |
| 650 |
* Note that Rset could be true even if the current itn |
| 651 |
* is infeasible. This allows r to be updated in |
| 652 |
* certain ways during temporary loss of feasibility. |
| 653 |
* |
| 654 |
* qnewtn in this version (Jan 1983) is the same as nonlin. |
| 655 |
* |
| 656 |
* parhes is true if r is not big enough to store a full |
| 657 |
* quasi-newton estimate of the projected hessian. |
| 658 |
* The null space is then z = ( z1 z2 ), and r estimates |
| 659 |
* the Hessian only in the subspace z1. A diagonal estimate |
| 660 |
* is used (also in r) for the subspace z2. |
| 661 |
* |
| 662 |
* fullz is true if the search direction for the superbasics is |
| 663 |
* computed using all of the reduced gradients. Otherwise, |
| 664 |
* the components corresponding to z2 are zero. |
| 665 |
* |
| 666 |
* grd2 is not used in this version (Aug 1986). |
| 667 |
* |
| 668 |
* xx Sep 1987: Anti-degeneracy m5chzr incorporated. |
| 669 |
* 29 Sep 1991: a, ha, ka etc. passed in. Argument list altered. |
| 670 |
* 04 Oct 1991: switch and difint added as arguments to m6srch. |
| 671 |
* 08 Apr 1992: hs(*) now has internal values. We have to be more |
| 672 |
* careful in setting jrstat, the state of a blocking |
| 673 |
* variable. |
| 674 |
* 15 Apr 1992: incres added as input parameter for linear Phase 3 |
| 675 |
* when djq = 0. |
| 676 |
* ------------------------------------------------------------------ |
| 677 |
|
| 678 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 679 |
common /m1file/ iread,iprint,isumm |
| 680 |
common /m2parm/ dparm(30),iparm(30) |
| 681 |
common /m5loc / lpi ,lpi2 ,lw ,lw2 , |
| 682 |
$ lx ,lx2 ,ly ,ly2 , |
| 683 |
$ lgsub ,lgsub2,lgrd ,lgrd2 , |
| 684 |
$ lr ,lrg ,lrg2 ,lxn |
| 685 |
common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj |
| 686 |
common /m5log1/ idebug,ierr,lprint |
| 687 |
common /m5log2/ jq1,jq2,jr1,jr2,lines1,lines2 |
| 688 |
common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2 |
| 689 |
logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed |
| 690 |
common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed |
| 691 |
common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi |
| 692 |
common /m5lp2 / invrq,invitn,invmod |
| 693 |
common /m5prc / nparpr,nmulpr,kprc,newsb |
| 694 |
common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm |
| 695 |
common /m7len / fobj ,fobj2 ,nnobj ,nnobj0 |
| 696 |
logical conv |
| 697 |
common /m7conv/ etash,etarg,lvltol,nfail,conv(4) |
| 698 |
common /m7phes/ rgmin1,rgnrm1,rgnrm2,jz1,jz2,labz,nfullz,mfullz |
| 699 |
common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg |
| 700 |
common /m8len / njac ,nncon ,nncon0,nnjac |
| 701 |
common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2) |
| 702 |
common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave |
| 703 |
* ------------------------------------------------------------------ |
| 704 |
|
| 705 |
intrinsic abs, max, min, sqrt |
| 706 |
|
| 707 |
logical feasbl, infsbl, linear, nonlin, |
| 708 |
$ fullz , parhes, qnewtn, Rset , |
| 709 |
$ debug , fonly , grdcon, grdobj, |
| 710 |
$ hitlow, move , onbnd , switch, unbndd, |
| 711 |
$ unbddf, unbddx, uncon , vertex |
| 712 |
|
| 713 |
parameter ( zero = 0.0d+0, one = 1.0d+0, |
| 714 |
$ point1 = 0.1d+0, point9 = 0.9d+0 ) |
| 715 |
|
| 716 |
character*19 msg(4:8) |
| 717 |
data msg/'max step too small.', |
| 718 |
$ 'step too small. ', |
| 719 |
$ 'no minimizer. ', |
| 720 |
$ 'too many functions.', |
| 721 |
$ 'uphill direction. '/ |
| 722 |
|
| 723 |
inform = 0 |
| 724 |
m1 = m + 1 |
| 725 |
tolz = eps0 |
| 726 |
feasbl = ninf .eq. 0 |
| 727 |
infsbl = .not. feasbl |
| 728 |
linear = infsbl .or. nn .eq. 0 |
| 729 |
nonlin = .not. linear |
| 730 |
qnewtn = nonlin |
| 731 |
Rset = r(1) .ne. zero |
| 732 |
unbndd = .false. |
| 733 |
bplus = point1 * plinfy |
| 734 |
toobig = sqrt(plinfy) |
| 735 |
oldfx = fsub |
| 736 |
if ( qnewtn ) then |
| 737 |
if (.not. Rset) then |
| 738 |
call m6rset( maxr, nr, ns, r, cond ) |
| 739 |
Rset = .true. |
| 740 |
end if |
| 741 |
end if |
| 742 |
|
| 743 |
if (nphs .eq. 3) then |
| 744 |
* --------------------------------------------------------------- |
| 745 |
* Phase 3. First, make sure we dont want to stay in Phase 4. |
| 746 |
* Since price can select nonbasics floating free between their |
| 747 |
* bounds with zero reduced cost, we have to check that dqj is |
| 748 |
* not zero. |
| 749 |
* --------------------------------------------------------------- |
| 750 |
if (linear) then |
| 751 |
if (djq .eq. zero) then |
| 752 |
djq = one |
| 753 |
if (incres) djq = - one |
| 754 |
rg(ns+1) = djq |
| 755 |
end if |
| 756 |
end if |
| 757 |
|
| 758 |
jq2 = jq |
| 759 |
djqmod = abs( djq ) |
| 760 |
ratio = rgnorm / djqmod |
| 761 |
|
| 762 |
if ( linear ) go to 100 |
| 763 |
if (ns .eq. 0 ) go to 100 |
| 764 |
if (nfail .gt. 0 ) go to 100 |
| 765 |
if (ratio .le. point9 ) go to 100 |
| 766 |
if (rgnorm .le. gtol(2)*pinorm) go to 100 |
| 767 |
|
| 768 |
tolrg = point9 * rgnorm |
| 769 |
jq2 = -jq |
| 770 |
nphs = 4 |
| 771 |
go to 200 |
| 772 |
|
| 773 |
* Add the superbasics selected during pricing. |
| 774 |
* The last newsb entries in kb and rg have been set by m5pric. |
| 775 |
|
| 776 |
100 nfullz = 0 |
| 777 |
jqstat = hs(jq) |
| 778 |
rgnorm = max( rgnorm, djqmod ) |
| 779 |
tolrg = etarg * djqmod |
| 780 |
|
| 781 |
* NOTE. The above line sets the level to which rgnorm |
| 782 |
* must be reduced in the current subspace (Phase 4) |
| 783 |
* before we consider moving off another constraint (Phase 3). |
| 784 |
* etarg between (0, 1) is set by the user at his/her own peril. |
| 785 |
* It is the Subspace tolerance in the SPECS file. |
| 786 |
|
| 787 |
do 130 j = 1, newsb |
| 788 |
ns = ns + 1 |
| 789 |
ms = m + ns |
| 790 |
kq = kb(ms) |
| 791 |
hs(kq) = 2 |
| 792 |
grd(ms) = zero |
| 793 |
if (kq .le. nn .and. feasbl) grd(ms) = gsub(kq) |
| 794 |
hrtype(ms) = 0 |
| 795 |
bbl(ms) = bl(kq) |
| 796 |
bbu(ms) = bu(kq) |
| 797 |
x(ms) = xn(kq) |
| 798 |
if (Rset) call m6radd( maxr, nr, ns, r ) |
| 799 |
130 continue |
| 800 |
|
| 801 |
else |
| 802 |
* --------------------------------------------------------------- |
| 803 |
* Phase 4. Exit if rgnorm is already very small. |
| 804 |
* --------------------------------------------------------------- |
| 805 |
if (rgnorm .le. point1 * toldj(3) * pinorm) go to 910 |
| 806 |
end if |
| 807 |
|
| 808 |
|
| 809 |
* ================================================================== |
| 810 |
* Get a search direction ys for the superbasics (in y(m1)...y(ms)) |
| 811 |
* and then a search direction y for the basics. |
| 812 |
* ================================================================== |
| 813 |
200 ms = m + ns |
| 814 |
nssave = ns |
| 815 |
vertex = ns .eq. 1 |
| 816 |
fullz = .true. |
| 817 |
parhes = nonlin .and. ns .gt. maxr |
| 818 |
nz1 = ns |
| 819 |
if (parhes) nz1 = maxr |
| 820 |
nz2 = ns - maxr |
| 821 |
lastr = maxr*(maxr+1)/2 |
| 822 |
|
| 823 |
* Compute the search direction for the superbasics. |
| 824 |
* rg2 is used for the search direction, but it is promptly |
| 825 |
* copied into y(m1) below. |
| 826 |
* If r is being updated, m7sdir saves w such that R(t)*w = rg. |
| 827 |
* w is used later by m6bfgs. |
| 828 |
|
| 829 |
mode = 0 |
| 830 |
if (nonlin) mode = 1 |
| 831 |
call m7sdir( mode, maxr, nr, ns, r, rg, rg2, z(lw), z, nwcore ) |
| 832 |
call dcopy ( ns, rg2, 1, y(m1), 1 ) |
| 833 |
|
| 834 |
if ( parhes ) then |
| 835 |
|
| 836 |
* Partial Hessian. Let all superbasics move for the first |
| 837 |
* mfullz iterations, but then move only the first set. |
| 838 |
|
| 839 |
mfullz = 3 |
| 840 |
nfullz = nfullz + 1 |
| 841 |
if (nfullz .gt. mfullz) then |
| 842 |
fullz = .false. |
| 843 |
labz = 2 |
| 844 |
j = m1 + maxr |
| 845 |
call dload ( nz2, zero, y(j), 1 ) |
| 846 |
end if |
| 847 |
end if |
| 848 |
|
| 849 |
* ================================================================== |
| 850 |
* The search direction for the superbasics, ys, is now in |
| 851 |
* y(m+1), ..., y(m+ns). |
| 852 |
* Find norms of xs and ys. |
| 853 |
* ================================================================== |
| 854 |
xsnrm1 = dasum( nz1, x(m1), 1 ) |
| 855 |
ysnrm1 = dasum( nz1, y(m1), 1 ) |
| 856 |
if (parhes) then |
| 857 |
xsnrm2 = dasum( nz2, x(m1+maxr), 1 ) |
| 858 |
ysnrm2 = dasum( nz2, y(m1+maxr), 1 ) |
| 859 |
else |
| 860 |
xsnrm2 = zero |
| 861 |
ysnrm2 = zero |
| 862 |
end if |
| 863 |
xsnorm = xsnrm1 + xsnrm2 |
| 864 |
ysnorm = ysnrm1 + ysnrm2 |
| 865 |
|
| 866 |
* Compute y2 = - S*ys and prepare to solve B*y = y2 |
| 867 |
* to get y, the search direction for the basics. |
| 868 |
* We first normalize y2 so the LU solver won't ignore |
| 869 |
* too many "small" elements while computing y. |
| 870 |
|
| 871 |
call dscal ( ns, (one / ysnorm), y(m1), 1 ) |
| 872 |
call dload ( m, zero, y2, 1 ) |
| 873 |
call m2aprd( 2, y(m1), ns, y2, m, |
| 874 |
$ ne, nka, a, ha, ka, |
| 875 |
$ z, nwcore ) |
| 876 |
|
| 877 |
* Solve B*y = y2 and then unnormalize all of y. |
| 878 |
|
| 879 |
call m2bsol( 2, m, y2, y, z, nwcore ) |
| 880 |
call dscal ( ms, ysnorm, y, 1 ) |
| 881 |
ynorm = dasum( m, y, 1 ) + ysnorm |
| 882 |
xnorm = dasum( m, x, 1 ) + xsnorm |
| 883 |
|
| 884 |
* ------------------------------------------------------------------ |
| 885 |
* Find the nearest constraint in direction x + theta*y (theta>=0). |
| 886 |
* With the present version of m5chzr (Sep 1987), theta is always |
| 887 |
* positive. |
| 888 |
* Exact is the step that takes x(jp) exactly onto bound. |
| 889 |
* It may be positive or slightly negative. (Not defined if unbndd.) |
| 890 |
* |
| 891 |
* For the linesearch, theta becomes stepmx, the largest |
| 892 |
* step that the linesearch is allowed to take. |
| 893 |
* |
| 894 |
* If exact is positive, or if we are at a vertex, we do a linesearch |
| 895 |
* to find a new step theta in the range 0 lt theta le stepmx. |
| 896 |
* Sometimes, this interval may be too small for the search to be |
| 897 |
* successful. One thing to remember is that increasing the |
| 898 |
* feasibility tolerance has the effect of increasing stepmx, |
| 899 |
* and therefore increasing the chance of a successful search. |
| 900 |
* |
| 901 |
* If exact isn't positive and we are not at a vertex, we change |
| 902 |
* theta and stepmx to zero and don't move. |
| 903 |
* |
| 904 |
* If onbnd is true, theta is a step that reaches a bound exactly. |
| 905 |
* x(jp) reaches the value bound. If the linesearch says to |
| 906 |
* take a constrained step, bound is used to put the new nonbasic |
| 907 |
* variable xn(jr) exactly on its bound. |
| 908 |
* |
| 909 |
* If unbndd is true, theta = stepmx. |
| 910 |
* ------------------------------------------------------------------ |
| 911 |
|
| 912 |
*xxx stepmn = eps1 / (one + ynorm) > 9 mar 1988: these are wrong if |
| 913 |
*xxx stepmx = 1.0d+12 / (one + ynorm) > ynorm is small (john stone). |
| 914 |
stepmx = 1.0d+12 / ynorm |
| 915 |
tolp = tolpiv * ynorm |
| 916 |
|
| 917 |
call m5chzr( ms , stepmx, plinfy, tolp , |
| 918 |
$ hrtype, bbl , bbu , x , y, |
| 919 |
$ hitlow, move , onbnd , unbndd, |
| 920 |
$ jp , bound , exact , theta ) |
| 921 |
|
| 922 |
if (.not. unbndd) then |
| 923 |
pivot = - y(jp) |
| 924 |
stepmx = theta |
| 925 |
jr = kb(jp) |
| 926 |
end if |
| 927 |
|
| 928 |
if (unbndd .or. vertex .or. exact .gt. zero) then |
| 929 |
* =============================================================== |
| 930 |
* A move looks possible. |
| 931 |
* If linear, do a normal update to x and skip the linesearch. |
| 932 |
* =============================================================== |
| 933 |
if ( linear ) then |
| 934 |
if ( unbndd ) go to 955 |
| 935 |
call daxpy ( ms, theta, y, 1, x, 1 ) |
| 936 |
call m5bsx ( 1, ms, nb, kb, x, xn ) |
| 937 |
go to 500 |
| 938 |
end if |
| 939 |
else |
| 940 |
* =============================================================== |
| 941 |
* Zero step. |
| 942 |
* The blocking variable x(jp) is currently on its bound or |
| 943 |
* slightly infeasible. We shall make it nonbasic and not move. |
| 944 |
* The iteration proceeds as a constrained step with theta = zero |
| 945 |
* and again we skip the linesearch. |
| 946 |
* |
| 947 |
* NOTE: If it is Phase 3 and the blocking variable is the one |
| 948 |
* we have just brought in, we would have to go back and try again |
| 949 |
* without it, to avoid "Cycling in Phase 3". |
| 950 |
* However, since r is expanded with a unit vector, |
| 951 |
* cycling of this kind can never happen. Hence, we no longer |
| 952 |
* guard against it. |
| 953 |
* =============================================================== |
| 954 |
theta = zero |
| 955 |
stepmx = zero |
| 956 |
onbnd = .false. |
| 957 |
go to 500 |
| 958 |
end if |
| 959 |
|
| 960 |
* ------------------------------------------------------------------ |
| 961 |
* Perform a linesearch to find a downhill point (x = x + theta*y). |
| 962 |
* switch tells m6srch whether there is an option to switch to |
| 963 |
* central differences to get a better search direction. |
| 964 |
* |
| 965 |
* m6srch returns the following values: |
| 966 |
* |
| 967 |
* inform =-1 (and ierr = 6) if the user wants to stop. |
| 968 |
* inform = 1 if the search is successful and theta < stepmx. |
| 969 |
* = 2 if the search is successful and theta = stepmx. |
| 970 |
* = 3 if a better point was found but too many functions |
| 971 |
* were needed (not sufficient decrease). |
| 972 |
* = 4 if stepmx < tolabs (too small to do a search). |
| 973 |
* = 5 if theta < alfsml (srchq only -- maybe want to switch |
| 974 |
* to central differences to get a better direction). |
| 975 |
* = 6 if the search found that there is no useful step. |
| 976 |
* The interval of uncertainty is less than 2*tolabs. |
| 977 |
* The minimizer is very close to theta = zero |
| 978 |
* or the gradients are not sufficiently accurate. |
| 979 |
* = 7 if there were too many function calls. |
| 980 |
* = 8 if the input parameters were bad |
| 981 |
* (stepmx le toltny or uphill). |
| 982 |
* ------------------------------------------------------------------ |
| 983 |
debug = itn .ge. iparm(2) |
| 984 |
grdcon = nncon .eq. 0 .or. lderiv .ge. 2 |
| 985 |
grdobj = nnobj .eq. 0 .or. lderiv .eq. 1 .or. lderiv .eq. 3 |
| 986 |
fonly = .not. (grdcon .and. grdobj) |
| 987 |
switch = lvldif .eq. 1 |
| 988 |
$ .and. |
| 989 |
$ ((.not. grdobj .and. knowng(1) .lt. nnobj) .or. |
| 990 |
$ (.not. grdcon .and. knowng(2) .lt. njac ) ) |
| 991 |
|
| 992 |
epsrf = dparm(3) |
| 993 |
damp = dparm(6) * (one + xnorm) / ynorm |
| 994 |
theta = min( one, damp ) |
| 995 |
|
| 996 |
call m6srch( ms, ns, n, nb, nn, itn, inform, |
| 997 |
$ debug, fonly, switch, |
| 998 |
$ ne, nka, a, ha, ka, |
| 999 |
$ theta, stepmx, difint(1), eps , epsrf, etash, |
| 1000 |
$ fsub , rgnorm, ynorm, xnorm, |
| 1001 |
$ gsub , grd, y, x, y2, xn, z, nwcore ) |
| 1002 |
|
| 1003 |
if (inform .lt. 0) return |
| 1004 |
if (inform .le. 3) go to 500 |
| 1005 |
|
| 1006 |
if (inform .eq. 4) then |
| 1007 |
* --------------------------------------------------------------- |
| 1008 |
* The linesearch says stepmx is too small. |
| 1009 |
* (Function precision affects when this happens.) |
| 1010 |
* |
| 1011 |
* If m5chzr set move false, stepmx is very small and we probably |
| 1012 |
* should not have attempted the linesearch. Rather than taking a |
| 1013 |
* zero step (which might lead to cycling if vertex is true) we |
| 1014 |
* force a step of stepmx regardless of the effect on the |
| 1015 |
* objective. The step will be constrained. |
| 1016 |
* onbnd has the correct value. |
| 1017 |
* |
| 1018 |
* Otherwise, we treat it as a linesearch failure |
| 1019 |
* and try for a better direction, possibly with more superbasics. |
| 1020 |
* --------------------------------------------------------------- |
| 1021 |
if ( unbndd ) go to 440 |
| 1022 |
if ( move ) go to 440 |
| 1023 |
|
| 1024 |
theta = stepmx |
| 1025 |
modefg = 2 |
| 1026 |
call daxpy ( ms, theta, y, 1, x, 1 ) |
| 1027 |
call m6fun ( 0, modefg, n, nb, ms, fsub, |
| 1028 |
$ ne, nka, a, ha, ka, |
| 1029 |
$ x, xn, z, nwcore ) |
| 1030 |
call m6fun ( 1, modefg, n, nb, ms, fsub, |
| 1031 |
$ ne, nka, a, ha, ka, |
| 1032 |
$ x, xn, z, nwcore ) |
| 1033 |
if (ierr .ne. 0) return |
| 1034 |
|
| 1035 |
call m6grd ( ms, nb, nn, gsub, grd, |
| 1036 |
$ ne, nka, a, ha, ka, |
| 1037 |
$ xn, z, nwcore ) |
| 1038 |
t = stepmx * ynorm |
| 1039 |
if (iprint .gt. 0) write(iprint, 1075) t |
| 1040 |
go to 500 |
| 1041 |
end if |
| 1042 |
|
| 1043 |
* ------------------------------------------------------------------ |
| 1044 |
* See if we should switch to central differences. |
| 1045 |
* ------------------------------------------------------------------ |
| 1046 |
if (inform .eq. 5) then |
| 1047 |
if (switch .and. lvltol .eq. 2) go to 920 |
| 1048 |
end if |
| 1049 |
|
| 1050 |
* ------------------------------------------------------------------ |
| 1051 |
* Trouble -- no function decrease. |
| 1052 |
* Try resetting the Hessian, deleting constraints, |
| 1053 |
* and finally refactorizing the basis. |
| 1054 |
* Give up after 8 consecutive failures. |
| 1055 |
* ------------------------------------------------------------------ |
| 1056 |
440 if (inform .ge. 4) then |
| 1057 |
if (iprint .gt. 0) |
| 1058 |
$ write(iprint, 1050) inform, msg(inform), itn, rgnorm |
| 1059 |
if (isumm .gt. 0) |
| 1060 |
$ write(isumm , 1050) inform, msg(inform), itn, rgnorm |
| 1061 |
end if |
| 1062 |
if (inform .eq. 5 .and. nfail .lt. 1) nfail = 1 |
| 1063 |
inform = 0 |
| 1064 |
|
| 1065 |
* This is the top of the loop through various recovery procedures. |
| 1066 |
|
| 1067 |
450 nfail = nfail + 1 |
| 1068 |
if ( prnt1 ) write(iprint, 1090) nfail |
| 1069 |
|
| 1070 |
if (nfail .eq. 1) then |
| 1071 |
if ( qnewtn .and. cond .ge. one/eps1 ) then |
| 1072 |
call m6rset( maxr, nr, ns, r, cond ) |
| 1073 |
go to 200 |
| 1074 |
end if |
| 1075 |
|
| 1076 |
else if (nfail .eq. 2) then |
| 1077 |
|
| 1078 |
* Switch to central differences if we have not yet done so |
| 1079 |
* and if we're trying to minimize accurately. |
| 1080 |
|
| 1081 |
if (switch .and. lvltol .eq. 2) go to 920 |
| 1082 |
|
| 1083 |
else if (nfail .le. 5) then |
| 1084 |
|
| 1085 |
* Ask for price (up to 3 times) if rgnorm is not too big. |
| 1086 |
|
| 1087 |
if (nfail .ge. 5 .and. nonopt .eq. 0) go to 475 |
| 1088 |
if (ns .ge. maxs ) go to 475 |
| 1089 |
t = 10.0 |
| 1090 |
if (lvltol .eq. 2) t = one |
| 1091 |
if (rgnorm .gt. t * pinorm ) go to 475 |
| 1092 |
go to 910 |
| 1093 |
|
| 1094 |
475 nfail = 5 |
| 1095 |
|
| 1096 |
else if (nfail .eq. 6) then |
| 1097 |
|
| 1098 |
* Switch to central differences even if we're not trying to |
| 1099 |
* minimize accurately. |
| 1100 |
|
| 1101 |
if (switch .and. lvltol .eq. 1) go to 920 |
| 1102 |
|
| 1103 |
else if (nfail .eq. 7) then |
| 1104 |
|
| 1105 |
* Request refactorization of the basis (inform = 5). |
| 1106 |
|
| 1107 |
if (invitn .gt. 0) go to 925 |
| 1108 |
|
| 1109 |
else if (nfail .eq. 8) then |
| 1110 |
|
| 1111 |
* Request a change of basis as a last resort (inform = 6). |
| 1112 |
|
| 1113 |
go to 930 |
| 1114 |
end if |
| 1115 |
|
| 1116 |
if (nfail .lt. 8) go to 450 |
| 1117 |
|
| 1118 |
* Can't think what to do now other than stop. |
| 1119 |
|
| 1120 |
if (ns .eq. maxs) go to 960 |
| 1121 |
ierr = 9 |
| 1122 |
return |
| 1123 |
|
| 1124 |
* ------------------------------------------------------------------ |
| 1125 |
* We got past the linesearch (or didn't have to do one). |
| 1126 |
* See if the step is unbounded, unconstrained, zero, or otherwise. |
| 1127 |
* ------------------------------------------------------------------ |
| 1128 |
500 inform = 0 |
| 1129 |
unbddf = abs( fsub ) .ge. dparm(1) |
| 1130 |
unbddx = theta*ynorm .ge. dparm(2) |
| 1131 |
unbndd = unbddf .or. unbddx |
| 1132 |
uncon = nonlin .and. theta .lt. stepmx |
| 1133 |
if (unbndd) go to 950 |
| 1134 |
|
| 1135 |
* ================================================================== |
| 1136 |
* Get the new reduced gradient rg2. |
| 1137 |
* ================================================================== |
| 1138 |
modr1 = 0 |
| 1139 |
modr2 = 0 |
| 1140 |
if (linear .or. theta .eq. zero) then |
| 1141 |
call dcopy ( ns, rg, 1, rg2, 1 ) |
| 1142 |
else |
| 1143 |
call dcopy ( m, grd, 1, y, 1 ) |
| 1144 |
call m5setp( 1, m, y, pi, z, nwcore ) |
| 1145 |
call m7rg ( m, ms, ns, grd, pi, rg2, rgnorm, |
| 1146 |
$ ne, nka, a, ha, ka, |
| 1147 |
$ z, nwcore ) |
| 1148 |
|
| 1149 |
* Update the reduced Hessian. |
| 1150 |
|
| 1151 |
if (qnewtn) then |
| 1152 |
nsmove = ns |
| 1153 |
if (parhes .and. .not. fullz) nsmove = maxr |
| 1154 |
call m6bfgs( maxr, nsmove, nr, r, rg, rg2, y(m1), z(lw), |
| 1155 |
$ theta, eps2, eps0, modr1 ) |
| 1156 |
end if |
| 1157 |
end if |
| 1158 |
|
| 1159 |
* ------------------------------------------------------------------ |
| 1160 |
* Update the active constraint set if necessary, and modify R, |
| 1161 |
* the Cholesky factorization of the approximate reduced Hessian. |
| 1162 |
* ------------------------------------------------------------------ |
| 1163 |
if ( uncon ) then |
| 1164 |
* =============================================================== |
| 1165 |
* The step is unconstrained. |
| 1166 |
* =============================================================== |
| 1167 |
pivot = zero |
| 1168 |
kmodlu = 0 |
| 1169 |
call dcopy ( ns, rg2, 1, rg, 1 ) |
| 1170 |
|
| 1171 |
else |
| 1172 |
if (jp .eq. 0) go to 955 |
| 1173 |
|
| 1174 |
* =============================================================== |
| 1175 |
* There is a blocking variable. |
| 1176 |
* It could be a fixed variable, whose new state must be 4. |
| 1177 |
* =============================================================== |
| 1178 |
if (bbl(jp) .eq. bbu(jp)) then |
| 1179 |
jrstat = 4 |
| 1180 |
else if (hitlow) then |
| 1181 |
jrstat = 0 |
| 1182 |
else |
| 1183 |
jrstat = 1 |
| 1184 |
end if |
| 1185 |
if ( onbnd ) xn(jr) = bound |
| 1186 |
|
| 1187 |
if (jp .le. m) then |
| 1188 |
* ============================================================ |
| 1189 |
* A variable in B hit a bound. |
| 1190 |
* Find a column in S to replace it. |
| 1191 |
* ============================================================ |
| 1192 |
|
| 1193 |
* Solve B(t)*y = e(jp). |
| 1194 |
|
| 1195 |
call dload ( m, zero, y2, 1 ) |
| 1196 |
y2(jp) = one |
| 1197 |
call m2bsol( 3, m, y2, y, z, nwcore ) |
| 1198 |
|
| 1199 |
* Select a superbasic to become basic. |
| 1200 |
* If Rset is true, y(S) is returned for use in m6bswp. |
| 1201 |
|
| 1202 |
call m7chzq( m, nb, ms, ns, Rset, jq, pivot, |
| 1203 |
$ ne, nka, a, ha, ka, |
| 1204 |
$ kb, bl, bu, x, y, z, nwcore ) |
| 1205 |
|
| 1206 |
if (jq .le. 0) go to 940 |
| 1207 |
hs(jr) = jrstat |
| 1208 |
jr1 = jr |
| 1209 |
jr2 = kb(jq) |
| 1210 |
jq1 = jr2 |
| 1211 |
kb(jp) = jr2 |
| 1212 |
bbl(jp) = bbl(jq) |
| 1213 |
bbu(jp) = bbu(jq) |
| 1214 |
grd(jp) = grd(jq) |
| 1215 |
x(jp) = x(jq) |
| 1216 |
kq = jq - m |
| 1217 |
|
| 1218 |
* Modify R to account for the change in basis. |
| 1219 |
* y2 is workspace. |
| 1220 |
|
| 1221 |
if (Rset .and. ns .gt. 1 .and. kq .le. maxr) then |
| 1222 |
call m6bswp( nz1, nr, r, y2, y(m1), kq, |
| 1223 |
$ eps0, eps2, modr2 ) |
| 1224 |
end if |
| 1225 |
|
| 1226 |
* Modify pi using y where B(t)*y = e(jp). |
| 1227 |
|
| 1228 |
t = rg2(kq) / pivot |
| 1229 |
call daxpy ( m, t, y, 1, pi, 1 ) |
| 1230 |
pinorm = dnorm1( m, pi, 1 ) |
| 1231 |
pinorm = max( pinorm, one ) |
| 1232 |
|
| 1233 |
* Set y2 for modifying L and U. |
| 1234 |
|
| 1235 |
hs(jq1) = 3 |
| 1236 |
call m2unpk( jq1, m, n, ne, nka, a, ha, ka, y2 ) |
| 1237 |
call m2bsol( 1, m, y2, y, z, nwcore ) |
| 1238 |
|
| 1239 |
else |
| 1240 |
* ============================================================ |
| 1241 |
* A variable in S hit a bound. |
| 1242 |
* ============================================================ |
| 1243 |
hs(jr) = jrstat |
| 1244 |
jr2 = jr |
| 1245 |
kmodlu = 0 |
| 1246 |
kq = jp - m |
| 1247 |
end if |
| 1248 |
|
| 1249 |
* =============================================================== |
| 1250 |
* If necessary, swap the largest reduced-gradient in Z2 into |
| 1251 |
* the front of Z2, so it will end up in the end of Z1. |
| 1252 |
* =============================================================== |
| 1253 |
call m7rg ( m, ms, ns, grd, pi, rg, rgnorm, |
| 1254 |
$ ne, nka, a, ha, ka, |
| 1255 |
$ z, nwcore ) |
| 1256 |
|
| 1257 |
rgdel = abs( rg(kq) ) |
| 1258 |
if ( parhes .and. kq .le. maxr ) then |
| 1259 |
call m6swap( m, maxr, nr, ns, ms, kb, bbl, bbu, |
| 1260 |
$ grd, r, rg, x ) |
| 1261 |
end if |
| 1262 |
|
| 1263 |
* Delete the kq-th superbasic, updating R if it exists. |
| 1264 |
|
| 1265 |
call m6rdel( m, maxr, nr, ns, ms, kb, bbl, bbu, |
| 1266 |
$ grd, r, rg, x, kq, Rset ) |
| 1267 |
ns = ns - 1 |
| 1268 |
ms = m + ns |
| 1269 |
nssave = ns |
| 1270 |
nfullz = 0 |
| 1271 |
|
| 1272 |
if (rgnorm .le. rgdel) then |
| 1273 |
rgnorm = zero |
| 1274 |
if (ns .gt. 0 ) then |
| 1275 |
imax = idamax( ns, rg, 1 ) |
| 1276 |
rgnorm = abs( rg(imax) ) |
| 1277 |
end if |
| 1278 |
end if |
| 1279 |
end if |
| 1280 |
|
| 1281 |
* ------------------------------------------------------------------ |
| 1282 |
* Estimate the condition of R(t)*R using the diagonals of R. |
| 1283 |
* ------------------------------------------------------------------ |
| 1284 |
if ( Rset .and. ns .gt. 0 ) then |
| 1285 |
call m6rcnd( maxr, nr, ns, r, dmax, dmin, cond ) |
| 1286 |
|
| 1287 |
if (infsbl) then |
| 1288 |
* Infeasible. If cond is pretty big, throw R away. |
| 1289 |
|
| 1290 |
if (cond .ge. one/eps2) r(1) = zero |
| 1291 |
else |
| 1292 |
* Feasible. |
| 1293 |
* If cond is hoplessly big, try a basis change. |
| 1294 |
* If cond is big but not quite that big, try modifying it. |
| 1295 |
|
| 1296 |
if (cond .ge. one/eps0) go to 930 |
| 1297 |
if (cond .ge. one/eps1) call m6rset( maxr, nr, ns, r, cond ) |
| 1298 |
end if |
| 1299 |
end if |
| 1300 |
|
| 1301 |
* ================================================================== |
| 1302 |
* Test for convergence in the current subspace. |
| 1303 |
* ================================================================== |
| 1304 |
nfail = 0 |
| 1305 |
parhes = nonlin .and. ns .gt. maxr |
| 1306 |
call m7sscv( m, maxr, maxs, ms, nr, ns, parhes, nxtphs, |
| 1307 |
$ fsub, oldfx, theta, xsnorm, xsnrm1, ysnorm, ysnrm1, |
| 1308 |
$ kb, bbl, bbu, grd, r, rg, x ) |
| 1309 |
|
| 1310 |
* If no room for more superbasics, try to stay in phase 4. |
| 1311 |
|
| 1312 |
if (nxtphs .eq. 3 .and. ns .eq. maxs) nxtphs = 4 |
| 1313 |
go to 990 |
| 1314 |
|
| 1315 |
* ------------------------------------------------------------------ |
| 1316 |
* Various exits. |
| 1317 |
* ------------------------------------------------------------------ |
| 1318 |
|
| 1319 |
* rgnorm is small -- return to Phase 3. |
| 1320 |
|
| 1321 |
910 nphs = 3 |
| 1322 |
inform = 1 |
| 1323 |
return |
| 1324 |
|
| 1325 |
* ================================================================== |
| 1326 |
* Switch to central differences. |
| 1327 |
* (Reset ninf to cause m5frmc to set dummy gradients.) |
| 1328 |
* ================================================================== |
| 1329 |
920 lvldif = 2 |
| 1330 |
nfail = 0 |
| 1331 |
ninf = 1 |
| 1332 |
inform = 4 |
| 1333 |
if (prnt1) write(iprint, 1040) |
| 1334 |
return |
| 1335 |
|
| 1336 |
* ================================================================== |
| 1337 |
* Request a basis factorization. |
| 1338 |
* ================================================================== |
| 1339 |
925 invrq = 22 |
| 1340 |
inform = 5 |
| 1341 |
return |
| 1342 |
|
| 1343 |
* ================================================================== |
| 1344 |
* Request a change of basis via m7fixb. |
| 1345 |
* ================================================================== |
| 1346 |
930 invrq = 23 |
| 1347 |
inform = 6 |
| 1348 |
return |
| 1349 |
|
| 1350 |
* ================================================================== |
| 1351 |
* m7chzq failed. |
| 1352 |
* ================================================================== |
| 1353 |
940 nfail = nfail + 1 |
| 1354 |
if (nfail .lt. 5) then |
| 1355 |
|
| 1356 |
* Treat as an unconstrained step. Then ask for price. |
| 1357 |
|
| 1358 |
jr1 = - jr1 |
| 1359 |
kmodlu = 0 |
| 1360 |
kmodpi = 0 |
| 1361 |
nxtphs = 3 |
| 1362 |
go to 990 |
| 1363 |
else |
| 1364 |
|
| 1365 |
* Fail. |
| 1366 |
|
| 1367 |
call m1page( 2 ) |
| 1368 |
if (iprint .gt. 0) write(iprint, 1200) |
| 1369 |
ierr = 11 |
| 1370 |
return |
| 1371 |
end if |
| 1372 |
|
| 1373 |
* ================================================================== |
| 1374 |
* Unbounded. |
| 1375 |
* ================================================================== |
| 1376 |
950 if (iprint .gt. 0) then |
| 1377 |
if (unbddf) write(iprint, 2000) fsub |
| 1378 |
if (unbddx) write(iprint, 2100) theta |
| 1379 |
end if |
| 1380 |
|
| 1381 |
955 if (iprint .gt. 0) write(iprint, 2200) nphs, kb(ms), rg(ns) |
| 1382 |
ierr = 2 |
| 1383 |
return |
| 1384 |
|
| 1385 |
* ================================================================== |
| 1386 |
* Too many superbasics. |
| 1387 |
* ================================================================== |
| 1388 |
960 call m1page( 2 ) |
| 1389 |
if (iprint .gt. 0) write(iprint, 1600) maxs |
| 1390 |
if (isumm .gt. 0) write(isumm , 1600) maxs |
| 1391 |
ierr = 5 |
| 1392 |
return |
| 1393 |
|
| 1394 |
* ------------------------------------------------------------------ |
| 1395 |
* Normal exit. |
| 1396 |
* ------------------------------------------------------------------ |
| 1397 |
990 if (nphs .eq. 4) djq = rgnorm |
| 1398 |
if (.not. fullz) djq = rgnrm1 |
| 1399 |
return |
| 1400 |
|
| 1401 |
1040 format(' Switch to central differences.') |
| 1402 |
1050 format(' Search exit', i3, ' -- ', a, |
| 1403 |
$ ' Itn =', i7, ' Norm rg =', 1p, e11.3) |
| 1404 |
1075 format(' Forcing a small step. alpha * norm(p) =', 1p, e9.1) |
| 1405 |
1090 format(' Retry', i3) |
| 1406 |
1200 format(' EXIT -- cannot find a superbasic to replace', |
| 1407 |
$ ' basic variable') |
| 1408 |
1600 format(' EXIT -- the superbasics limit is too small...', i7) |
| 1409 |
2000 format(/ ' XXX Linesearch has exceeded the Unbounded function', |
| 1410 |
$ ' value. fsub =', 1p, e15.5) |
| 1411 |
2100 format(/ ' XXX Linesearch has exceeded the Unbounded step size.', |
| 1412 |
$ 4x, ' step =', 1p, e13.3) |
| 1413 |
2200 format(/ ' XXX Unbounded in Phase', i3, |
| 1414 |
$ ' Last SB =', i7, 1p, e13.3) |
| 1415 |
|
| 1416 |
* end of m7rgit |
| 1417 |
end |
| 1418 |
|
| 1419 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1420 |
|
| 1421 |
subroutine m7sdir( mode, maxr, nr, ns, r, rg, p, w, z, nwcore ) |
| 1422 |
|
| 1423 |
implicit double precision (a-h,o-z) |
| 1424 |
double precision r(nr), rg(ns), p(ns), w(ns), z(nwcore) |
| 1425 |
|
| 1426 |
* ------------------------------------------------------------------ |
| 1427 |
* m7sdir computes a search direction p for the superbasic |
| 1428 |
* variables, using the current reduced gradient rg. |
| 1429 |
* |
| 1430 |
* mode method |
| 1431 |
* 0 steepest descent: p = - rg |
| 1432 |
* 1 quasi-Newton: R(t)*R p = - rg |
| 1433 |
* ------------------------------------------------------------------ |
| 1434 |
|
| 1435 |
parameter ( one = 1.0d+0 ) |
| 1436 |
|
| 1437 |
if (mode .eq. 0) then |
| 1438 |
|
| 1439 |
* Steepest descent (used when infeasible). |
| 1440 |
|
| 1441 |
do 20 j = 1, ns |
| 1442 |
p(j) = - rg(j) |
| 1443 |
20 continue |
| 1444 |
else |
| 1445 |
|
| 1446 |
* Quasi-Newton. We must save w satisfying R(t)*w = rg. |
| 1447 |
|
| 1448 |
call dcopy ( ns, rg, 1, p, 1 ) |
| 1449 |
call m6rsol( 2, maxr, nr, ns, r, p ) |
| 1450 |
call dcopy ( ns, p, 1, w, 1 ) |
| 1451 |
call m6rsol( 1, maxr, nr, ns, r, p ) |
| 1452 |
call dscal ( ns, (- one), p, 1 ) |
| 1453 |
end if |
| 1454 |
|
| 1455 |
* end of m7sdir |
| 1456 |
end |
| 1457 |
|
| 1458 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1459 |
|
| 1460 |
subroutine m7sscv( m, maxr, maxs, ms, nr, ns, parhes, nxtphs, |
| 1461 |
$ fx, oldfx, theta, xsnorm,xsnrm1, ysnorm,ysnrm1, |
| 1462 |
$ kb, bbl, bbu, grd, r, rg, x ) |
| 1463 |
|
| 1464 |
implicit double precision (a-h,o-z) |
| 1465 |
logical parhes |
| 1466 |
integer kb(ms) |
| 1467 |
double precision bbl(ms), bbu(ms) |
| 1468 |
double precision grd(ms), r(nr), rg(maxs), x(ms) |
| 1469 |
|
| 1470 |
* ------------------------------------------------------------------ |
| 1471 |
* m7sscv (subspace convergence) decides whether or not |
| 1472 |
* optimization should continue on the current subspace. |
| 1473 |
* On exit, nxtphs = 4 means it should, |
| 1474 |
* nxtphs = 3 means it should not. |
| 1475 |
* |
| 1476 |
* parhes is false if this is a linear iteration or if there |
| 1477 |
* is sufficient storage in r for the full projected Hessian. |
| 1478 |
* |
| 1479 |
* parhes is true if r is just a partial Hessian for the |
| 1480 |
* first maxr superbasic variables. The superbasics are then |
| 1481 |
* in two sets Z1 (containing nz1 = maxr variables) |
| 1482 |
* and Z2 (containing nz2 = ns - maxr variables). |
| 1483 |
* |
| 1484 |
* The null-space matrix is similarly partitioned as Z = ( Z1 Z2 ). |
| 1485 |
* |
| 1486 |
* The normal convergence test is first applied to Z1. If it looks |
| 1487 |
* like we are near enough to an optimum in that restricted |
| 1488 |
* subspace, we find the largest reduced gradient in Z2 (index k2). |
| 1489 |
* If this is less than tolrg we exit with nxtphs = 3 (to ask for |
| 1490 |
* price). Otherwise we make room in Z1 for the corresponding |
| 1491 |
* variable by the moving the superbasic with the smallest reduced |
| 1492 |
* gradient in Z1 (index k1) to the end of Z2. |
| 1493 |
* ------------------------------------------------------------------ |
| 1494 |
|
| 1495 |
logical conv |
| 1496 |
common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy |
| 1497 |
common /m7conv/ etash,etarg,lvltol,nfail,conv(4) |
| 1498 |
common /m7phes/ rgmin1,rgnrm1,rgnrm2,jz1,jz2,labz,nfullz,mfullz |
| 1499 |
common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg |
| 1500 |
|
| 1501 |
intrinsic abs |
| 1502 |
parameter ( one = 1.0d+0 ) |
| 1503 |
|
| 1504 |
rgnrm1 = rgnorm |
| 1505 |
if (parhes) then |
| 1506 |
xsnorm = xsnrm1 |
| 1507 |
ysnorm = ysnrm1 |
| 1508 |
k1 = idamax( maxr, rg, 1 ) |
| 1509 |
rgnrm1 = abs( rg(k1) ) |
| 1510 |
end if |
| 1511 |
|
| 1512 |
deltax = theta * ysnorm |
| 1513 |
deltaf = abs( fx - oldfx ) |
| 1514 |
conv(1) = deltax .le. xtol(lvltol) * (one + xsnorm ) |
| 1515 |
conv(2) = deltaf .le. ftol(lvltol) * (one + abs( fx )) |
| 1516 |
conv(3) = rgnrm1 .le. tolrg |
| 1517 |
conv(4) = rgnrm1 .le. 0.1 * tolrg .or. |
| 1518 |
$ rgnrm1 .le. gtol(2) * pinorm |
| 1519 |
|
| 1520 |
nxtphs = 4 |
| 1521 |
if (conv(1).and.conv(2).and.conv(3) .or. conv(4)) nxtphs = 3 |
| 1522 |
if (.not. parhes ) go to 900 |
| 1523 |
if (nxtphs .eq. 4) go to 900 |
| 1524 |
|
| 1525 |
* Swap the largest reduced gradient in Z2 into the front of Z2 |
| 1526 |
* and see if it is significantly large. |
| 1527 |
|
| 1528 |
call m6swap( m, maxr, nr, ns, ms, kb, bbl, bbu, grd, r, rg, x ) |
| 1529 |
k2 = maxr + 1 |
| 1530 |
rgnrm2 = abs( rg(k2) ) |
| 1531 |
if (ns .lt. maxs .and. rgnrm2 .le. tolrg) go to 900 |
| 1532 |
|
| 1533 |
* Find the smallest component of Z1(t)*g. |
| 1534 |
|
| 1535 |
rgmin1 = abs( rg(1) ) |
| 1536 |
k1 = 1 |
| 1537 |
do 200 k = 1, maxr |
| 1538 |
if (rgmin1 .ge. abs( rg(k) )) then |
| 1539 |
rgmin1 = abs( rg(k) ) |
| 1540 |
k1 = k |
| 1541 |
end if |
| 1542 |
200 continue |
| 1543 |
|
| 1544 |
if (rgmin1 .lt. rgnrm2) then |
| 1545 |
|
| 1546 |
* Save the relevant values. |
| 1547 |
|
| 1548 |
nxtphs = 4 |
| 1549 |
nfullz = 0 |
| 1550 |
lastr = maxr*(maxr + 1)/2 |
| 1551 |
ldiag1 = k1 * (k1 + 1)/2 |
| 1552 |
rdiag1 = r(ldiag1) |
| 1553 |
rg1 = rg(k1) |
| 1554 |
k = m + k1 |
| 1555 |
jz1 = kb(k) |
| 1556 |
jz2 = kb(m + k2) |
| 1557 |
bl1 = bbl(k) |
| 1558 |
bu1 = bbu(k) |
| 1559 |
grd1 = grd(k) |
| 1560 |
x1 = x(k) |
| 1561 |
|
| 1562 |
* Delete the k1-th variable from Z1, and shift the remaining |
| 1563 |
* superbasics in Z1 and Z2 one place to the left. |
| 1564 |
|
| 1565 |
call m6rdel( m, maxr, nr, ns, ms, |
| 1566 |
$ kb, bbl, bbu, grd, r, rg, x, k1, .true. ) |
| 1567 |
|
| 1568 |
* Put the old k1-th superbasic in at the very end. |
| 1569 |
|
| 1570 |
lr = lastr + (ns - maxr) |
| 1571 |
r(lr) = rdiag1 |
| 1572 |
rg(ns) = rg1 |
| 1573 |
kb(ms) = jz1 |
| 1574 |
bbl(ms) = bl1 |
| 1575 |
bbu(ms) = bu1 |
| 1576 |
grd(ms) = grd1 |
| 1577 |
x(ms) = x1 |
| 1578 |
end if |
| 1579 |
|
| 1580 |
900 return |
| 1581 |
|
| 1582 |
* end of m7sscv |
| 1583 |
end |
Parent Directory
| 
Revision Log