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

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

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