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

Contents of /trunk/minos54/mi50lp.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (17 years, 9 months ago) by aw0a
File size: 104708 byte(s)
Setting up web subdirectory in repository
1 ************************************************************************
2 *
3 * file mi50lp fortran.
4 *
5 * m5bsx m5chzr m5dgen m5frmc m5hs m5log m5lpit
6 * m5pric m5rc m5setp m5setx m5solv
7 *
8 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9
10 subroutine m5bsx ( mode, ms, nb, kb, x, xn )
11
12 implicit double precision (a-h,o-z)
13 integer kb(ms)
14 double precision x(ms), xn(nb)
15
16 * ------------------------------------------------------------------
17 * m5bsx copies basic and superbasic components from x into xn
18 * or vice versa, depending on whether mode = 1 or 2.
19 * ------------------------------------------------------------------
20
21 if (mode .eq. 1) then
22 do 10 k = 1, ms
23 j = kb(k)
24 xn(j) = x(k)
25 10 continue
26 else
27 do 30 k = 1, ms
28 j = kb(k)
29 x(k) = xn(j)
30 30 continue
31 end if
32
33 * end of m5bsx
34 end
35
36 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37
38 subroutine m5chzr( ms , stepmx, plinfy, tolpiv,
39 $ hrtype, bbl , bbu , x , y,
40 $ hitlow, move , onbnd , unbndd,
41 $ jp , bound , exact , alpha )
42
43 implicit double precision (a-h,o-z)
44 integer*4 hrtype(ms)
45 double precision bbl(ms), bbu(ms), x(ms), y(ms)
46 logical hitlow , move , onbnd, unbndd
47
48 * ------------------------------------------------------------------
49 * m5chzr finds a step alpha such that the point x + alpha*y
50 * reaches one of the bounds on x.
51 *
52 * In this version of chuzr, when x is infeasible, the number of
53 * infeasibilities will never increase. If the number stays the
54 * same, the sum of infeasibilities will decrease.
55 * If the number decreases by one or more,
56 * the sum of infeasibilities will usually decrease also, but
57 * occasionally it will increase after the step alpha is taken.
58 * (Convergence is still assured because the number has decreased.)
59 *
60 * Two possible steps are computed as follows:
61 *
62 * alphaf = the maximum step that can be taken without violating
63 * one of the bounds that are currently satisfied.
64 *
65 * alphai = the maximum (nonzero) step that has the property of
66 * reaching a bound that is currently violated,
67 * subject to the pivot being reasonably close to the
68 * maximum pivot among infeasible variables.
69 * (alphai is not defined if x is feasible.)
70 *
71 * alphai is needed occasionally when infeasible, to prevent
72 * going unnecessarily far when alphaf is quite large. It will
73 * always come into effect when x is about to become feasible.
74 * The sum of infeasibilities will decrease initially as alpha
75 * increases from zero, but may start increasing for larger steps.
76 * choosing a large alphai allows several elements of x to
77 * become feasible at the same time.
78 *
79 * In the end, we take alpha = alphaf if x is feasible, or if
80 * alphai > alphap (where alphap is the perturbed step from pass 1).
81 * Otherwise, we take alpha = alphai.
82 *
83 * Input parameters
84 * ----------------
85 * ms is m + 1 for m5lpit, m + ns for m7rgit.
86 * stepmx defines what should be treated as an unbounded step.
87 * plinfy provides insurance for detecting unboundedness.
88 * if alpha reaches a bound as large as plinfy, it is
89 * classed as an unbounded step.
90 * tolpiv is a tolerance to exclude negligible elements of y.
91 * featol (in common) is the current feasibility tolerance used by
92 * m5frmc. typically in the range 0.5*tolx to 0.99*tolx,
93 * where tolx is the featol specified by the user.
94 * tolinc (in common) is used to determine stepmn (see below),
95 * the minimum positive step.
96 * hrtype is set by m5frmc as follows:
97 * hrtype(j) = -2 if x(j) .lt. bl(j) - featol
98 * = 0 if x(j) is feasible
99 * = +2 if x(j) .gt. bu(j) + featol
100 * bbl the lower bounds on the basic and superbasic variables.
101 * bbu the upper bounds on ditto.
102 * x the values of ditto.
103 * y the search direction.
104 *
105 *
106 * Output parameters
107 * -----------------
108 * hitlow = true if a lower bound restricted alpha.
109 * = false otherwise.
110 * move = true if exact ge stepmn (defined at end of code).
111 * onbnd = true if alpha = exact.
112 * this means that the step alpha moves x(jp)
113 * exactly onto one of its bounds, namely bound.
114 * = false if the exact step would be too small
115 * ( exact lt stepmn ).
116 * (with these definitions, move = onbnd).
117 * unbndd = true if alpha = stepmx. jp may possibly be zero.
118 * the parameters hitlow, move, onbnd, bound and exact
119 * should not be used.
120 * jp = the index (if any) such that x(jp) reaches a bound.
121 * bound = the bound value bbl(jp) or bbu(jp) corresponding
122 * to hitlow.
123 * exact = the step that would take x(jp) exactly onto bound.
124 * alpha = an allowable, positive step.
125 * if unbndd is true, alpha = stepmx.
126 * otherwise, alpha = max( stepmn, exact ).
127 *
128 *
129 *
130 *
131 * Original version written November 1981.
132 * The two-pass approach used follows Paula Harris (1973).
133 *
134 * September 1986: Modified to deal with x + alpha*y only, and to
135 * return a few extra parameters.
136 *
137 * 20 aug 1987: EXPAND procedure implemented to deal with
138 * degeneracy in a more rigorous way. The step alphaf is chosen
139 * the same way as Harris, but since this may be negative, this
140 * version insists on returning a positive step, alpha.
141 * Two features make this possible:
142 *
143 * 1. featol increases slightly each iteration.
144 *
145 * 2. The blocking variable, when made nonbasic, is allowed to
146 * retain the value x(jp) + alpha * y(jp), even if this is
147 * not exactly on the blocking bound.
148 *
149 * 15 feb 1988: Modified to prevent a very small pivot being
150 * selected in connection with alphai. For infeasible variables
151 * moving towards their bound, we now require the chosen pivot
152 * to be at least gamma times as large as the biggest available.
153 * This still gives us freedom in pass 2.
154 * gamma = 0.1 and 0.01 seemed to inhibit phase 1 somewhat.
155 * gamma = 0.001 seems to be safe.
156 *
157 * Note: if the weight on the objective is positive, there is still a
158 * danger of small pivots in phase 1. We have to let an "unbounded"
159 * exit occur. m5solv will set wtobj = zero and try again.
160 * ------------------------------------------------------------------
161
162 common /m5step/ featol, tolx0,tolinc,kdegen,ndegen,
163 $ itnfix, nfix(2)
164
165 intrinsic abs, max
166 logical blockf, blocki
167
168 parameter ( zero = 0.0d+00, gamma = 0.001d+00 )
169
170 * ------------------------------------------------------------------
171 * First pass.
172 * For feasible variables, find the step alphap that reaches the
173 * nearest perturbed (expanded) bound. alphap will be slight larger
174 * than the step to the nearest true bound.
175 * For infeasible variables, find the maximum pivot pivmxi.
176 * ------------------------------------------------------------------
177
178 delta = featol
179 alphap = stepmx
180 pivmxi = zero
181
182 do 200 j = 1, ms
183 pivot = y(j)
184 pivabs = abs( pivot )
185 if (pivabs .le. tolpiv) go to 200
186 jtype = hrtype(j)
187 if (pivot .gt. zero ) go to 150
188
189 * x is decreasing.
190 * Test for smaller alphap if lower bound is satisfied.
191
192 if (jtype .lt. 0) go to 200
193 res = x(j) - bbl(j) + delta
194 if (alphap*pivabs .gt. res) alphap = res / pivabs
195
196 * Test for bigger pivot if upper bound is violated.
197
198 if (jtype .gt. 0) pivmxi = max( pivmxi, pivabs )
199 go to 200
200
201 * x is increasing.
202 * Test for smaller alphap if upper bound is satisfied.
203
204 150 if (jtype .gt. 0) go to 200
205 res = bbu(j) - x(j) + delta
206 if (alphap*pivabs .gt. res) alphap = res / pivabs
207
208 * Test for bigger pivot if lower bound is violated.
209
210 if (jtype .lt. 0) pivmxi = max( pivmxi, pivabs )
211 200 continue
212
213 * ------------------------------------------------------------------
214 * Second pass.
215 * For feasible variables, recompute steps without perturbation.
216 * Choose the largest pivot element subject to the step being
217 * no greater than alphap.
218 * For infeasible variables, find the largest step subject to the
219 * pivot element being no smaller than gamma * pivmxi.
220 * ------------------------------------------------------------------
221
222 alphai = zero
223 pivmxf = zero
224 pivmxi = gamma * pivmxi
225 jhitf = 0
226 jhiti = 0
227
228 do 400 j = 1, ms
229 pivot = y(j)
230 pivabs = abs( pivot )
231 if (pivabs .le. tolpiv) go to 400
232 jtype = hrtype(j)
233 if (pivot .gt. zero ) go to 350
234
235 * x is decreasing.
236 * Test for bigger pivot if lower bound is satisfied.
237
238 if (jtype .lt. 0 ) go to 400
239 if (pivabs .le. pivmxf) go to 340
240 res = x(j) - bbl(j)
241 if (alphap*pivabs .lt. res) go to 340
242 pivmxf = pivabs
243 jhitf = j
244
245 * Test for bigger alphai if upper bound is violated.
246
247 340 if (jtype .eq. 0 ) go to 400
248 if (pivabs .lt. pivmxi) go to 400
249 res = x(j) - bbu(j)
250 if (alphai*pivabs .ge. res) go to 400
251 alphai = res / pivabs
252 jhiti = j
253 go to 400
254
255 * x is increasing.
256 * Test for bigger pivot if upper bound is satisfied.
257
258 350 if (jtype .gt. 0 ) go to 400
259 if (pivabs .le. pivmxf) go to 360
260 res = bbu(j) - x(j)
261 if (alphap*pivabs .lt. res) go to 360
262 pivmxf = pivabs
263 jhitf = j
264
265 * Test for bigger alphai if lower bound is violated.
266
267 360 if (jtype .eq. 0 ) go to 400
268 if (pivabs .lt. pivmxi) go to 400
269 res = bbl(j) - x(j)
270 if (alphai*pivabs .ge. res) go to 400
271 alphai = res / pivabs
272 jhiti = j
273 400 continue
274
275 * ------------------------------------------------------------------
276 * See if a feasible and/or infeasible variable blocks.
277 * ------------------------------------------------------------------
278 blockf = jhitf .gt. 0
279 blocki = jhiti .gt. 0
280 unbndd = .not. ( blockf .or. blocki )
281 if ( unbndd ) go to 900
282 if ( .not. blockf ) go to 500
283
284 * ------------------------------------------------------------------
285 * A variable hits a bound for which it is currently feasible.
286 * the corresponding step alphaf is not used, so no need to get it,
287 * but we know that alphaf .le. alphap, the step from pass 1.
288 * ------------------------------------------------------------------
289 jp = jhitf
290 pivot = y(jp)
291 hitlow = pivot .lt. zero
292
293 * If there is a choice between alphaf and alphai, it is probably
294 * best to take alphai (so we can kick the infeasible variable jhiti
295 * out of the basis).
296 * However, we can't if alphai is bigger than alphap.
297
298 if ( .not. blocki ) go to 600
299 if (alphai .gt. alphap) go to 600
300
301 * ------------------------------------------------------------------
302 * An infeasible variable reaches its violated bound.
303 * ------------------------------------------------------------------
304 500 jp = jhiti
305 pivot = y(jp)
306 hitlow = pivot .gt. zero
307
308 * ------------------------------------------------------------------
309 * Try to step exactly onto bound, but make sure the exact step
310 * is sufficiently positive. (exact will be alphaf or alphai.)
311 * Since featol increases by tolinc each iteration, we know that
312 * a step as large as stepmn (below) will not cause any feasible
313 * variables to become infeasible (where feasibility is measured
314 * by the current featol).
315 * ------------------------------------------------------------------
316 600 if ( hitlow ) then
317 bound = bbl(jp)
318 else
319 bound = bbu(jp)
320 end if
321 unbndd = abs( bound ) .ge. plinfy
322 if ( unbndd ) go to 900
323
324 stepmn = tolinc / abs( pivot )
325 exact = (bound - x(jp)) / pivot
326 alpha = max( stepmn, exact )
327 onbnd = alpha .eq. exact
328 move = exact .ge. stepmn
329 if (.not. move) ndegen = ndegen + 1
330 return
331
332 * ------------------------------------------------------------------
333 * Unbounded.
334 * ------------------------------------------------------------------
335 900 alpha = stepmx
336 move = .true.
337 onbnd = .false.
338
339 * end of m5chzr
340 end
341
342 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
343
344 subroutine m5dgen( mode, m, n, nb, ms, inform,
345 $ ne, nka, a, ha, ka,
346 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
347
348 implicit double precision (a-h,o-z)
349 integer*4 ha(ne), hs(nb)
350 integer ka(nka), kb(ms)
351 double precision a(ne), bl(nb), bu(nb)
352 double precision x(ms), xn(nb), y(m), y2(m), z(nwcore)
353
354 * ------------------------------------------------------------------
355 * m5dgen performs most of the manoeuvres associated with degeneracy.
356 * The degeneracy-resolving strategy operates in the following way.
357 *
358 * Over a cycle of iterations, the feasibility tolerance featol
359 * increases slightly (from tolx0 to tolx1 in steps of tolinc).
360 * This ensures that all steps taken will be positive.
361 *
362 * After kdegen consecutive iterations, nonbasic variables within
363 * featol of their bounds are set exactly on their bounds and the
364 * basic variables are recomputed to satisfy ax = b.
365 * featol is then reduced to tolx0 for the next cycle of iterations.
366 *
367 *
368 * If mode = 1, m5dgen initializes the parameters in
369 * common block m5step:
370 *
371 * featol is the current feasibility tolerance.
372 * tolx0 is the minimum feasibility tolerance.
373 * tolx1 is the maximum feasibility tolerance.
374 * tolinc is the increment to featol.
375 * kdegen is the expand frequency (specified by the user).
376 * it is the frequency of resetting featol to tolx0.
377 * ndegen counts the number of degenerate steps (incremented
378 * by m5chzr).
379 * itnfix is the last iteration at which a mode 2 or 3 entry
380 * set nonbasics onto their bound.
381 * nfix(j) counts the number of times a mode 3 entry has
382 * set nonbasics onto their bound,
383 * where j=1 if infeasible, j=2 if feasible.
384 *
385 * tolx0 and tolx1 are both close to the feasibility tolerance tolx
386 * specified by the user. (They must both be less than tolx.)
387 *
388 *
389 * If mode = 2, m5dgen has been called after a cycle of kdegen
390 * iterations. Nonbasic xn(j)s are examined to see if any are
391 * outside their bounds. (It will never be as much as featol).
392 * inform returns how many. Deviations as small as tolz
393 * (e.g. 1.0d-11) are not counted.
394 * If inform is positive, the basic variables are recomputed.
395 * It is assumed that m5solv will then continue iterations.
396 *
397 *
398 * If mode = 3, m5dgen is being called after a subproblem has been
399 * judged optimal, infeasible or unbounded. Nonbasic xn(j)s are
400 * examined as above.
401 *
402 * First version: August 1987.
403 * 22 Sep 1987: itnfix, nfix(j) and maxfix introduced to allow
404 * more than one mode 3 call at the end of a run.
405 * 14 Oct 1991: a, ha, ka passed in to help with minoss.
406 * 08 Apr 1992: Internal values of hs(*) now used.
407 * 04 Jun 1992: mode 2 and 3 fixed to allow for hs(j) = 4.
408 * ------------------------------------------------------------------
409
410 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
411 common /m1file/ iread,iprint,isumm
412 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
413 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
414 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
415 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
416 common /m5step/ featol, tolx0,tolinc,kdegen,ndegen,
417 $ itnfix, nfix(2)
418 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
419 common /m8len / njac ,nncon ,nncon0,nnjac
420 * ------------------------------------------------------------------
421
422 parameter ( zero = 0.0d+0 )
423
424 inform = 0
425 if (mode .eq. 1) then
426
427 * mode = 1.
428 * Initialize at the start of each major iteration.
429 * kdegen is the expand frequency and
430 * tolx is the feasibility tolerance
431 * (specified by the user). They are not changed.
432
433 ndegen = 0
434 itnfix = 0
435 nfix(1) = 0
436 nfix(2) = 0
437 tolx0 = tolx * 0.5
438 tolx1 = tolx * 0.99
439 if (kdegen .lt. 99999999) then
440 tolinc = (tolx1 - tolx0) / kdegen
441 else
442 tolinc = 0.0
443 end if
444 featol = tolx0
445 else
446
447 * mode = 2 or 3.
448 * Initialize local variables maxfix and tolz.
449
450 maxfix = 2
451 tolz = eps1
452 if (mode .eq. 3) then
453
454 * mode = 3.
455 * Return with inform = 0 if the last call was at the
456 * same itn, or if there have already been maxfix calls
457 * with the same state of feasibility.
458
459 if (itnfix .eq. itn ) return
460 if (ninf .gt. 0 ) j = 1
461 if (ninf .eq. 0 ) j = 2
462 if (nfix(j).ge. maxfix) return
463 nfix(j) = nfix(j) + 1
464 end if
465
466 * Set nonbasics on their nearest bound if they are within
467 * the current featol of that bound. hs(j) = 0, 1 or 4.
468 * To avoid go to's we test on diff for all hs(j).
469
470 itnfix = itn
471
472 do 250 j = 1, nb
473 diff = zero
474 js = hs(j)
475
476 if (js .eq. 0) then
477 bnd = bl(j)
478 diff = bnd - xn(j)
479 else if (js .eq. 1 .or. js .eq. 4) then
480 bnd = bu(j)
481 diff = xn(j) - bnd
482 end if
483
484 if (diff .gt. zero) then
485 if (diff .gt. tolz) inform = inform + 1
486 xn(j) = bnd
487 end if
488 250 continue
489
490 * Reset featol to its minimum value.
491
492 featol = tolx0
493 if (inform .gt. 0) then
494
495 * Reset the basic variables.
496 * Set ninf positive to make sure m5frmc tests for feasibility.
497
498 call m5setx( 1, m, n, nb, ms, kb,
499 $ ne, nka, a, ha, ka,
500 $ bl, bu, x, xn, y, y2, z, nwcore )
501 ninf = 1
502 if (prnt1 .or. nncon .eq. 0) then
503 if (iprint .gt. 0) write(iprint, 1000) itn, inform
504 if (isumm .gt. 0) write(isumm , 1000) itn, inform
505 end if
506 end if
507 end if
508
509 return
510
511 1000 format(' Itn', i7, ' --', i7,
512 $ ' nonbasics set on bound, basics recomputed')
513
514 * end of m5dgen
515 end
516
517 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
518
519 subroutine m5frmc( n, nb, nn, ns, ms, maxs,
520 $ lcrash, first, fsub, featol, objadd,
521 $ ne, nka, a, ha, ka,
522 $ bl, bu, bbl, bbu, hrtype, hs, kb,
523 $ grd2, x, xn, z, nwcore )
524
525 implicit double precision (a-h,o-z)
526 logical first
527 integer*4 ha(ne), hrtype(ms), hs(nb)
528 integer ka(nka), kb(ms)
529 double precision a(ne), bl(nb), bu(nb)
530 double precision bbl(ms), bbu(ms),
531 $ grd2(ms), x(ms), xn(nb), z(nwcore)
532
533 * ------------------------------------------------------------------
534 * m5frmc sets up a vector in grd2 to be used to compute pi.
535 * It also defines hrtype to be used in m5chzr.
536 * After m2bfac, or at a newly feasible point, the subproblem
537 * objective function is evaluated.
538 *
539 * lcrash is used at a first feasible point to print a suitable msg.
540 *
541 * 10 Apr 1992: objadd added as an input parameter.
542 * 04 Jun 1992: lcrash added as an input parameter.
543 * ------------------------------------------------------------------
544
545 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
546 common /m1file/ iread,iprint,isumm
547 common /m3scal/ sclobj,scltol,lscale
548 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
549 common /m5log1/ idebug,ierr,lprint
550 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
551 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
552 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
553 common /m5lp2 / invrq,invitn,invmod
554 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
555 common /m8len / njac ,nncon ,nncon0,nnjac
556 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
557
558 logical kludge, feasbl, infsbl,
559 $ linear, nonlin, lincon, nlncon
560
561 parameter ( zero = 0.0d+0, one = 1.0d+0 )
562
563 * ------------------------------------------------------------------
564
565 feasbl = ninf .eq. 0
566 linear = nn .eq. 0
567 if ( first ) feasbl = .false.
568 infsbl = .not. feasbl
569 nonlin = .not. linear
570
571 call dload ( ms, zero, grd2, 1 )
572 if (iobj .ne. 0) grd2(iobj) = - minimz * sclobj
573
574 * Revert to the simplex method if the number of superbasics
575 * has just decreased to zero.
576 * Note that it is ok to do phase 1 simplex on nonlinear problems
577 * as long as we are at a vertex.
578
579 if (ns .eq. 0) then
580 if ( infsbl ) then
581 nphs = 1
582 else if ( linear ) then
583 nphs = 2
584 end if
585 end if
586
587 * Exit if the previous iteration was feasible and the
588 * basis has not just been refactorized.
589
590 if (feasbl .and. invitn .gt. 0) return
591
592 * ------------------------------------------------------------------
593 * Find the current number and sum of infeasibilities.
594 * ------------------------------------------------------------------
595 lincon = nncon .eq. 0
596 nlncon = nncon .gt. 0
597 numinf = 0
598 suminf = zero
599
600 do 100 k = 1, ms
601 hrtype(k) = 0
602 xk = x(k)
603 res = bbl(k) - xk
604 if (res .le. featol) go to 50
605 grd2(k) = - one
606 hrtype(k) = - 2
607 go to 60
608
609 50 res = xk - bbu(k)
610 if (res .le. featol) go to 100
611 grd2(k) = one
612 hrtype(k) = 2
613
614 60 numinf = numinf + 1
615 suminf = suminf + res
616 100 continue
617
618 sinf = suminf
619 ninf = numinf
620
621 if (numinf .gt. 0) then
622 * ---------------------------------------------------------------
623 * Infeasible.
624 * ---------------------------------------------------------------
625
626 * If first iteration, set nphs.
627
628 if (first) then
629 nphs = 1
630 if (ns .gt. 0) nphs = 4
631 end if
632
633 * Set grd2(iobj) to allow for a composite objective.
634
635 if (nphs .eq. 2) nphs = 1
636 kmodpi = 1
637 fsub = sinf
638 fobj = zero
639 if (iobj .ne. 0) grd2(iobj) = - minimz * wtobj * sclobj
640
641 * Print something if the basis has just been refactorized.
642
643 if (prnt1 .and. invitn .eq. 0)
644 $ write(iprint, 1010) itn, numinf, suminf
645
646 else
647 * ---------------------------------------------------------------
648 * Feasible.
649 * ---------------------------------------------------------------
650
651 * Reset nphs if the previous itn was infeasible.
652
653 if (first .or. infsbl) then
654 nphs = 2
655 if ( nonlin ) nphs = 3
656 if (ns .gt. 0) nphs = 4
657 end if
658
659 * Set up the feasible objective value.
660
661 kmodpi = 1
662 flin = zero
663 if (iobj .ne. 0) flin = - x(iobj) * sclobj
664 fsub = minimz * flin
665
666 if (lcrash .eq. 3) then
667
668 * Only the linear equalities (E rows) have been satisfied.
669 * We don't want to evaluate the functions yet.
670 * Print something and quit.
671
672 if (iprint .gt. 0) then
673 if (prnt1 .or. lincon) then
674 write(iprint, 1000)
675 write(iprint, 1015) itn, fsub
676 end if
677 end if
678 if (isumm .gt. 0) then
679 if (summ1 .or. lincon) write(isumm , 1015) itn, fsub
680 end if
681 return
682 end if
683
684 if ( nonlin ) then
685 modefg = 2
686 call m6fun ( 0, modefg, n, nb, ms, fsub,
687 $ ne, nka, a, ha, ka,
688 $ x, xn, z, nwcore )
689 call m6fun ( 1, modefg, n, nb, ms, fsub,
690 $ ne, nka, a, ha, ka,
691 $ x, xn, z, nwcore )
692 if (ierr .ne. 0) return
693 end if
694
695 wtobj = zero
696 obj = minimz * fsub + objadd
697 objtru = flin + fobj + objadd
698
699 * Print something unless Print level = 0 and there are
700 * nonlinear constraints.
701
702 if (iprint .gt. 0) then
703 if (prnt1 .or. (lincon .and. infsbl)) then
704 if (infsbl) write(iprint, 1000)
705 if (lincon) write(iprint, 1020) itn, obj
706 if (nlncon) write(iprint, 1030) itn, objtru, obj
707 if (infsbl) write(iprint, 1000)
708 end if
709 end if
710
711 if (isumm .gt. 0) then
712 if (infsbl) then
713 if (lincon ) write(isumm , 1020) itn, obj
714 if (nlncon .and. summ1) write(isumm , 1040) itn, obj
715 end if
716 end if
717 end if
718
719 return
720
721 1000 format(' ')
722 1010 format(' Itn', i7, ' -- infeasible. Num =', i5,
723 $ ' Sum =', 1p, e18.9)
724 1015 format(' Itn', i7, ' -- linear E rows feasible. Obj =',
725 $ 1p, e18.9)
726 1020 format(' Itn', i7, ' -- feasible solution. Objective =',
727 $ 1p, e18.9)
728 1030 format(' Itn', i7, ' -- feasible subproblem. True obj =',
729 $ 1p, e17.9, 4x, ' Auglag obj =', e17.9)
730 1040 format(' Itn', i7, ' -- feasible subproblem.',
731 $ 1p, 4x, ' Auglag obj =', e17.9)
732
733 * end of m5frmc
734 end
735
736 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
737
738 subroutine m5hs ( mode, nb, bl, bu, hs, xn )
739
740 implicit double precision (a-h,o-z)
741 character*8 mode
742 double precision bl(nb), bu(nb)
743 integer*4 hs(nb)
744 double precision xn(nb)
745
746 * ------------------------------------------------------------------
747 * m5hs is called from m5solv and m8setj.
748 * if mode = 'Internal', m5hs sets hs(j) = -1 or 4 for certain
749 * nonbasic variables. This allows m5pric to operate more
750 * efficiently. The internal values of hs are now as follows:
751 *
752 * hs(j) = -1 Nonbasic between bounds (bl < x < bu )
753 * hs(j) = 0 Nonbasic on lower bound (bl-tol < x <= bl )
754 * hs(j) = 1 Nonbasic on upper bound (bu <= x < bu+tol)
755 * hs(j) = 2 Superbasic
756 * hs(j) = 3 Basic
757 * hs(j) = 4 Nonbasic and fixed (bl = x = bu )
758 *
759 * where 0 <= tol < the feasibility tolerance.
760 *
761 * if mode = 'External', m5hs changes -1 or 4 values to hs(j) = 0,
762 * ready for basis saving and the outside world.
763 *
764 * 08 Apr 1992: First version.
765 * ------------------------------------------------------------------
766
767 if (mode .eq. 'Internal') then
768 * ---------------------------------------------------------------
769 * Change nonbasic hs(j) to internal values (including 4 and -1).
770 * This may change existing internal values if bl and bu have been
771 * changed -- e.g. at the start of each major iteration.
772 * ---------------------------------------------------------------
773 do 100 j = 1, nb
774 if (hs(j) .le. 1) then
775 if (bl(j) .eq. bu(j)) then
776 hs(j) = 4
777 else if (xn(j) .le. bl(j)) then
778 hs(j) = 0
779 else if (xn(j) .ge. bu(j)) then
780 hs(j) = 1
781 else
782 hs(j) = -1
783 end if
784 end if
785 100 continue
786
787 else
788 * ---------------------------------------------------------------
789 * Change hs to external values.
790 * Some nonbasic hs(j) may be 4 or -1. Change them to 0.
791 * ---------------------------------------------------------------
792 do 200 j = 1, nb
793 if (hs(j) .eq. 4 .or. hs(j) .eq. -1) hs(j) = 0
794 200 continue
795 end if
796
797 * end of m5hs
798 end
799
800 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
801
802 subroutine m5log ( m, maxs, mbs, n, nb, nn, ns, nx, fsub, objadd,
803 $ ne, nka, a, ha, ka,
804 $ bl, bu, hs, kb, x, xn, y, z, nwcore )
805
806 implicit double precision (a-h,o-z)
807 double precision a(ne), bl(nb), bu(nb)
808 integer*4 ha(ne), hs(nb)
809 integer ka(nka), kb(mbs)
810 double precision x(nx), xn(nb), y(nx), z(nwcore)
811
812 * ------------------------------------------------------------------
813 * m5log prints the iteration log.
814 * Normally the only parameters used are nn, ns and fsub.
815 * The others are there to allow monitoring of various items for
816 * experimental purposes.
817 * mbs = m + maxs .ge. m + ns, and nx = max( mbs, nn ).
818 * The array y(nx) is available for workspace.
819 *
820 * The print controls are as follows:
821 *
822 * prnt0 is true if print level = 0.
823 * For lp and lc, a brief log is output every k minor iterations,
824 * where k is the log frequency.
825 * For nlc, a brief log is output by m8setj every major iteration.
826 *
827 * prnt1 is true if print level > 0.
828 * A fuller log is output every k minor iterations.
829 *
830 * summ0 and summ1 are the same as prnt0 and prnt1, but are false
831 * if there is no summary file. summary frequency defines k.
832 *
833 * newhed is true if a new heading is required for some reason other
834 * than frequency (e.g., after a basis factorization).
835 *
836 * lines1 and lines2 count how many lines have been output to the
837 * print and summary files respectively since the last heading.
838 * They too may force a new heading.
839 *
840 * 10 Apr 1992: objadd added as an input parameter.
841 * ------------------------------------------------------------------
842
843 logical conv
844 common /m1file/ iread,iprint,isumm
845 common /m2lu3 / lenl,lenu,ncp,lrow,lcol
846 common /m3scal/ sclobj,scltol,lscale
847 common /m5freq/ kchk,kinv,ksav,klog,ksumm,i1freq,i2freq,msoln
848 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
849 common /m5log1/ idebug,ierr,lprint
850 common /m5log2/ jq1,jq2,jr1,jr2,lines1,lines2
851 common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2
852 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
853 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
854 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
855 common /m5lp2 / invrq,invitn,invmod
856 common /m5prc / nparpr,nmulpr,kprc,newsb
857 common /m7conv/ etash,etarg,lvltol,nfail,conv(4)
858 common /m7phes/ rgmin1,rgnrm1,rgnrm2,jz1,jz2,labz,nfullz,mfullz
859 common /m8len / njac ,nncon ,nncon0,nnjac
860 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
861 common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave
862
863 intrinsic mod
864 logical first , newset, phead , pline, hshort, lshort
865 parameter ( zero = 0.0d+0 )
866 character*4 label(2), l
867 data label(1)/' '/, label(2)/'r'/
868 * ------------------------------------------------------------------
869
870 first = minits .le. 1
871 hshort = ns .eq. 0 .and. nn .eq. 0
872 lshort = ns .eq. 0
873 it = mod( itn, 1000000 )
874 ms = m + ns
875 if (ninf .gt. 0) then
876 if (iobj .gt. 0) then
877 obj = - x(iobj) * sclobj
878 else
879 obj = zero
880 end if
881 else
882 obj = fsub
883 end if
884 obj = minimz * obj + objadd
885
886 if (prnt0 .and. nncon .eq. 0) then
887 * ---------------------------------------------------------------
888 * Terse print.
889 * If there are nonlinear constraints, this is handled by m8setj.
890 * ---------------------------------------------------------------
891 lu = lenl + lenu
892 newset = lines1 .ge. 40
893 pline = first .or. mod( itn, klog ) .eq. 0
894 phead = pline .and. ( newhed .or. newset )
895
896 if ( phead ) then
897 newhed = .false.
898 lines1 = 0
899 if ( hshort ) then
900 write(iprint, 1000)
901 else
902 write(iprint, 1100)
903 end if
904 end if
905
906 if ( pline ) then
907 lines1 = lines1 + 1
908 if ( lshort ) then
909 write(iprint, 1200) it, djq, ninf, sinf, obj, lu
910 else
911 write(iprint, 1200) it, djq, ninf, sinf, obj, lu,
912 $ nfobj(1), ns, cond
913 end if
914 end if
915 end if
916
917 if (summ0 .and. nncon .eq. 0) then
918 * ---------------------------------------------------------------
919 * Output to the summary file.
920 * ---------------------------------------------------------------
921 newset = lines2 .ge. 10
922 pline = first .or. mod( itn, ksumm ) .eq. 0
923 phead = pline .and. ( first .or. newset )
924
925 if ( phead ) then
926 lines2 = 0
927 if ( hshort ) then
928 write(isumm , 2000)
929 else
930 write(isumm , 2100)
931 end if
932 end if
933
934 if ( pline ) then
935 lines2 = lines2 + 1
936 if ( lshort ) then
937 write(isumm , 2200) it, djq, ninf, sinf, obj
938 else
939 write(isumm , 2200) it, djq, ninf, sinf, obj, nfobj(1),ns
940 end if
941 end if
942 end if
943
944 if ( prnt1 ) then
945 * ---------------------------------------------------------------
946 * Detailed print.
947 * ---------------------------------------------------------------
948 l = label(labz)
949 newset = lines1 .ge. 40
950 pline = first .or. mod( itn, klog ) .eq. 0
951 phead = pline .and. ( newhed .or. newset )
952
953 if ( phead ) then
954 newhed = .false.
955 lines1 = 0
956 if ( hshort ) then
957 write(iprint, 3000)
958 else
959 write(iprint, 3100)
960 end if
961 end if
962
963 if ( pline ) then
964 lines1 = lines1 + 1
965 if ( ninf .gt. 0 ) then
966 sumobj = sinf
967 else
968 sumobj = obj
969 end if
970
971 if ( lshort ) then
972 write(iprint, 3200) it, l, nphs, kprc, djq, jq2, jr2,jr1,
973 $ theta, pivot, ninf, sumobj, lenl, lenu, ncp
974 else
975 write(iprint, 3200) it, l, nphs, kprc, djq, jq2, jr2,jr1,
976 $ theta, pivot, ninf, sumobj, lenl, lenu, ncp,
977 $ nfobj(1), nfcon(1), ns, modr1, modr2, cond, conv
978 end if
979
980 * Special output if Hessian dimension is less than the current
981 * number of superbasics, and if variable jz2 has just been
982 * moved from set z2 into set z1 in place of variable jz1.
983
984 if ( jz2 .gt. 0 ) then
985 write(iprint, 3300) rgnrm2, jz2, jz1, rgmin1
986 end if
987 end if
988 end if
989
990 if ( summ1 ) then
991 * ---------------------------------------------------------------
992 * Output to the summary file.
993 * ---------------------------------------------------------------
994 newset = lines2 .ge. 10
995 pline = first .or. mod( itn, ksumm ) .eq. 0
996 phead = pline .and. ( first .or. newset )
997
998 if ( phead ) then
999 lines2 = 0
1000 if ( hshort ) then
1001 write(isumm , 4000)
1002 else
1003 write(isumm , 4100)
1004 end if
1005 end if
1006
1007 if ( pline ) then
1008 lines2 = lines2 + 1
1009 if ( lshort ) then
1010 write(isumm , 4200) it, djq, ninf, sinf, obj
1011 else
1012 write(isumm , 4200) it, djq, ninf, sinf, obj,
1013 $ nfobj(1), nfcon(1), ns
1014 end if
1015 end if
1016 end if
1017
1018 * ------------------------------------------------------------------
1019 * Debug output.
1020 * ------------------------------------------------------------------
1021 if (idebug .eq. 100) then
1022 write(iprint, 5000) (kb(k), x(k), k = 1, ms)
1023 end if
1024
1025 * Exit.
1026
1027 return
1028
1029 1000 format(/ ' Itn dj ninf sinf objective',
1030 $ ' LU')
1031 1100 format(/ ' Itn rg ninf sinf objective',
1032 $ ' LU nobj nsb cond(H)')
1033 1200 format(1p, i7, e9.1, i6, e10.3, e16.8, i7, i6, i5, e8.1)
1034 2000 format(/ ' Itn dj ninf sinf objective')
1035 2100 format(/ ' Itn rg ninf sinf objective',
1036 $ ' nobj nsb')
1037 2200 format(1p, i7, e9.1, i6, e10.3, e16.8, i6, i5)
1038 3000 format(/ ' Itn ph pp dj +sbs -sbs -bs',
1039 $ ' step pivot ninf sinf,objective L U ncp')
1040 3100 format(/ ' Itn ph pp rg +sbs -sbs -bs',
1041 $ ' step pivot ninf sinf,objective L U ncp',
1042 $ ' nobj ncon nsb Hmod cond(H) conv')
1043 3200 format(1p, i7, 1x, a1, i1, i3, e9.1, 3i6,
1044 $ e8.1, e9.1, i5, e16.8, 2i6, i4,
1045 $ 2i6, i5, i3, i2, e8.1, 1x, 4l1)
1046 3300 format(1p, 17x, e9.1, 2i5, e9.1)
1047 4000 format(/ ' Itn dj ninf sinf objective')
1048 4100 format(/ ' Itn rg ninf sinf objective',
1049 $ ' nobj ncon nsb')
1050 4200 format(1p, i7, e9.1, i6, e10.3, e16.8, 2i6, i5)
1051 5000 format(/ ' BS and SB values...' / (5(i7, g17.8)))
1052
1053 * end of m5log
1054 end
1055
1056 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1057
1058 subroutine m5lpit( m, m1, n, nb, incres,
1059 $ ne, nka, a, ha, ka,
1060 $ hrtype, hs, kb,
1061 $ bl, bu, bbl, bbu, x, xn, y, y2, z, nwcore )
1062
1063 implicit double precision (a-h,o-z)
1064 logical incres
1065 integer*4 ha(ne), hrtype(m1), hs(nb)
1066 integer ka(nka), kb(m1)
1067 double precision a(ne), bl(nb), bu(nb)
1068 double precision bbl(m1), bbu(m1),
1069 $ x(m1), xn(nb), y(m1), y2(m), z(nwcore)
1070
1071 * ------------------------------------------------------------------
1072 * m5lpit performs an iteration of the primal simplex method.
1073 * jq is the variable entering the basis and djq is its reduced cost.
1074 *
1075 * 00 Aug 1987 This version allows the variable leaving the basis
1076 * to be up to featol outside its bounds.
1077 * 26 Mar 1992 incres now passed in as a parameter from m5pric.
1078 * We used to say incres = djq .lt. zero to mean that
1079 * the new variable x(jq) wants to increase, but that
1080 * didn't allow a rigorous way of eliminating nonbasics
1081 * that are between their bounds with djq = zero.
1082 * 08 Apr 1992 Internal values of hs(*) mean hs(jq) and hs(jr)
1083 * must be set more carefully.
1084 * ------------------------------------------------------------------
1085
1086 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1087 common /m1file/ iread,iprint,isumm
1088 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1089 common /m5log1/ idebug,ierr,lprint
1090 common /m5log2/ jq1,jq2,jr1,jr2,lines1,lines2
1091 common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2
1092 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1093 common /m5lp2 / invrq,invitn,invmod
1094 common /m5step/ featol, tolx0,tolinc,kdegen,ndegen,
1095 $ itnfix, nfix(2)
1096 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1097
1098 logical hitlow, move, onbnd, unbndd
1099 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1100 * ------------------------------------------------------------------
1101
1102 jq2 = jq
1103 jr2 = jq
1104 hrtype(m1) = 0
1105 bbl(m1) = bl(jq)
1106 bbu(m1) = bu(jq)
1107 x(m1) = xn(jq)
1108
1109 * Unpack column jq into y2 and solve B*y = y2.
1110 * y2 will be altered, and is needed later to modify L and U.
1111
1112 call m2unpk( jq, m, n, ne, nka, a, ha, ka, y2 )
1113 call m2bsol( 2, m, y2, y, z, nwcore )
1114
1115 * Select a variable to be dropped from B.
1116 * m5chzr uses the m+1-th element of hrtype, bbl, bbu, x and y.
1117
1118 if (incres) then
1119 call dscal ( m, (- one), y, 1 )
1120 y(m1) = one
1121 else
1122 y(m1) = - one
1123 end if
1124 stepmx = 1.0d+11
1125
1126 call m5chzr( m1 , stepmx, plinfy, tolpiv,
1127 $ hrtype, bbl , bbu , x , y,
1128 $ hitlow, move , onbnd , unbndd,
1129 $ jp , bound , exact , theta )
1130
1131 if (unbndd) go to 800
1132
1133 * See if there is a basis change.
1134
1135 jr = kb(jp)
1136 if (jp .eq. m1) then
1137
1138 * Variable jq reaches its opposite bound.
1139
1140 if (incres) then
1141 hs(jq) = 1
1142 else
1143 hs(jq) = 0
1144 end if
1145 pivot = zero
1146 kmodlu = 0
1147 if (ninf .eq. 0) kmodpi = 0
1148 else
1149
1150 * Variable jq replaces jr, the jp-th variable of B.
1151 * If jr is a fixed variable, its new state is hs(jr) = 4.
1152
1153 jq1 = jq
1154 jr1 = jr
1155 hs(jq) = 3
1156 if (bbl(jp) .eq. bbu(jp)) then
1157 hs(jr) = 4
1158 else if (hitlow) then
1159 hs(jr) = 0
1160 else
1161 hs(jr) = 1
1162 end if
1163 bbl(jp) = bbl(m1)
1164 bbu(jp) = bbu(m1)
1165 pivot = - y(jp)
1166 end if
1167
1168 * Update the basic variables x and copy them into xn.
1169
1170 call daxpy ( m1, theta, y, 1, x, 1 )
1171 call m5bsx ( 1, m1, nb, kb, x, xn )
1172 kb(jp) = jq
1173 x(jp) = x(m1)
1174 if (onbnd) xn(jr) = bound
1175 go to 900
1176
1177 * The solution is apparently unbounded.
1178
1179 800 if (iprint .gt. 0) then
1180 if (incres) then
1181 write(iprint, 1000) jq
1182 else
1183 write(iprint, 1100) jq
1184 end if
1185 end if
1186 ierr = 2
1187
1188 * Exit.
1189
1190 900 return
1191
1192 1000 format(' Variable', i6, ' can increase indefinitely')
1193 1100 format(' Variable', i6, ' can decrease indefinitely')
1194
1195 * end of m5lpit
1196 end
1197
1198 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1199
1200 subroutine m5pric( m, mbs, n, nb, nn, nn0, ns, maxr, maxs, incres,
1201 $ ne, nka, a, ha, ka,
1202 $ hs, kb, bl, bu, gsub,
1203 $ pi, rc, rg )
1204
1205 implicit double precision (a-h,o-z)
1206 logical incres
1207 integer*4 ha(ne), hs(nb)
1208 integer ka(nka), kb(mbs)
1209 double precision a(ne), bl(nb), bu(nb)
1210 double precision gsub(nn0), pi(m), rc(nb), rg(maxs)
1211
1212 * ------------------------------------------------------------------
1213 * m5pric selects a nonbasic variable to enter the basis,
1214 * using the reduced gradients dj = gsub(j) - pi(t)*a(j).
1215 *
1216 * This version does partial pricing on both structurals and slacks.
1217 * Also, multiple price no longer cancels partial price.
1218 * Dynamic tolerances are used if partial price is in effect.
1219 *
1220 * Partial pricing here means sectional pricing, because the
1221 * columns of A and I are both sliced up into nparpr sections
1222 * of equal size. (The last section of each may be a little bigger,
1223 * since nparpr is unlikely to divide evenly into n or m.)
1224 *
1225 * Input gsub = gradient for nonlinear variables
1226 * (for the subproblem objective function).
1227 * pi = pricing vector.
1228 * kprc = the no. of the section where m5pric last found
1229 * a useful dj.
1230 * (kprc = 0 at the start of each major iteration.)
1231 * toldj(1-2) hold the current told if partial pricing, for
1232 * phase 1 and 2 respectively. told is used to
1233 * determine if a dj is significant.
1234 * toldj(3) holds the specified optimality tolerance.
1235 * biggst keeps track of the biggest dj found during the
1236 * last scan of all sections of ( A I ).
1237 *
1238 * Output kprc = the last section scanned.
1239 * nonopt = the no. of useful djs found in that section.
1240 * jq = best column found.
1241 * djq = best dj.
1242 * newsb = no. of new superbasics selected, if nonopt > 0.
1243 * toldj(1-2) save the current told if partial pricing.
1244 * incres says if variable jq should increase or decrease.
1245 * kb(ns+1), ... kb(ns+newsb) = superbasics selected,
1246 * rg(ns+1), ... rg(ns+newsb) = their reduced gradients.
1247 *
1248 * In the code below,
1249 * the next section of A contains npr1 structurals (j1+1 thru k1),
1250 * the next section of I contains npr2 slacks (j2+1 thru k2).
1251 * If nparpr is rather large, either npr1 or npr2 could be zero,
1252 * but not both.
1253 *
1254 * Original version written September 1981.
1255 *
1256 * 00 sep 1986: Modified to allow nonbasic variables to have
1257 * any value xn(j) between bl(j) and bu(j). Variables that are
1258 * strictly between their bounds are allowed to go either way.
1259 *
1260 *x 00 aug 1987: Modified to let nonbasic variables go either way
1261 *x only if they are greater than bl(j) + featol
1262 *x and less than bu(j) - featol,
1263 *x where featol (in common) is the current feasibility tolerance.
1264 *x
1265 *x 01 sep 1987: Modified to ignore nonbasics whose bounds are
1266 *x featol or less apart. this prevents foolish data such as
1267 *x 0 .le. xn(j) .le. 0.0000001
1268 *x from causing lots of bound swaps. m4chek has done its best to
1269 *x keep xn(j) at zero in such cases. The final solution may appear
1270 *x to be non-optimal if xn(j) should really be at the other bound.
1271 *
1272 * 17 sep 1987: Previous two mods undone. The preceding one allows
1273 * nonbasic xn(j)'s to be frozen at either 0 or 0.0000001
1274 * (in that example), and it may then be impossible to make
1275 * certain basic variables feasible. pilotja strikes again!
1276 *
1277 * We now price as follows. Fixed variables are ignored as always.
1278 * Nonbasics that are on or outside a bound are allowed to move only
1279 * one way, towards the opposite bound.
1280 * Nonbasics strictly between their bounds are allowed to move
1281 * either way.
1282 *
1283 * In the latter case, the worry is that rounding error in previous
1284 * steps (xn = xn + theta*y) might have put a blocking variable
1285 * very slightly inside one of its bounds. Letting it move towards
1286 * that bound will give a very small (and essentially wasted) step.
1287 * To guard against this, m5chzr chooses a step that puts the
1288 * blocking variable exactly onto a bound if possible, and now
1289 * returns the logical onbnd and the relevant value bound so that
1290 * m5lpit and m7rgit will know whether to set xn(jr) = bound exactly.
1291 *
1292 * 16 Sep 1988: Partial pricing strategy altered to reduce told more
1293 * systematically. biggst keeps track of the biggest dj encountered
1294 * during a scan from section 1 to section nparpr. In general, told
1295 * is held fixed until the problem is optimal to within that
1296 * tolerance. told is then reduced to reduce * told.
1297 * (This strategy abandoned April 1992.)
1298 *
1299 * 31 Jan 1992: m5rc used to compute reduced costs.
1300 * 26 Mar 1992: incres added as output parameter.
1301 * Floating nonbasics treated more rigorously
1302 * (i.e. nonbasics with djq = 0 or very small
1303 * that we want to move towards a bound).
1304 * 08 Apr 1992: Internal values of hs(*) now used (-1 and 4).
1305 * xn, bl and bu are no longer needed.
1306 * 09 Apr 1992: Initialize toldj(1-2) = plinfy in m5solv, so told
1307 * is always set on the first iteration of Phase 1 or 2.
1308 * 09 Apr 1992: We have had reduce = 0.2 for a while.
1309 * Tried 0.1, 0.25 and 0.5. Decided to go with 0.25.
1310 * It still gives a moderate number of reductions
1311 * to told, and did a bit better on STAIR and ETAMACRO.
1312 * 11 Aug 1992: Went back to 0.2.
1313 * ------------------------------------------------------------------
1314
1315 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1316 common /m1file/ iread,iprint,isumm
1317 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1318 common /m5log1/ idebug,ierr,lprint
1319 common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2
1320 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1321 common /m5prc / nparpr,nmulpr,kprc,newsb
1322 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1323 common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg
1324
1325 intrinsic abs, max, min
1326 logical feasbl, mulprc, normal
1327 parameter ( zero = 0.0d+0, reduce = 0.2d+0 )
1328 * ------------------------------------------------------------------
1329
1330 maxmp = min( nmulpr, maxs - ns )
1331 feasbl = ninf .eq. 0
1332 mulprc = maxmp .gt. 1
1333 normal = .not. mulprc
1334 djmax = - plinfy
1335 djq = zero
1336 jq = 0
1337 jfree = 0
1338 ms = m + ns
1339 ncan = 0
1340 newsb = 1
1341 nonopt = 0
1342 nprc = 0
1343 nparp = nparpr
1344 npr1 = n / nparp
1345 npr2 = m / nparp
1346 if (max( npr1, npr2 ) .le. 0) nparp = 1
1347
1348 * Set the tolerance for a significant dj.
1349
1350 tolmin = toldj(3) * pinorm
1351 lvldj = 1
1352 if ( feasbl ) lvldj = 2
1353 told = toldj(lvldj)
1354 if (nparp .eq. 1) told = tolmin
1355
1356 * Set pointers to the next section of A and I.
1357 * nprc counts how many sections have been scanned in this call.
1358 * kprc keeps track of which one to start with.
1359
1360 100 nprc = nprc + 1
1361 kprc = kprc + 1
1362 if (kprc .gt. nparp) kprc = 1
1363
1364 npr1 = n / nparp
1365 j1 = (kprc - 1)*npr1
1366 k1 = j1 + npr1
1367 if (kprc .eq. nparp) k1 = n
1368 npr1 = max( 0, k1-j1 )
1369
1370 npr2 = m / nparp
1371 j2 = n + (kprc - 1)*npr2
1372 k2 = j2 + npr2
1373 if (kprc .eq. nparp) k2 = nb
1374 npr2 = max( 0, k2-j2 )
1375
1376 * ------------------------------------------------------------------
1377 * Main loops for partial pricing (or full pricing).
1378 * Compute reduced costs rc(*)
1379 * for the kprc-th section of structurals
1380 * and the kprc-th section of slacks.
1381 * ------------------------------------------------------------------
1382 call m5rc ( j1+1, k1, feasbl,
1383 $ m, n, nn, nn0,
1384 $ ne, nka, a, ha, ka,
1385 $ hs, gsub, pi, rc )
1386
1387 do 200 j = j2+1, k2
1388 rc(j) = - pi(j - n)
1389 200 continue
1390
1391 * ------------------------------------------------------------------
1392 * Main loop for testing rc(*).
1393 * dj is rc(j), the reduced cost.
1394 * d is -dj or +dj, depending on which way x(j) can move.
1395 * We are looking for the largest d (which will be positive).
1396 * ------------------------------------------------------------------
1397 np = npr1 + npr2
1398 j = j1
1399 jslk = npr1 + 1
1400
1401 if ( normal ) then
1402 * ---------------------------------------------------------------
1403 * Regular case -- no multiple pricing.
1404 * ---------------------------------------------------------------
1405 do 400 jj = 1, np
1406 if (jj .eq. jslk) j = j2
1407 j = j + 1
1408 js = hs(j)
1409
1410 if (js .le. 1) then
1411 dj = rc(j)
1412
1413 if (js .eq. 0) then
1414 * xj is allowed to increase.
1415 d = - dj
1416 else if (js .eq. 1) then
1417 * xj is allowed to decrease.
1418 d = dj
1419 else
1420 * xj is free to move either way.
1421 * Remember him as jfree in case he is the only one.
1422 d = abs( dj )
1423 jfree = j
1424 end if
1425
1426 * See if this dj is significant.
1427 * Also see if it is the biggest dj so far.
1428
1429 if (d .gt. told) nonopt = nonopt + 1
1430 if (djmax .lt. d) then
1431 djmax = d
1432 djq = dj
1433 jq = j
1434 kpsav = kprc
1435 end if
1436 end if
1437 400 continue
1438
1439 else
1440 * ---------------------------------------------------------------
1441 * Multiple pricing. The first part of the loop is the same.
1442 * ---------------------------------------------------------------
1443 do 500 jj = 1, np
1444 if (jj .eq. jslk) j = j2
1445 j = j + 1
1446 js = hs(j)
1447
1448 if (js .le. 1) then
1449 dj = rc(j)
1450
1451 if (js .eq. 0) then
1452 * xj is allowed to increase.
1453 d = - dj
1454 else if (js .eq. 1) then
1455 * xj is allowed to decrease.
1456 d = dj
1457 else
1458 * xj is free to move either way.
1459 * Remember him as jfree in case he is the only one.
1460 d = abs( dj )
1461 jfree = j
1462 end if
1463
1464 * See if this dj is significant.
1465 * Also see if it is the biggest dj so far.
1466
1467 if (d .gt. told) nonopt = nonopt + 1
1468 if (djmax .lt. d) then
1469 djmax = d
1470 djq = dj
1471 jq = j
1472 kpsav = kprc
1473 end if
1474
1475 * ========================================
1476 * Extra stuff for multiple pricing.
1477 * We have to build the list of candidates.
1478 * ========================================
1479 if (d .gt. tolmin) then
1480
1481 * Search the list backwards,
1482 * starting with smallest rg values.
1483
1484 do 410 i = ncan, 1, -1
1485 msi = ms + i
1486 nsi = ns + i
1487
1488 * See if this list element will stay put.
1489 * If not, make sure we have space to move it.
1490
1491 if (d .le. abs( rg(nsi) )) go to 420
1492
1493 if (i .lt. maxmp) then
1494 * This element lost the contest but is
1495 * staying in the list.
1496 kb(msi + 1) = kb(msi)
1497 rg(nsi + 1) = rg(nsi)
1498 end if
1499 410 continue
1500
1501 * This is the biggest dj so far.
1502
1503 i = 0
1504
1505 * If room, put the new candidate into the list.
1506
1507 420 if (i .lt. maxmp) then
1508 i = i + 1
1509 kb(ms + i) = j
1510 rg(ns + i) = dj
1511 ncan = min( ncan + 1, maxmp )
1512 end if
1513 end if
1514 end if
1515 500 continue
1516 end if
1517 * ------------------------------------------------------------------
1518 * End of loop looking for biggest dj in the kprc-th section.
1519 * ------------------------------------------------------------------
1520
1521 if (nonopt .eq. 0) then
1522 if (nparp .gt. 1) then
1523 * ============================================================
1524 * No significant dj has been found. (All are less than told.)
1525 * Price the next section, if any remain.
1526 * ============================================================
1527 if (nprc .lt. nparp) go to 100
1528
1529 * ============================================================
1530 * All sections have been scanned. Reduce told
1531 * and grab the best dj if it is bigger than tolmin.
1532 * ============================================================
1533 if (djmax .gt. tolmin) then
1534 nonopt = 1
1535 kprc = kpsav
1536 told = max( reduce * djmax, tolmin )
1537 toldj(lvldj) = told
1538 if (itn .gt. 0) wtobj = reduce * wtobj
1539 if (idebug .ge. 1) then
1540 if (iprint .gt. 0)
1541 $ write(iprint, 1000) itn, told, pinorm, wtobj
1542 if (isumm .gt. 0)
1543 $ write(isumm , 1000) itn, told, pinorm, wtobj
1544 end if
1545 end if
1546 end if
1547 end if
1548
1549 * Finish up multiple pricing.
1550 * The new dj's are now in descending order, with the biggest
1551 * in rg(ns+1). Truncate the list if any are much smaller.
1552
1553 if ( mulprc ) then
1554 d = reduce * djmax
1555 do 550 j = 2, ncan
1556 if (abs( rg(ns+j) ) .ge. d) newsb = newsb + 1
1557 550 continue
1558 end if
1559
1560 * -----------------------------------------------------------------
1561 * Finish if we found someone nonoptimal (nonopt .gt. 0)
1562 * or if there's a nonbasic floating free
1563 * between its bounds (jfree .eq. 0)
1564 * or if the problem is nonlinear (nn .gt. 0)
1565 * or there are some superbasics (ns .gt. 0).
1566 * -----------------------------------------------------------------
1567 incres = djq .lt. zero
1568 if (nonopt .gt. 0) go to 900
1569 if (jfree .eq. 0) go to 900
1570 if ( feasbl ) then
1571 if (nn .gt. 0) go to 900
1572 if (ns .gt. 0) go to 900
1573 end if
1574
1575 * -----------------------------------------------------------------
1576 * jfree > 0.
1577 * We prefer not to end an LP problem (or an infeasible problem)
1578 * with some nonbasic variables floating free between their bounds
1579 * (hs(j) = -1). Their dj's will be essentially zero
1580 * but we might as well jam them into the basis.
1581 *
1582 * We leave true free variables alone -- they will probably be zero.
1583 * (To be rigorous, we should test if they ARE zero and move them
1584 * towards zero if not. This has yet to be done...)
1585 * -----------------------------------------------------------------
1586 do 600 j = 1, nb
1587 if (hs(j) .eq. -1) then
1588 b1 = bl(j)
1589 b2 = bu(j)
1590 if (b1 .gt. -plinfy .or. b2 .lt. plinfy) then
1591
1592 * We just found a floating variable with finite bounds.
1593 * Ask for a move towards the bound nearest zero.
1594
1595 incres = abs(b1) .ge. abs(b2)
1596 nonopt = 1
1597 jq = j
1598 djq = rc(j)
1599 go to 900
1600 end if
1601 end if
1602 600 continue
1603
1604 * Exit.
1605
1606 900 kb(ms + 1) = jq
1607 rg(ns + 1) = djq
1608 return
1609
1610 1000 format(' Itn', i7, ' -- toldj =', 1p, e8.1,
1611 $ ' Norm pi =', e8.1, ' wtobj = ', e8.1)
1612
1613 * end of m5pric
1614 end
1615
1616 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1617
1618 subroutine m5rc ( j1, j2, feasbl,
1619 $ m, n, nn, nn0,
1620 $ ne, nka, a, ha, ka,
1621 $ hs, gsub, pi, rc )
1622
1623 implicit double precision (a-h,o-z)
1624 logical feasbl
1625 integer*4 ha(ne), hs(n)
1626 integer ka(nka)
1627 double precision a(ne)
1628 double precision gsub(nn0), pi(m), rc(n)
1629
1630 * ------------------------------------------------------------------
1631 * m5rc computes reduced costs rc(j) for nonbasic columns of A
1632 * in the range j = j1 to j2. It is called by m5pric.
1633 *
1634 * The loop for computing dj for each column could conceivably be
1635 * optimized on some machines. However, there are seldom more than
1636 * 5 or 10 entries in a column.
1637 *
1638 * Note that we could skip fixed variables by passing in the bounds
1639 * and testing if bl(j) .eq. bu(j), but these are relatively rare.
1640 * But see comment for 08 Apr 1992.
1641 *
1642 * 31 Jan 1992: First version
1643 * 08 Apr 1992: Internal values of hs(j) are now used, so fixed
1644 * variables (hs(j) = 4) are skipped as we would like.
1645 * ------------------------------------------------------------------
1646
1647 parameter ( zero = 0.0d+0 )
1648
1649 do 300 j = j1, j2
1650 if (hs(j) .le. 1) then
1651 dj = zero
1652
1653 do 200 i = ka(j), ka(j+1) - 1
1654 ir = ha(i)
1655 dj = dj + pi(ir) * a(i)
1656 200 continue
1657
1658 rc(j) = - dj
1659 end if
1660 300 continue
1661
1662 * Include the nonlinear objective gradient if relevant.
1663
1664 if (feasbl .and. j1 .le. nn) then
1665 do 350 j = j1, min( j2, nn )
1666 if (hs(j) .le. 1) rc(j) = rc(j) + gsub(j)
1667 350 continue
1668 end if
1669
1670 * end of m5rc
1671 end
1672
1673 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1674
1675 subroutine m5setp( mode, m, rhs, pi, z, nwcore )
1676
1677 implicit double precision (a-h,o-z)
1678 double precision pi(m), rhs(m), z(nwcore)
1679
1680 * ------------------------------------------------------------------
1681 * m5setp solves for pi or finds pinorm (or does both).
1682 * mode = 1 means solve B'pi = rhs and get pinorm.
1683 * mode = 2 means solve B'pi = rhs.
1684 * mode = 3 means just get pinorm.
1685 *
1686 * Beware -- in mode 1 and 2, rhs is altered by m2bsol.
1687 *
1688 * 10 Mar 1992: Until now we have required pinorm ge 1.
1689 * However, many users have very small objectives,
1690 * and the effect is to cause premature termination.
1691 * Using pimin < 1 should help fix this.
1692 * 16 Mar 1992: mode used more fully as described above.
1693 * ------------------------------------------------------------------
1694
1695 common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg
1696
1697 intrinsic max
1698 parameter ( pimin = 1.0d-2 )
1699
1700 if (mode .eq. 1 .or. mode .eq. 2) then
1701 call m2bsol( 3, m, rhs, pi, z, nwcore )
1702 end if
1703
1704 if (mode .eq. 1 .or. mode .eq. 3) then
1705 pinorm = dnorm1( m, pi, 1 )
1706 pinorm = max( pinorm, pimin )
1707 end if
1708
1709 * end of m5setp
1710 end
1711
1712 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1713
1714 subroutine m5setx( mode, m, n, nb, ms, kb,
1715 $ ne, nka, a, ha, ka,
1716 $ bl, bu, x, xn, y, y2, z, nwcore )
1717
1718 implicit double precision (a-h,o-z)
1719 integer*4 ha(ne)
1720 integer ka(nka), kb(ms)
1721 double precision a(ne), bl(nb), bu(nb)
1722 double precision x(ms), xn(nb), y(m), y2(m), z(nwcore)
1723
1724 * ------------------------------------------------------------------
1725 * m5setx performs the following functions:
1726 *
1727 * If mode = 0, the slacks are set to satisfy Ax + s = rhs
1728 * if possible. (Thus s = rhs - Ax, but s may
1729 * have to be moved inside its bounds.)
1730 * Called (perhaps) from the top of m5solv.
1731 *
1732 * If mode = 1, the basic components of x are computed to satisfy
1733 * Ax + s = rhs; that is (A I)*xn = rhs.
1734 * Then a row check is performed to see how well
1735 * (A I)*xn = rhs is satisfied.
1736 * y is set to be the row residuals, y = rhs - Ax - s,
1737 * and the row error is norm(y).
1738 * Called from m2bfac, m5dgen, m5solv.
1739 *
1740 * If mode = 2 or more, a row check is performed.
1741 * Called from m5solv.
1742 *
1743 * 14 Oct 1991: a, ha, ka added as parameters to help with minoss.
1744 * 06 Dec 1991: mode 0 added. n, bl, bu added as parameters.
1745 * ------------------------------------------------------------------
1746
1747 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1748 common /m1file/ iread,iprint,isumm
1749 common /m5log1/ idebug,ierr,lprint
1750 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1751 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1752 common /m8len / njac ,nncon ,nncon0,nnjac
1753 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
1754 $ lblslk,lbuslk,lxlam ,lrhs ,
1755 $ lgcon ,lgcon2,lxdif ,lxold
1756
1757 intrinsic abs
1758 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1759
1760 * ------------------------------------------------------------------
1761 * Compute row residuals y = rhs - A*x.
1762 * ------------------------------------------------------------------
1763 if (nncon .gt. 0) call dcopy ( nncon, z(lrhs), 1, y, 1 )
1764 if (nncon .lt. m) call dload ( m-nncon, zero, y(nncon+1), 1 )
1765 call m2aprd( 5, xn, n, y, m,
1766 $ ne, nka, a, ha, ka,
1767 $ z, nwcore )
1768
1769 if (mode .eq. 0) then
1770 * ---------------------------------------------------------------
1771 * Set the slacks.
1772 * s = y in principle, but we have to make s feasible.
1773 * Any elements CLOSE to a bound are moved exactly onto the bound.
1774 * ---------------------------------------------------------------
1775 tolb = eps1
1776 do 10 i = 1, m
1777 j = n + i
1778 b1 = bl(j)
1779 b2 = bu(j)
1780 s = max( y(i), b1 )
1781 s = min( s , b2 )
1782 if ((s - b1) .gt. (b2 - s)) b1 = b2
1783 if (abs(s - b1) .le. tolb) s = b1
1784 xn(j) = s
1785 10 continue
1786
1787 else
1788 * ---------------------------------------------------------------
1789 * Do a row check, perhaps after recomputing the basic x.
1790 * ---------------------------------------------------------------
1791
1792 * Get the full row residuals y = rhs - (A I)*xn.
1793
1794 call daxpy ( m, (- one), xn(n+1), 1, y, 1 )
1795
1796 if (mode .eq. 1) then
1797 * ===============================================
1798 * Extract the basic and superbasic x from xn.
1799 * See if iterative refinement is worth doing.
1800 * ===============================================
1801 call m5bsx ( 2, ms, nb, kb, x, xn )
1802 ir = idamax( m, y, 1 )
1803 rmax = abs( y(ir) )
1804 ir = idamax( m, x, 1 )
1805 xnorm = abs( x(ir) )
1806 rowerr = rmax / (one + xnorm)
1807
1808 if (rowerr .gt. eps0) then
1809 * ===============================================
1810 * Compute a correction to basic x from B*y2 = y.
1811 * Set basic x = x + y2.
1812 * Store the new basic x into xn.
1813 * ===============================================
1814 call m2bsol( 2, m, y, y2, z, nwcore )
1815 call daxpy ( m, one, y2, 1, x, 1 )
1816 call m5bsx ( 1, m , nb, kb, x, xn )
1817
1818 * Compute y = rhs - (A I)*xn again for the new xn.
1819
1820 if (nncon .gt. 0) call dcopy ( nncon, z(lrhs), 1, y, 1 )
1821 if (nncon .lt. m) call dload ( m-nncon, zero,
1822 $ y(nncon+1), 1 )
1823 call m2aprd( 5, xn, n, y, m,
1824 $ ne, nka, a, ha, ka,
1825 $ z, nwcore )
1826 call daxpy ( m, (- one), xn(n+1), 1, y, 1 )
1827 end if
1828 end if
1829
1830 * Find the norm of x(BS), the basic and superbasic x.
1831 * Allow for the rare case x = 0.
1832
1833 xnorm = dnorm1( ms, x, 1 )
1834 if (xnorm .le. eps4) xnorm = one
1835
1836 * Find the maximum row residual.
1837
1838 ir = idamax( m, y, 1 )
1839 rmax = abs( y(ir) )
1840 rowerr = rmax / (one + xnorm)
1841 if (rowerr .gt. tolrow) ierr = 10
1842 if (rowerr .gt. tolrow .or. idebug .ge. 2) then
1843 if (iprint .gt. 0) write(iprint, 1000) itn, rmax, ir, xnorm
1844 end if
1845 end if
1846 return
1847
1848 1000 format(' Itn', i7, ' -- row check',
1849 $ '. Max residual =', 1p, e8.1, ' on row', i5,
1850 $ '. Norm x =', e9.2)
1851
1852 * end of m5setx
1853 end
1854
1855 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1856
1857 subroutine m5solv( m, maxr, maxs, mbs, n, nb, nn, nn0, nr,
1858 $ lcrash, ns, nscl, nx, objadd,
1859 $ ne, nka, a, ha, ka,
1860 $ hrtype, hs, kb, ascale, bl, bu,
1861 $ bbl, bbu, fsub, gsub, grd, grd2,
1862 $ pi, r, rc, rg, rg2,
1863 $ x, xn, y, y2, z, nwcore )
1864
1865 implicit double precision (a-h,o-z)
1866 integer*4 ha(ne), hrtype(mbs), hs(nb)
1867 integer ka(nka), kb(mbs)
1868 double precision a(ne), ascale(nscl), bl(nb), bu(nb)
1869 double precision bbl(mbs), bbu(mbs),
1870 $ gsub(nn0), grd(mbs), grd2(mbs),
1871 $ pi(m), r(nr), rc(nb), rg(maxs), rg2(maxs),
1872 $ x(nx), xn(nb), y(nx), y2(nx), z(nwcore)
1873
1874 * ------------------------------------------------------------------
1875 * m5solv solves the current problem. A basis is assumed to be
1876 * specified by ns, hs(*), xn(*) and the superbasic parts of kb(*).
1877 * In particular, there must be ns values hs(j) = 2, and the
1878 * corresponding j's must be listed in kb(m+1) thru kb(m+ns).
1879 * The ordering in kb matches the projected Hessian R (if any).
1880 *
1881 * GOTFAC and GOTHES are input and output variables.
1882 * On entry they say whether to preserve LU and R. For the first
1883 * cycle they will normally be false.
1884 * On exit they say whether LU and R exist, for possible use on
1885 * on later cycles.
1886 * gotlam no longer exists. Instead, initial estimates of the
1887 * Lagrange multipliers for the nonlinear constraints are
1888 * always extracted from pi.
1889 *
1890 * On exit, if ierr .lt. 30 it is safe to save the final
1891 * basis files and print the solution. Otherwise, a fatal error
1892 * condition exists and numerous items will be undefined.
1893 * The last basis map saved (if any) retains the only useful
1894 * information.
1895 *
1896 * 31 Jan 1992: rc(*) introduced as workspace for m5pric.
1897 * 06 Mar 1992: Major iteration 1 now just satisfies the linear
1898 * constraints and bounds. It may take as long as it
1899 * likes -- the minor iterations limit does not apply.
1900 * 17 Mar 1992: m7fixb implemented to try and fix an ill-conditioned
1901 * basis by swapping a column B with a column of S.
1902 * In this version, we consider it only at the start of
1903 * a major iteration.
1904 * fixb says when we want to try.
1905 * 08 Apr 1992: m5hs implemented to switch nonbasic hs(j) between
1906 * external and internal values to help speed up m5pric.
1907 * 10 Apr 1992: objadd added as an input parameter.
1908 * 04 Jun 1992: Crash level 3 implemented.
1909 * 07 Oct 1992: Minor iterations limit now restricts each major itn
1910 * to nminor FEASIBLE iterations.
1911 * ------------------------------------------------------------------
1912
1913 logical conv,restrt
1914 logical alone, AMPL, GAMS, MINT, page1, page2
1915 common /m1env / alone, AMPL, GAMS, MINT, page1, page2
1916 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1917 common /m1file/ iread,iprint,isumm
1918 common /m2file/ iback,idump,iload,imps,inewb,insrt,
1919 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
1920 common /m2lu3 / lenl,lenu,ncp,lrow,lcol
1921 common /m2parm/ dparm(30),iparm(30)
1922 common /m3mps2/ lcnam1,lrnam1,lcnam2,lrnam2,lxs,lxl,lfree
1923 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1924 common /m3scal/ sclobj,scltol,lscale
1925 common /m5freq/ kchk,kinv,ksav,klog,ksumm,i1freq,i2freq,msoln
1926 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1927 common /m5log1/ idebug,ierr,lprint
1928 common /m5log2/ jq1,jq2,jr1,jr2,lines1,lines2
1929 common /m5log3/ djq,theta,pivot,cond,nonopt,jp,jq,modr1,modr2
1930 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
1931 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
1932 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1933 common /m5lp2 / invrq,invitn,invmod
1934 common /m5prc / nparpr,nmulpr,kprc,newsb
1935 common /m5step/ featol, tolx0,tolinc,kdegen,ndegen,
1936 $ itnfix, nfix(2)
1937 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1938 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
1939 common /m7loc / lgobj ,lgobj2
1940 common /m7cg1 / cgbeta,itncg,msgcg,modcg,restrt
1941 common /m7conv/ etash,etarg,lvltol,nfail,conv(4)
1942 common /m7phes/ rgmin1,rgnrm1,rgnrm2,jz1,jz2,labz,nfullz,mfullz
1943 common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg
1944 common /m8len / njac ,nncon ,nncon0,nnjac
1945 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
1946 $ lblslk,lbuslk,lxlam ,lrhs ,
1947 $ lgcon ,lgcon2,lxdif ,lxold
1948 common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor
1949 common /m8al2 / radius,rhsmod,modpen,modrhs
1950 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
1951 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
1952 common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave
1953 logical gotbas,gotfac,gothes,gotscl
1954 common /cycle1/ gotbas,gotfac,gothes,gotscl
1955 common /cycle2/ objtru,suminf,numinf
1956
1957 character*12 istate
1958 logical first , fixb , incres,
1959 $ feasbl, infsbl, jstfes,
1960 $ linear, nonlin, lincon, nlncon,
1961 $ major1, convgd, optsub
1962
1963 intrinsic abs, min, max, mod
1964
1965 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1966
1967 * ------------------------------------------------------------------
1968
1969 * Initialize scalars (in alphabetical order).
1970
1971 fixb = .false.
1972 linear = nn .eq. 0
1973 nonlin = nn .gt. 0
1974 lincon = nncon .eq. 0
1975 nlncon = nncon .gt. 0
1976
1977 cond = zero
1978 flin = zero
1979 fobj = zero
1980 fsub = zero
1981 fxprev = zero
1982
1983 ierr = 0
1984 invitn = 0
1985 invrq = 0
1986 itn = 0
1987 lines1 = 0
1988 lines2 = 0
1989 majits = 0
1990 minits = 0
1991 minfes = 0
1992 modwt = 0
1993 m1 = m + 1
1994 nfac = 0
1995 nfail = 0
1996 ninf = 0
1997 nphs = 1
1998 nsamef = 0
1999 msamef = max( 2*nb, 200 )
2000 msamef = min( msamef, 1000 )
2001 nssave = ns
2002
2003 objtru = zero
2004 penpar = dparm(7)
2005 pinorm = one
2006 rgnorm = zero
2007 rgtest = zero
2008 sinf = zero
2009 tolz = eps0
2010 xnorm = one
2011 do 30 i = 1, 4
2012 conv(i) = .false.
2013 30 continue
2014
2015 newhed = .true.
2016 if (.not. gothes) r(1) = zero
2017
2018 if ( lincon ) then
2019
2020 * Set nonbasic hs(j) to the correct internal values.
2021
2022 call m5hs ( 'Internal', nb, bl, bu, hs, xn )
2023 else
2024 * ( nlncon )
2025 infsub = 0
2026 lcstat = 0
2027 nomove = 0
2028 nreduc = 0
2029 gotfac = .false.
2030 optsub = .false.
2031 end if
2032
2033 * ------------------------------------------------------------------
2034 * Start of a Major Iteration.
2035 * "first" says we've just started, and tells m5frmc to set nphs.
2036 * If the problem is linearly constrained, there is only one
2037 * major iteration (only one subproblem to be solved).
2038 * ------------------------------------------------------------------
2039 90 majits = majits + 1
2040 major1 = majits .eq. 1
2041 first = .true.
2042 ierr = 0
2043 kprc = 0
2044 lvldif = 1
2045 lvltol = 1
2046 ms = m + ns
2047 nfail = 0
2048 nfullz = 0
2049 nonopt = 0
2050 toldj(1) = plinfy
2051 toldj(2) = plinfy
2052 tolrg = zero
2053 call m5dgen( 1, m, n, nb, ms, inform,
2054 $ ne, nka, a, ha, ka,
2055 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
2056
2057 if (nlncon) then
2058 call m8setj( lcstat, m, n, nb, nn, nncon, nnjac, njac, nscl,
2059 $ lcrash, ns,
2060 $ nconvg, nomove, nreduc, optsub, convgd,
2061 $ ne, nka, a, ha, ka, hs, ascale, bl, bu,
2062 $ z(lblslk),z(lbuslk),z(lfcon), z(lfcon2),
2063 $ z(lfdif), z(lfold), z(lgcon), z(lgcon2), z(lxlam),
2064 $ z(lrhs ), z(lxdif), z(lxold),
2065 $ pi, xn, y, z, nwcore )
2066
2067 minits = 0
2068 minfes = 0
2069 fixb = ns .gt. 0
2070 call m1envt( 999 )
2071 if (ierr .ne. 0 ) go to 900
2072 if ( convgd ) go to 800
2073 if (majits .gt. nmajor) go to 825
2074 if (itn .ge. itnlim) go to 830
2075 if (nomove .ge. 4 ) go to 850
2076 optsub = .false.
2077 end if
2078
2079 *-- if (major1) then
2080
2081 * Set the slacks s using the current structurals x.
2082 * If lincon, try and satisfy Ax + s = 0.
2083 * If nlncon, try and satisfy Ax + s = rhs.
2084
2085 *-- call m5setx( 0, m, n, nb, ms, kb,
2086 *-- $ ne, nka, a, ha, ka,
2087 *-- $ bl, bu, x, xn, y, y2, z, nwcore )
2088 *-- end if
2089
2090 * ------------------------------------------------------------------
2091 * Factorize the basis ( B = L*U ). If gotfac is true on entry to
2092 * m5solv, the first call to m2bfac will try to use existing factors.
2093 * ------------------------------------------------------------------
2094 100 if (first) then
2095 * relax
2096 else
2097 * Safeguard against system error: If there have been no
2098 * iterations since the last factorize, something is haywire.
2099 if (invitn .eq. 0) go to 865
2100 end if
2101
2102 105 call m2bfac( gotfac, nfac, m, mbs, n, nb, nr, nn, ns,
2103 $ lcrash, fsub, objadd,
2104 $ ne, nka, a, ha, ka,
2105 $ kb, hs, bl, bu, bbl, bbu,
2106 $ r, grd, x, xn, y, y2, z, nwcore )
2107 if (ierr .gt. 0) go to 860
2108
2109 nssave = ns
2110 lusiz0 = lenl + lenu
2111 lumax = 2*lusiz0
2112 kmodpi = 1
2113 restrt = .true.
2114 if (mod(lprint,10) .gt. 0) newhed = .true.
2115
2116 * Under certain circumstances, see if B is ill-conditioned
2117 * and perhaps improve it by swapping in a superbasic.
2118
2119 110 if ( fixb ) then
2120 fixb = .false.
2121 ms = m + ns
2122 call m7fixb( m, maxr, ms, n, nb, nr, ns, nx, inform,
2123 $ ne, nka, a, ha, ka,
2124 $ hs, kb, bl, bu, bbl, bbu,
2125 $ r, x, y, y2, z, nwcore )
2126
2127 * Refactorize if m7fixb tried to update B = L*U and failed.
2128
2129 if (inform .eq. 5) go to 105
2130 end if
2131
2132 * ------------------------------------------------------------------
2133 * MAIN ITERATION LOOP.
2134 * ------------------------------------------------------------------
2135 200 ierr = 0
2136 ms = m + ns
2137 ninf0 = ninf
2138 sinf0 = sinf
2139
2140 * Increment featol every iteration.
2141 * Every kdegen iterations, reset featol and move nonbasic variables
2142 * onto their bounds if they are currently very close.
2143
2144 featol = featol + tolinc
2145 if (mod( minits, kdegen ) .eq. 0) then
2146 if (minits .gt. 0) then
2147 call m5dgen( 2, m, n, nb, ms, inform,
2148 $ ne, nka, a, ha, ka,
2149 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
2150 end if
2151 end if
2152
2153 * ------------------------------------------------------------------
2154 * Check feasibility.
2155 * ------------------------------------------------------------------
2156 call m5frmc( n, nb, nn, ns, ms, maxs,
2157 $ lcrash, first, fsub, featol, objadd,
2158 $ ne, nka, a, ha, ka,
2159 $ bl, bu, bbl, bbu, hrtype, hs, kb,
2160 $ grd2, x, xn, z, nwcore )
2161
2162 if (ierr .ne. 0) go to 900
2163 feasbl = ninf .eq. 0
2164 infsbl = .not. feasbl
2165 jstfes = feasbl .and. (first .or. ninf0 .gt. 0)
2166 first = .false.
2167
2168 if (infsbl) then
2169
2170 * Reduce wtobj if sinf has increased noticeably.
2171
2172 if (wtobj .gt. zero) then
2173 if (itn .gt. 0) then
2174 if (sinf .gt. sinf0 * 1.1) go to 655
2175 end if
2176 end if
2177 end if
2178
2179 * ------------------------------------------------------------------
2180 * See if the next phase of Crash is needed.
2181 * ------------------------------------------------------------------
2182 if (jstfes .and. lcrash .gt. 0) then
2183 if (lcrash .eq. 1) then
2184
2185 * All constraints have now been satisfied.
2186 * Major1 will be terminated below after pi is set.
2187
2188 lcrash = 0
2189
2190 else if (lcrash .eq. 2 .or. lcrash .eq. 4) then
2191
2192 * Linear constraints have been satisfied.
2193 * Major1 will be terminated and Major2 will
2194 * call Crash on nonlinear rows.
2195
2196 if (lincon) lcrash = 0
2197 if (nlncon) lcrash = 5
2198
2199 else if (lcrash .eq. 3) then
2200 if (nncon .lt. m) then
2201
2202 * So far the linear inequality constraints (if any)
2203 * have been ignored. Do another Crash and continue Major1.
2204 * m2amat sets row-types.
2205
2206 lcrash = 4
2207 call m5hs ( 'External', nb, bl, bu, hs, xn )
2208 call m2amat( 2, m, n, nb,
2209 $ ne, nka, a, ha, ka,
2210 $ bl, bu, hrtype )
2211 call m2crsh( lcrash, m, n, nb, nn,
2212 $ ne, nka, a, ha, ka,
2213 $ kb, hs, hrtype, bl, bu, xn, z, nwcore )
2214 call m5hs ( 'Internal', nb, bl, bu, hs, xn )
2215 ninf = 1
2216 go to 105
2217
2218 else
2219 lcrash = 5
2220 end if
2221 end if
2222 end if
2223
2224 * ------------------------------------------------------------------
2225 * Compute pi.
2226 * ------------------------------------------------------------------
2227 if (kmodpi .gt. 0) then
2228 if (feasbl .and. nonlin) then
2229 call m6grd ( ms, nb, nn, gsub, grd2,
2230 $ ne, nka, a, ha, ka,
2231 $ xn, z, nwcore )
2232 if (ierr .ne. 0) go to 900
2233 end if
2234 call dcopy ( ms, grd2, 1, grd, 1 )
2235 call m5setp( 1, m, grd2, pi, z, nwcore )
2236 end if
2237
2238 * ------------------------------------------------------------------
2239 * If we have just got feasible, reset a few things.
2240 * If there are nonlinear constraints and it's Major1,
2241 * we have satisfied the linear constraints for the first time.
2242 * That is enough for now --- time to relinearize.
2243 * ------------------------------------------------------------------
2244 if (jstfes) then
2245 lvltol = 1
2246 tolrg = zero
2247 if (nlncon .and. major1) go to 500
2248 end if
2249
2250 * ------------------------------------------------------------------
2251 * Compute rg = reduced gradient.
2252 * ------------------------------------------------------------------
2253 rgnorm = zero
2254 if (nphs .gt. 2) then
2255 if (ns .gt. 0) then
2256 call m7rg ( m, ms, ns, grd, pi, rg, rgnorm,
2257 $ ne, nka, a, ha, ka,
2258 $ z, nwcore )
2259 end if
2260 if (tolrg .eq. zero) tolrg = etarg * rgnorm
2261 end if
2262
2263 * ------------------------------------------------------------------
2264 * Proceed with the next iteration.
2265 * ------------------------------------------------------------------
2266 250 kmodlu = 1
2267 kmodpi = 1
2268 modr1 = 0
2269 modr2 = 0
2270 jq1 = 0
2271 jq2 = 0
2272 jr1 = 0
2273 jr2 = 0
2274 jz2 = 0
2275 labz = 1
2276 nxtphs = nphs
2277 if (jobj .gt. 0) flin = - xn(jobj) * sclobj + objadd
2278 objtru = flin + fobj
2279
2280 if (nphs .le. 3) then
2281
2282 * Phase 1, 2 or 3.
2283 * PRICE:
2284 * Select one or more nonbasic variables to become superbasic.
2285
2286 if (ns .eq. maxs) go to 840
2287
2288 call m5pric( m, mbs, n, nb, nn, nn0, ns, maxr, maxs, incres,
2289 $ ne, nka, a, ha, ka,
2290 $ hs, kb, bl, bu, gsub,
2291 $ pi, rc, rg )
2292
2293 if (nonopt .eq. 0) go to 400
2294 if (newsb .gt. 1) nphs = 3
2295 if ( nlncon ) nconvg = 0
2296 lvldif = 1
2297 lvltol = 1
2298 nsubsp = 0
2299 end if
2300
2301 * Test for excess iterations or time.
2302
2303 260 if (lincon) then
2304 if (itn .ge. itnlim) go to 830
2305 if (mod(itn,10) .eq. 0) call m1envt( 999 )
2306 else
2307 if (itn .ge. itnlim) go to 500
2308
2309 * The first major iteration may take as long as it likes
2310 * to satisfy the linear constraints. Later major iterations
2311 * are terminated by the minor iterations limit.
2312 * 07 Oct 1992: Let all majors take as long as they like
2313 * to get feasible. Use minfes instead of minits.
2314
2315 if (minfes .ge. nminor) go to 500
2316 end if
2317 if (ierr .gt. 0) go to 880
2318
2319 * ==================================================================
2320 * Perform a minor iteration (either LP or reduced-gradient method).
2321 * ==================================================================
2322 280 if (nphs .le. 2) then
2323 call m5lpit( m, m1, n, nb, incres,
2324 $ ne, nka, a, ha, ka,
2325 $ hrtype, hs, kb,
2326 $ bl, bu, bbl, bbu, x, xn, y, y2, z, nwcore )
2327 else
2328 if (nphs .eq. 4 .and. ns .le. 0) then
2329 nphs = 3
2330 go to 250
2331 end if
2332
2333 call m7rgit( m, maxr, maxs, mbs, n, nb, incres,
2334 $ nn, nn0, nr, ns, nx, inform, nxtphs,
2335 $ ne, nka, a, ha, ka,
2336 $ hrtype, hs, kb, bl, bu, bbl, bbu,
2337 $ fsub, gsub, grd, grd2,
2338 $ pi, r, rg, rg2, x, xn, y, y2, z, nwcore )
2339
2340 * For certain values of inform we want to try again.
2341 * If inform = 6, m7fixb has been requested (to fix R).
2342
2343 if (inform .eq. 6) then
2344 fixb = .true.
2345 if (invitn .gt. 0) go to 105
2346 go to 110
2347 end if
2348
2349 * If inform = 5, m2bfac has been requested
2350 * (part of error recovery).
2351 * If inform = 4, we want a better gsub and pi using
2352 * central differences (non-derivative case only).
2353 * If inform = 1, we want to call m5pric.
2354
2355 if (inform .eq. 5) go to 100
2356 if (inform .eq. 4) go to 200
2357 if (inform .eq. 1) go to 250
2358
2359 * Otherwise, m7rgit exits with a good pi and grd
2360 * (unless ierr gt 0).
2361 * Set kmodpi to indicate that we don't need a new pi.
2362
2363 kmodpi = 0
2364 end if
2365
2366 * ==================================================================
2367 * Test for error condition and/or frequency interrupts.
2368 * ==================================================================
2369 300 ms = m + ns
2370 if (ierr .eq. 2) go to 460
2371 if (ierr .gt. 0) go to 680
2372 itn = itn + 1
2373 invitn = invitn + 1
2374 minits = minits + 1
2375 if (ninf .eq. 0) minfes = minfes + 1
2376 if (linear .and. iobj .ne. 0) fsub = - minimz * x(iobj) * sclobj
2377
2378 * Print the iteration log.
2379
2380 call m5log ( m, maxs, mbs, n, nb, nn, ns, nx, fsub, objadd,
2381 $ ne, nka, a, ha, ka,
2382 $ bl, bu, hs, kb, x, xn, y, z, nwcore )
2383
2384 * Check for the dreaded infinite loop
2385 * by seeing if the objective has stayed the same much too long.
2386
2387 nsamef = nsamef + 1
2388 fx = sinf
2389 if ( feasbl ) fx = fsub
2390 if (abs( fxprev - fx ) .gt. (one + abs( fx )) * eps1) nsamef = 0
2391 if (nsamef .ge. msamef) go to 835
2392
2393 * Press on.
2394
2395 fxprev = fx
2396 nphs = nxtphs
2397 if (nlncon) then
2398 if (infsbl .or. rgnorm .gt. rgtest) nconvg = 0
2399 end if
2400
2401 * Save a basis map (frequency controlled).
2402 * We have to convert hs(*) to External and then back to Internal.
2403
2404 if (mod(itn,ksav) .eq. 0) then
2405 if (inewb .gt. 0 .and. itn .lt. itnlim) then
2406 call m5hs ( 'External', nb, bl, bu, hs, xn )
2407 call m4stat( 0, istate )
2408 call m4newb( 1, inewb, m, n, nb, nn, ns, ms, nscl, fsub,
2409 $ kb, hs, ascale, bl, bu, x, xn, istate )
2410 if (iback .gt. 0)
2411 $ call m4newb( 1, iback, m, n, nb, nn, ns, ms, nscl, fsub,
2412 $ kb, hs, ascale, bl, bu, x, xn, istate )
2413 call m5hs ( 'Internal', nb, bl, bu, hs, xn )
2414 end if
2415 end if
2416
2417 * Refactorize the basis if it has been modified a lot.
2418
2419 if (invmod .ge. kinv - 1) then
2420 invrq = 1
2421 go to 100
2422 end if
2423
2424 lusize = lenl + lenu
2425 if (invmod .ge. 20 .and. lusize .gt. lumax) then
2426 bgrwth = lusize
2427 bold = lusiz0
2428 bgrwth = bgrwth / bold
2429 if (prnt1) write(iprint, 1030) bgrwth
2430 invrq = 2
2431 go to 100
2432 end if
2433
2434 * Update the LU factors of the basis if requested.
2435
2436 if (kmodlu .eq. 1) call m2bsol( 4, m, y2, y, z, nwcore )
2437 if (invrq .ne. 0) go to 100
2438
2439 * Check row error (frequency controlled).
2440
2441 if (mod(invitn,kchk) .eq. 0) then
2442 call m5setx( 2, m, n, nb, ms, kb,
2443 $ ne, nka, a, ha, ka,
2444 $ bl, bu, x, xn, y, y2, z, nwcore )
2445 if (ierr .gt. 0) then
2446 invrq = 10
2447 go to 100
2448 end if
2449 end if
2450
2451 go to 200
2452
2453 * ------------------------------------------------------------------
2454 * END OF MAIN ITERATION LOOP.
2455 * ------------------------------------------------------------------
2456
2457 * The solution is apparently optimal -- check rgtols, etc.
2458
2459 400 if ( infsbl ) then
2460 * =============================
2461 * Optimal but still infeasible.
2462 * =============================
2463 if (wtobj .ne. zero ) go to 655
2464
2465 * See if any nonbasics have to be set back on their bound.
2466
2467 call m5dgen( 3, m, n, nb, ms, inform,
2468 $ ne, nka, a, ha, ka,
2469 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
2470 if (inform .gt. 0 ) go to 200
2471
2472 if (ns .eq. 0 ) go to 410
2473 if (lvltol .eq. 2 ) go to 410
2474 lvltol = 2
2475 tolrg = toldj(3) * pinorm
2476 if (rgnorm .le. tolrg) go to 410
2477
2478 nphs = 4
2479 if (prnt1) write(iprint, 1400) tolrg, lvltol
2480 go to 280
2481
2482 * The subproblem is infeasible.
2483 * Stop if the constraints are linear
2484 * or if this is the first major iteration
2485 * (which could not satisfy the linear constraints).
2486
2487 410 lcstat = 1
2488 if (lincon .or. major1) go to 810
2489 if (modrhs .eq. 0) infsub = infsub + 1
2490 if (infsub .gt. 5) go to 810
2491
2492 * Slightly relax the bounds on the linearized constraints.
2493
2494 call m8augl( 5, m, n, nb, ns, inform,
2495 $ ne, nka, a, ha, ka,
2496 $ hs, bl, bu, xn, z, nwcore )
2497
2498 if (inform .eq. 0) then
2499
2500 * Reset nonbasic hs(j) to internal values.
2501 * Reset the basic x's and check the row error.
2502
2503 kmodpi = 1
2504 call m5hs ( 'Internal', nb, bl, bu, hs, xn )
2505 call m5setx( 1, m, n, nb, ms, kb,
2506 $ ne, nka, a, ha, ka,
2507 $ bl, bu, x, xn, y, y2, z, nwcore )
2508 if (ierr .gt. 0) then
2509 invrq = 10
2510 go to 100
2511 end if
2512 go to 200
2513 end if
2514
2515 * The rhs has been relaxed often enough for this subproblem.
2516 * Keep going if it's the first time or something is happening.
2517
2518 if (infsub .le. 1 .or. minits .gt. 0) go to 500
2519 go to 810
2520 else
2521 * =============================================
2522 * Feasible and maybe optimal.
2523 * m5pric says there are no nonbasics to bring in.
2524 * We can stop if rgnorm is suitably small.
2525 * For Partial Completion (ncom = 0), a looser
2526 * measure of small will do.
2527 * =============================================
2528 if (ns .eq. 0) go to 420
2529 if (ncom .eq. 0) then
2530 tolmin = 1.0d-2 * pinorm
2531 else
2532 tolmin = toldj(3) * pinorm
2533 end if
2534 if (rgnorm .le. tolmin) go to 420
2535 if (lvltol .eq. 2 ) go to 420
2536 if (nfail .gt. 0 .and.
2537 $ nfullz .le. 1 ) go to 420
2538
2539 nfail = 0
2540 nphs = 4
2541 tolrg = 0.1 * min( tolrg, rgnorm )
2542 if (tolrg .le. tolmin) then
2543 lvltol = 2
2544 tolrg = tolmin
2545 end if
2546 if (prnt1) write(iprint, 1400) tolrg, lvltol
2547 go to 280
2548
2549 * The rgtols are small enough.
2550 * See if any nonbasics have to be set back on their bound.
2551
2552 420 call m5dgen( 3, m, n, nb, ms, inform,
2553 $ ne, nka, a, ha, ka,
2554 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
2555 if (inform .gt. 0) go to 200
2556
2557 * Check residuals before terminating.
2558
2559 djq = minimz * djq
2560 if (prnt1) write(iprint, 1410) djq, jq, rgnorm, pinorm
2561
2562 if (invitn .gt. 0) then
2563 call m5setx( 3, m, n, nb, ms, kb,
2564 $ ne, nka, a, ha, ka,
2565 $ bl, bu, x, xn, y, y2, z, nwcore )
2566 if (ierr .gt. 0) then
2567 invrq = 10
2568 go to 100
2569 end if
2570 end if
2571
2572 * The subproblem is optimal.
2573
2574 optsub = .true.
2575 lcstat = 0
2576 if ( nlncon ) then
2577 if (prnt1 ) write(iprint, 1510) majits, minits, itn
2578 if (summ1 ) write(isumm , 1511) minits, itn
2579 rgtest = 100.0 * toldj(3) * pinorm
2580 if (modrhs .eq. 0) infsub = 0
2581 if (modrhs .gt. 0) lcstat = 1
2582 nphs = 4
2583 go to 90
2584 end if
2585 go to 800
2586 end if
2587
2588 * The subproblem is unbounded.
2589 * Reduce wtobj or increase the penalty parameter.
2590
2591 460 if (infsbl .and. wtobj .gt. zero) go to 650
2592 call m5dgen( 3, m, n, nb, ms, inform,
2593 $ ne, nka, a, ha, ka,
2594 $ hs, kb, bl, bu, x, xn, y, y2, z, nwcore )
2595 if (inform .gt. 0) go to 200
2596
2597 lcstat = 2
2598 if ( nlncon ) then
2599 call m8augl( 4, m, n, nb, ns, inform,
2600 $ ne, nka, a, ha, ka,
2601 $ hs, bl, bu, xn, z, nwcore )
2602 if (inform .eq. 0) then
2603 kmodpi = 1
2604 ninf = 1
2605 go to 200
2606 end if
2607 end if
2608 go to 820
2609
2610 * Minor iterations terminated.
2611
2612 500 lcstat = 3
2613 if (prnt1) write(iprint, 1500) majits, minits, itn
2614 if (summ1) write(isumm , 1501) minits, itn
2615 nphs = 4
2616 go to 90
2617
2618 * Reset or reduce wtobj.
2619
2620 650 wtobj = zero
2621 655 modwt = modwt + 1
2622 wtobj = 0.1 * wtobj
2623 if (modwt .gt. 5) wtobj = zero
2624 if (prnt1 .or. idebug .ge. 1) write(iprint, 1200) wtobj
2625 go to 200
2626
2627 * Error condition.
2628 * If it was a linesearch failure, keep the major iterations
2629 * going if this one did at least one minor iteration.
2630 * Otherwise, print something helpful and terminate.
2631 * t1 is something that ought to be small at a solution.
2632 * t2 is what we consider 'reasonably small'.
2633
2634 680 if (ierr .ne. 9) go to 900
2635 t1 = rgnorm / pinorm
2636 t2 = 0.1
2637 if ( nlncon ) then
2638 if (minits .gt. 0) go to 500
2639 t1 = max( t1, virel )
2640 end if
2641
2642 if (iprint .gt. 0) write(iprint, 1104)
2643 if (isumm .gt. 0) write(isumm , 1104)
2644 if (lderiv .gt. 0 .and. .not. GAMS) then
2645 if (iprint .gt. 0) write(iprint, 1105)
2646 if (isumm .gt. 0) write(isumm , 1105)
2647 end if
2648 if (t1 .lt. t2) then
2649 if (iprint .gt. 0) write(iprint, 1114)
2650 if (isumm .gt. 0) write(isumm , 1114)
2651 else
2652 if (iprint .gt. 0) write(iprint, 1115)
2653 if (isumm .gt. 0) write(isumm , 1115)
2654 end if
2655 go to 850
2656
2657 * ------------------------------------------------------------------
2658 * Exit.
2659 * m1page( ) decides whether to do a page eject before the EXIT msg.
2660 * m1page(2) also calls m1envt(1) to print '=1' for GAMS
2661 * ------------------------------------------------------------------
2662
2663 * Optimal.
2664
2665 800 call m1page( 2 )
2666 rgtest = rgnorm / pinorm
2667 if (rgtest .lt. 0.1) then
2668 if (iprint .gt. 0) write(iprint, 1800)
2669 if (isumm .gt. 0) write(isumm , 1800)
2670 else
2671 if (iprint .gt. 0) write(iprint, 1802)
2672 if (isumm .gt. 0) write(isumm , 1802)
2673 ierr = 13
2674 end if
2675 go to 900
2676
2677 * Infeasible.
2678
2679 810 call m1page( 2 )
2680 if (iprint .gt. 0) write(iprint, 1810)
2681 if (isumm .gt. 0) write(isumm , 1810)
2682 ierr = 1
2683 go to 900
2684
2685 * Unbounded.
2686
2687 820 call m1page( 2 )
2688 if (iprint .gt. 0) write(iprint, 1820)
2689 if (isumm .gt. 0) write(isumm , 1820)
2690 go to 900
2691
2692 * Too many iterations.
2693
2694 825 call m1envt( 1 )
2695 if (iprint .gt. 0) write(iprint, 1825)
2696 call m1envt( 2 )
2697 830 call m1page( 2 )
2698 if (iprint .gt. 0) write(iprint, 1830)
2699 if (isumm .gt. 0) write(isumm , 1830)
2700 ierr = 3
2701 go to 900
2702
2703 * The objective hasn't changed for ages.
2704
2705 835 call m1page( 2 )
2706 if (iprint .gt. 0) write(iprint, 1835) msamef
2707 if (isumm .gt. 0) write(isumm , 1835) msamef
2708 ierr = 4
2709 go to 900
2710
2711 * maxs too small.
2712
2713 840 call m1page( 2 )
2714 if (iprint .gt. 0) write(iprint, 1840) maxs
2715 if (isumm .gt. 0) write(isumm , 1840) maxs
2716 ierr = 5
2717 go to 900
2718
2719 * Linesearch failure.
2720
2721 850 call m1page( 2 )
2722 if (iprint .gt. 0) write(iprint, 1850)
2723 if (isumm .gt. 0) write(isumm , 1850)
2724 ierr = 9
2725 go to 900
2726
2727 * Fatal error after LU factorization.
2728
2729 860 if (ierr .ne. 10) go to 900
2730 call m1page( 2 )
2731 if (iprint .gt. 0) write(iprint, 1860)
2732 if (isumm .gt. 0) write(isumm , 1860)
2733 go to 900
2734
2735 * m2bfac called twice in a row.
2736
2737 865 call m1page( 2 )
2738 if (iprint .gt. 0) write(iprint, 1865)
2739 if (isumm .gt. 0) write(isumm , 1865)
2740 ierr = 12
2741 go to 900
2742
2743 * Resource interrupt (too much time, etc.)
2744
2745 880 call m1envt( 1 )
2746 if (iprint .gt. 0) write(iprint, 1880)
2747 if (isumm .gt. 0) write(isumm , 1880)
2748 ierr = 19
2749 go to 900
2750
2751 * ------------------------------------------------------------------
2752 * Set output variables and print a summary of the final solution.
2753 * M1ENVT( 2 ) prints '=2' for GAMS.
2754 * ------------------------------------------------------------------
2755
2756 900 call m1envt( 2 )
2757 gotfac = ierr .le. 9
2758 gothes = r(1) .ne. zero
2759 degen = 100.0 * ndegen / max( itn, 1 )
2760 infsbl = ninf .gt. 0
2761 suminf = sinf
2762 numinf = ninf
2763 rgtest = rgnorm / pinorm
2764 xnorm = dnorm1( nb, xn, 1 )
2765
2766 * Count basic nonlinear variables (not used for anything).
2767
2768 nnb = 0
2769 do 910 j = 1, nn
2770 if (hs(j) .eq. 3) nnb = nnb + 1
2771 910 continue
2772
2773 * Restore nonlinear constraints bounds, in case they were relaxed.
2774 * Maybe reduce majits for printing purposes.
2775
2776 if (nlncon) then
2777 call m8augl( 6, m, n, nb, ns, inform,
2778 $ ne, nka, a, ha, ka,
2779 $ hs, bl, bu, xn, z, nwcore )
2780
2781 if (minits .eq. 0) majits = majits - 1
2782 end if
2783
2784 * Set nonbasic hs(j) to external values.
2785 * Save the final basis file.
2786
2787 call m5hs ( 'External', nb, bl, bu, hs, xn )
2788 if (inewb .gt. 0 .and. ierr .lt. 20) then
2789 k = ierr + 1
2790 call m4stat( k, istate )
2791 call m4newb( 2, inewb, m, n, nb, nn, ns, ms, nscl, fsub,
2792 $ kb, hs, ascale, bl, bu, x, xn, istate )
2793 end if
2794
2795 * Print statistics.
2796
2797 if (iprint .gt. 0) then
2798 write(iprint, 1900) name, itn, objtru
2799 if (infsbl) write(iprint, 1910) ninf, sinf
2800 if (nonlin) write(iprint, 1920) majits, flin, penpar, fobj
2801 if (nonlin) write(iprint, 1950) nfobj(1), nfcon(1)
2802 if (nonlin .and. lderiv .lt. 3)
2803 $ write(iprint, 1960) (nfobj(i), nfcon(i), i=2,4)
2804 if (nonlin .or. ns .gt. 0)
2805 $ write(iprint, 1970) ns, rgnorm, nnb, rgtest
2806 write(iprint, 1975) ndegen, degen
2807 end if
2808
2809 if (isumm .gt. 0) then
2810 write(isumm , 1900) name, itn, objtru
2811 if (infsbl) write(isumm , 1910) ninf, sinf
2812 if (nonlin) write(isumm , 1920) majits, flin, penpar, fobj
2813 if (nonlin) write(isumm , 1950) nfobj(1), nfcon(1)
2814 if (nonlin .and. lderiv .lt. 3)
2815 $ write(isumm , 1960) (nfobj(i), nfcon(i), i=2,4)
2816 if (nonlin .or. ns .gt. 0)
2817 $ write(isumm , 1970) ns, rgnorm, nnb, rgtest
2818 write(isumm , 1975) ndegen, degen
2819 end if
2820
2821 950 return
2822
2823
2824 1030 format( ' LU file has increased by a factor of', f6.1)
2825 1104 format(/ ' XXX The functions are non-smooth or very nonlinear')
2826 1105 format( ' XXX MAKE SURE THE GRADIENTS ARE CODED CORRECTLY')
2827 1114 format( ' XXX The final point is REASONABLY NEAR AN OPTIMUM' /)
2828 1115 format( ' XXX The final point is NOT CLOSE TO AN OPTIMUM' /)
2829 1200 format(/ ' Weight on obj reduced to', 1p, e12.2)
2830 1400 format( ' tolrg reduced to', 1p, e11.3, 5x, ' lvltol =', i2)
2831 1410 format(/ ' Biggest dj =', 1p, e11.3, ' (variable', i7, ')',
2832 $ ' norm rg =', e11.3, ' norm pi =', e11.3)
2833 1500 format(/ ' End of major itn', i4,
2834 $ ' - minor itns terminated at', i4,
2835 $ ' - total itns =', i7)
2836 1501 format(' Minor itns terminated at', i4,
2837 $ ' - total itns =', i7)
2838 1510 format(/ ' End of major itn', i4,
2839 $ ' - optimal subproblem at minor itn', i4,
2840 $ ' - total itns =', i7)
2841 1511 format(' Optimal subproblem at minor itn', i4,
2842 $ ' - total itns =', i7)
2843
2844 1800 format(' EXIT -- optimal solution found')
2845 1802 format(' EXIT -- near-optimal solution found'
2846 $ /' XXX WARNING: reduced gradient is large --'
2847 $ ' solution is not really optimal')
2848 1810 format(' EXIT -- the problem is infeasible')
2849 1820 format(' EXIT -- the problem is unbounded ',
2850 $ ' (or badly scaled).')
2851 1825 format(' XXX Major iterations terminated before convergence')
2852 1830 format(' EXIT -- too many iterations')
2853 1835 format(' EXIT -- the objective has not changed',
2854 $ ' for the last', i8, ' iterations')
2855 1840 format(' EXIT -- the superbasics limit is too small:', i7)
2856 1850 format(' EXIT -- the current point cannot be improved upon')
2857 1860 format(' EXIT -- numerical error. General constraints',
2858 $ ' cannot be satisfied accurately')
2859 1865 format(' EXIT -- basis factorization requested',
2860 $ ' twice in a row')
2861 1880 format(' EXIT -- resource interrupt')
2862
2863 1900 format(/ 1x, 2a4
2864 $ / ' No. of iterations', i20,
2865 $ 2x, ' Objective value', 1p, e22.10)
2866 1910 format( ' No. of infeasibilities', i15,
2867 $ 2x, ' Sum of infeas', 1p, e24.10)
2868 1920 format( ' No. of major iterations', i14,
2869 $ 2x, ' Linear objective', 1p, e21.10
2870 $ / ' Penalty parameter', 0p, f20.6,
2871 $ 2x, ' Nonlinear objective', 1p, e18.10)
2872 1950 format( ' No. of calls to funobj', i15,
2873 $ 2x, ' No. of calls to funcon', i15)
2874 1960 format( ' Calls with mode=2 (f, known g)', i7,
2875 $ 2x, ' Calls with mode=2 (f, known g)', i7
2876 $ / ' Calls for forward differencing', i7,
2877 $ 2x, ' Calls for forward differencing', i7
2878 $ / ' Calls for central differencing', i7,
2879 $ 2x, ' Calls for central differencing', i7)
2880 1970 format( ' No. of superbasics', i19,
2881 $ 2x, ' Norm of reduced gradient', 1p, e13.3
2882 $ / ' No. of basic nonlinears', i14,
2883 $ 2x, ' Norm rg / Norm pi', e20.3)
2884 1975 format( ' No. of degenerate steps', i14,
2885 $ 2x, ' Percentage', f27.2)
2886
2887 * end of m5solv
2888 end

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