/[ascend]/trunk/minos54/mi70nobj.f
ViewVC logotype

Contents of /trunk/minos54/mi70nobj.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (19 years, 7 months ago) by aw0a
File size: 55792 byte(s)
Setting up web subdirectory in repository
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

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22