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

Contents of /trunk/minos54/mi40bfil.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (show annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (19 years, 7 months ago) by aw0a
File size: 60835 byte(s)
Setting up web subdirectory in repository
1 ************************************************************************
2 *
3 * File mi40bfil fortran.
4 *
5 * m4getb m4chek m4id m4name m4inst m4load m4oldb
6 * m4savb m4dump m4newb m4pnch m4rc m4infs
7 * m4rept m4soln m4solp m4stat
8 *
9 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
10
11 subroutine m4getb( ncycle, istart, m, mbs, n, nb, nn, nname, nscl,
12 $ lcrash, ns,
13 $ ne, nka, a, ha, ka,
14 $ hrtype, hs, kb, ascale, bl, bu,
15 $ pi, xn, y, y2, name1, name2, z, nwcore )
16
17 implicit double precision (a-h,o-z)
18 integer*4 ha(ne), hrtype(mbs), hs(nb)
19 integer ka(nka), kb(mbs), name1(nname), name2(nname)
20 double precision a(ne), ascale(nscl), bl(nb), bu(nb)
21 double precision pi(m), xn(nb), y(m), y2(m), z(nwcore)
22
23 * ------------------------------------------------------------------
24 * m4getb is called with ncycle = 0 before the first cycle.
25 * A basis file is loaded (if any exists).
26 *
27 * m4getb is called with ncycle > 0 before every cycle.
28 * The Jacobian is evaluated and stored in a (unscaled).
29 * If gotscl is false, m2scal is called to compute scales.
30 * If relevant, the scales are applied to a, bl, bu, xn, pi, fcon.
31 * If gotbas is false, m2crsh is called to initialize hs.
32 * (kb, y, y2 are used as workspace by m2crsh.)
33 *
34 * In both cases, lcrash is an output parameter to tell m5solv
35 * if further calls to crash are needed.
36 *
37 * 14 May 1992: pi(1:nncon) assumed to be initialized on entry.
38 * xn passed to m2crsh.
39 * 04 Jun 1992: lcrash added as an output parameter.
40 * 09 Jul 1992: istart added as an input parameter.
41 * Used when ncycle = 0 to load a basis file only if
42 * it is a cold start.
43 * ------------------------------------------------------------------
44
45 common /m1file/ iread,iprint,isumm
46 common /m2file/ iback,idump,iload,imps,inewb,insrt,
47 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
48 common /m2parm/ dparm(30),iparm(30)
49 common /m3scal/ sclobj,scltol,lscale
50 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
51 common /m5log1/ idebug,ierr,lprint
52 common /m8len / njac ,nncon ,nncon0,nnjac
53 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
54 $ lblslk,lbuslk,lxlam ,lrhs ,
55 $ lgcon ,lgcon2,lxdif ,lxold
56 common /m8diff/ difint(2),gdummy,lderiv,lvldif,knowng(2)
57 common /m8func/ nfcon(4),nfobj(4),nprob,nstat1,nstat2
58 common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave
59 logical gotbas,gotfac,gothes,gotscl
60 common /cycle1/ gotbas,gotfac,gothes,gotscl
61
62 intrinsic max, min
63 equivalence ( iparm(1), icrash )
64 logical gotjac
65 parameter ( zero = 0.0d+0, one = 1.0d+0 )
66
67 lcrash = 0
68
69 if (ncycle .eq. 0) then
70 * ---------------------------------------------------------------
71 * This call is made before the cycle loop.
72 * ---------------------------------------------------------------
73
74 * Initialize lvldif = 1, nfcon(*) = 0, nfobj(*) = 0.
75 * Also set linear pi(i) = 0 in case they are printed.
76 * We have to do this before scales are applied anyway.
77 * One day the basis files might load nonlinear pi(i).
78
79 lvldif = 1
80 call iload1( 4, 0, nfcon, 1 )
81 call iload1( 4, 0, nfobj, 1 )
82 if (nncon .lt. m) call dload ( m-nncon, zero, pi(nncon+1), 1 )
83
84 * Load a basis file if one exists and istart = 0 (Cold start).
85
86 if (istart .eq. 0) then
87 if (ioldb .gt. 0) then
88 call m4oldb( m, n, nb, ns, hs, bl, bu, xn )
89 else if (insrt .gt. 0) then
90 call m4inst( m, n, nb, ns, hs, bl, bu, xn, name1, name2 )
91 else if (iload .gt. 0) then
92 call m4load( m, n, nb, ns, hs, bl, bu, xn, name1, name2 )
93 end if
94 nssave = ns
95 end if
96
97 * Make sure the nonlinear variables are within bounds
98 * (for the gradient checkers if nothing else).
99
100 do 160 j = 1, nn
101 xn(j) = max( xn(j), bl(j) )
102 xn(j) = min( xn(j), bu(j) )
103 160 continue
104
105 else
106 * ---------------------------------------------------------------
107 * ncycle > 0. This call is made every cycle.
108 * ---------------------------------------------------------------
109 if (nncon .gt. 0) then
110
111 * Evaluate the Jacobian at the initial xn
112 * and store it in a(*) for the first linearization.
113 * We may have to copy some constant (unscaled) Jacobian
114 * elements from gcon2 into gcon.
115 * Also, we have to disable scaling temporarily, since if this
116 * is the first cycle, the scales have not yet been computed.
117
118 if (lscale .eq. 2 .and. lderiv .ge. 2) then
119 call dcopy ( njac, z(lgcon2), 1, z(lgcon), 1 )
120 end if
121 lssave = lscale
122 lscale = 0
123 gotjac = .false.
124 call m8ajac( gotjac, nncon, nnjac, njac,
125 $ ne, nka, a, ha, ka,
126 $ z(lfcon), z(lfcon2), z(lgcon), z(lgcon2),
127 $ xn, y, z, nwcore )
128 lscale = lssave
129 if (ierr .ne. 0) go to 900
130 end if
131
132 * Compute scales from a, bl, bu (unless we already have them).
133 * Then apply them to a, bl, bu, pi, xn and fcon.
134
135 if (lscale .gt. 0) then
136 if (.not. gotscl) then
137 gotscl = .true.
138 call m2scal( m, n, nb, ne, nka, nn, nncon, nnjac,
139 $ hrtype, ha, ka, a, ascale, bl, bu, y, y2 )
140 end if
141
142 call m2scla( 1, m, n, nb, ne, nka,
143 $ ha, ka, a, ascale, bl, bu, pi, xn )
144
145 if (lscale .eq. 2 .and. nncon .gt. 0) then
146 call dddiv ( nncon, ascale(n+1), 1, z(lfcon), 1 )
147 end if
148 end if
149
150 * Initialize xlam from pi. These are Lagrange multipliers
151 * for the nonlinear constraints.
152 * Do it at the very start, or on later cycles if the present
153 * solution is feasible. (ninf = 0 in both cases.
154 * Keep the previous xlam if solution is infeasible.)
155 * First change the sign of pi if maximizing.
156 *
157 * If pi(1:nncon) seems ridiculously big, assume that it
158 * has not been correctly initialized and just use xlam = 0.
159
160 if (ninf .eq. 0 .and. nncon .gt. 0) then
161 if (minimz .lt. 0) then
162 call dscal ( nncon, -one, pi, 1 )
163 end if
164 toobig = 1.0d+10
165 imax = idamax( nncon, pi, 1 )
166 xlmax = abs( pi(imax) )
167 if (xlmax .lt. toobig) then
168 call dcopy ( nncon, pi, 1, z(lxlam), 1 )
169 else
170 call dload ( nncon, zero, z(lxlam), 1 )
171 if (iprint .gt. 0) write(iprint, 1000)
172 if (isumm .gt. 0) write(isumm , 1000)
173 end if
174 end if
175
176 * ---------------------------------------------------------------
177 * If there was no basis file, find an initial basis via Crash.
178 * This works best with A scaled (hence much of the complication).
179 * xn(1:n) is input. xn(n+1:nb) is initialized by m2crsh.
180 * ---------------------------------------------------------------
181 if (.not. gotbas) then
182 gotbas = .true.
183 if (icrash .gt. 0 .and. icrash .le. 3) lcrash = icrash
184
185 call m2crsh( lcrash, m, n, nb, nn,
186 $ ne, nka, a, ha, ka,
187 $ kb, hs, hrtype, bl, bu, xn, z, nwcore )
188 end if
189 end if
190
191 * Exit.
192
193 900 return
194
195 1000 format(' XXX pi(1:nncon) is big (perhaps not initialized).'
196 $ / ' XXX Will set lambda = 0 for nonlinear rows.')
197 * end of m4getb
198 end
199
200 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
201
202 subroutine m4chek( m, maxs, mbs, n, nb, ns,
203 $ hs, kb, bl, bu, xn )
204
205 implicit double precision (a-h,o-z)
206 integer*4 hs(nb)
207 integer kb(mbs)
208 double precision bl(nb), bu(nb)
209 double precision xn(nb)
210
211 * ------------------------------------------------------------------
212 * m4chek takes hs and xn and checks they contain reasonable values.
213 * The entries hs(j) = 2 are used to set ns and nssave and possibly
214 * the list of superbasic variables kb(m+1) thru kb(m+ns).
215 * Scaling, if any, has taken place by this stage.
216 *
217 * If gotbas and gothes are both true, nssave and the superbasic kb's
218 * are assumed to be set. It must be a Hot start, or ncycle > 1.
219 * ------------------------------------------------------------------
220
221 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
222 common /m1file/ iread,iprint,isumm
223 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
224 common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave
225 logical gotbas,gotfac,gothes,gotscl
226 common /cycle1/ gotbas,gotfac,gothes,gotscl
227
228 intrinsic abs, max, min
229 logical setkb
230 parameter ( zero = 0.0d+0, tolb = 1.0d-4 )
231
232 * Make sure hs(j) = 0, 1, 2 or 3 only.
233
234 do 5 j = 1, nb
235 js = hs(j)
236 if (js .lt. 0) hs(j) = 0
237 if (js .ge. 4) hs(j) = js - 4
238 5 continue
239
240 setkb = .not. (gotbas .and. gothes)
241
242 * ------------------------------------------------------------------
243 * Make sure the objective is basic and free.
244 * Then count the basics and superbasics, making sure they don't
245 * exceed m and maxs respectively. Also, set ns and possibly
246 * kb(m+1) thru kb(m+ns) to define the list of superbasics.
247 * Mar 1988: Loop 100 now goes backwards to make sure we grab obj.
248 * Apr 1992: Backwards seems a bit silly in the documentation.
249 * We now go forward through the slacks,
250 * then forward through the columns.
251 * ------------------------------------------------------------------
252 10 nbasic = 0
253 ns = 0
254 if (iobj .gt. 0) then
255 jobj = n + iobj
256 hs(jobj) = 3
257 bl(jobj) = - plinfy
258 bu(jobj) = plinfy
259 end if
260
261 * If too many basics or superbasics, make them nonbasic.
262 * Do slacks first to make sure we grab the objective slack.
263
264 j = n
265
266 do 100 jj = 1, nb
267 j = j + 1
268 if (j .gt. nb) j = 1
269 js = hs(j)
270 if (js .eq. 2) then
271 ns = ns + 1
272 if (ns .le. maxs) then
273 if ( setkb ) kb(m + ns) = j
274 else
275 hs(j) = 0
276 end if
277
278 else if (js .eq. 3) then
279 nbasic = nbasic + 1
280 if (nbasic .gt. m) hs(j) = 0
281 end if
282 100 continue
283
284 * Proceed if the superbasic kbs were reset, or if ns seems to
285 * agree with nssave from the previous cycle.
286 * Otherwise, give up trying to save the projected Hessian, and
287 * reset the superbasic kbs after all.
288
289 if (setkb) then
290 * ok
291 else if (ns .ne. nssave) then
292 setkb = .true.
293 gothes = .false.
294 if (iprint .gt. 0) write(iprint, 1000) ns, nssave
295 if (isumm .gt. 0) write(isumm , 1000) ns, nssave
296 go to 10
297 end if
298
299 * Check the number of basics.
300
301 ns = min( ns, maxs )
302 nssave = ns
303 if (nbasic .ne. m ) then
304 gothes = .false.
305 if (iprint .gt. 0) write(iprint, 1100) nbasic, m
306 if (isumm .gt. 0) write(isumm , 1100) nbasic, m
307 end if
308
309 * -----------------------------------------------------------
310 * On all cycles, set each nonbasic xn(j) to be exactly on its
311 * nearest bound if it is within tolb of that bound.
312 * -----------------------------------------------------------
313 bplus = 0.1*plinfy
314 do 300 j = 1, nb
315 xj = xn(j)
316 if (abs( xj ) .ge. bplus) xj = zero
317 if (hs(j) .le. 1 ) then
318 b1 = bl(j)
319 b2 = bu(j)
320 xj = max( xj, b1 )
321 xj = min( xj, b2 )
322 if ((xj - b1) .gt. (b2 - xj)) b1 = b2
323 if (abs(xj - b1) .le. tolb) xj = b1
324 hs(j) = 0
325 if (xj .gt. bl(j)) hs(j) = 1
326 end if
327 xn(j) = xj
328 300 continue
329
330 return
331
332 1000 format(/ ' WARNING:', i6, ' superbasics in hs(*);',
333 $ ' previously ns =', i6, '. Hessian not saved')
334 1100 format(/ ' WARNING:', i7, ' basics specified;',
335 $ ' preferably should have been', i7)
336
337 * end of m4chek
338 end
339
340 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
341
342 subroutine m4id ( j, m, n, nb, nname, name1, name2, id1, id2 )
343
344 integer name1(nname), name2(nname)
345
346 * ------------------------------------------------------------------
347 * m4id returns a name id1-id2 for the j-th variable.
348 * If nname = nb, the name is already in name1, name2.
349 * Otherwise nname = 1. Some generic column or row name is cooked up
350 * ------------------------------------------------------------------
351
352 character*7 f1
353 character*5 f2
354 character*1 cname
355 character*1 rname
356 character*8 gname
357 data f1 /'(a1,i7)'/
358 data f2 /'(2a4)'/
359 data cname /'x'/
360 data rname /'r'/
361
362 if (nname .eq. nb) then
363 id1 = name1(j)
364 id2 = name2(j)
365 else if (j .le. n) then
366 write(gname, f1) cname, j
367 read (gname, f2) id1, id2
368 else
369 i = j - n
370 write(gname, f1) rname, i
371 read (gname, f2) id1, id2
372 end if
373
374 * end of m4id
375 end
376
377 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
378
379 subroutine m4name( n, name1, name2, id1, id2,
380 $ ncard, notfnd, maxmsg, j1, j2, jmark, jfound )
381
382 implicit double precision (a-h,o-z)
383 integer name1(n), name2(n)
384
385 * ------------------------------------------------------------------
386 * m4name searches for name1-name2 in arrays name1-2(j), j = j1, j2.
387 * jmark will probably speed the search on the next entry.
388 * Used by subroutines m3mpsc, m4inst, m4load.
389 *
390 * Left-justified alphanumeric data is being tested for a match.
391 * On Burroughs B6700-type systems, one could replace .eq. by .is.
392 * ------------------------------------------------------------------
393
394 common /m1file/ iread,iprint,isumm
395
396 do 50 j = jmark, j2
397 if (id1 .eq. name1(j) .and. id2 .eq. name2(j)) go to 100
398 50 continue
399
400 do 60 j = j1, jmark
401 if (id1 .eq. name1(j) .and. id2 .eq. name2(j)) go to 100
402 60 continue
403
404 * Not found.
405
406 jfound = 0
407 jmark = j1
408 notfnd = notfnd + 1
409 if (notfnd .le. maxmsg) then
410 if (iprint .gt. 0) write(iprint, 1000) ncard, id1, id2
411 end if
412 return
413
414 * Got it.
415
416 100 jfound = j
417 jmark = j
418 return
419
420 1000 format(' XXX Line', i6, ' -- name not found:', 8x, 2a4)
421
422 * end of m4name
423 end
424
425 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
426
427 subroutine m4inst( m, n, nb, ns,
428 $ hs, bl, bu, xn, name1, name2 )
429
430 implicit double precision (a-h,o-z)
431 integer*4 hs(nb)
432 integer name1(nb), name2(nb)
433 double precision bl(nb), bu(nb)
434 double precision xn(nb)
435
436 * ------------------------------------------------------------------
437 * This impression of INSERT reads a file produced by m4pnch.
438 * It is intended to read files similar to those produced by
439 * standard MPS systems. It recognizes SB as an additional key.
440 * Also, values are extracted from columns 25--36.
441 *
442 * 17 May 1992: John Stone (Ketron) mentioned trouble if rows and
443 * columns have the same name. The quick fix is to
444 * search column names from the beginning always,
445 * rather than from position jmark. Just set jmark = 1.
446 * ------------------------------------------------------------------
447
448 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
449 common /m1file/ iread,iprint,isumm
450 common /m2file/ iback,idump,iload,imps,inewb,insrt,
451 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
452 common /m3mps3/ aijtol,bstruc(2),mlst,mer,
453 $ aijmin,aijmax,na0,line,ier(20)
454 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
455
456 intrinsic abs
457 integer id(5)
458 character*4 key
459 character*4 lll , lul , lxl ,
460 $ lxu , lsb , lend
461 data lll /' LL '/, lul /' UL '/, lxl /' XL '/,
462 $ lxu /' XU '/, lsb /' SB '/, lend/'ENDA'/
463
464 bplus = 0.9*plinfy
465 if (iprint .gt. 0) write(iprint, 1999) insrt
466 if (isumm .gt. 0) write(isumm , 1999) insrt
467 read (insrt , 1000) id
468 if (iprint .gt. 0) write(iprint, 2000) id
469 l1 = n + 1
470
471 * Make logicals basic.
472
473 do 20 j = l1, nb
474 hs(j) = 3
475 20 continue
476
477 ignord = 0
478 nbs = 0
479 ns = 0
480 notfnd = 0
481 ncard = 0
482 * jmark = 1
483 lmark = l1
484 ndum = n + 100000
485
486 * Read names until ENDATA
487
488 do 300 nloop = 1, ndum
489 read(insrt, 1020) key, name1c, name2c, name1r, name2r, xj
490 if (key .eq. lend) go to 310
491
492 * Look for name1. It may be a column or a row,
493 * since a superbasic variable could be either.
494 * 17 May 1992: Set jmark = 1 so columns are searched first.
495 * This avoids trouble when columns and rows have the same name.
496
497 ncard = nloop
498 jmark = 1
499 call m4name( nb, name1, name2, name1c, name2c,
500 $ ncard, notfnd, mer, 1, nb, jmark, j )
501 if ( j .le. 0) go to 300
502 if (hs(j) .gt. 1) go to 290
503
504 if (key .eq. lxl .or. key .eq. lxu) then
505 * ------------------------------------------------------------
506 * XL, XU (exchange card) -- make col j basic, row l nonbasic.
507 * ------------------------------------------------------------
508
509 * Look for name2. It has to be a row.
510
511 call m4name( nb, name1, name2, name1r, name2r,
512 $ ncard, notfnd, mer, l1, nb, lmark, l )
513 if (l .le. 0 ) go to 300
514 if (l .eq. jobj) go to 290
515 if (hs(l) .ne. 3) go to 290
516
517 nbs = nbs + 1
518 hs(j) = 3
519 if (key .eq. lxl) then
520 hs(l) = 0
521 if (bl(l) .gt. -bplus) xn(l) = bl(l)
522 else
523 hs(l) = 1
524 if (bu(l) .lt. bplus) xn(l) = bu(l)
525 end if
526
527 * ---------------------------------------------------------------
528 * else LL, UL, SB -- only j and xj are relevant.
529 * ---------------------------------------------------------------
530 else if (key .eq. lll) then
531 hs(j) = 0
532 else if (key .eq. lul) then
533 hs(j) = 1
534 else if (key .eq. lsb) then
535 hs(j) = 2
536 ns = ns + 1
537 else
538 go to 290
539 end if
540
541 * Save xj.
542
543 if (abs(xj) .lt. bplus) xn(j) = xj
544 go to 300
545
546 * Card ignored.
547
548 290 ignord = ignord + 1
549 if (iprint .gt. 0 .and. ignord .le. mer) then
550 write(iprint, 2010) ncard, key, name1c,name2c,name1r,name2r
551 end if
552 300 continue
553
554 310 ignord = ignord + notfnd
555 if (iprint .gt. 0) write(iprint, 2050) ncard, ignord, nbs, ns
556 if (isumm .gt. 0) write(isumm , 2050) ncard, ignord, nbs, ns
557 if (insrt .ne. iread) rewind insrt
558 return
559
560 1000 format(14x, 2a4, 2x, 3a4)
561 1005 format(2a4)
562 1020 format(3a4, 2x, 2a4, 2x, e12.5)
563 1999 format(/ ' INSERT file to be input from file', i4)
564 2000 format(/ ' NAME', 10x, 2a4, 2x, 3a4)
565 2010 format(' XXX Line', i6, ' ignored:', 8x, 3a4, 2x, 2a4)
566 2050 format(/ ' No. of lines read ', i6, ' Lines ignored', i6
567 $ / ' No. of basics specified', i6, ' Superbasics ', i6)
568
569 * end of m4inst
570 end
571
572 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
573
574 subroutine m4load( m, n, nb, ns,
575 $ hs, bl, bu, xn, name1, name2 )
576
577 implicit double precision (a-h,o-z)
578 integer*4 hs(nb)
579 integer name1(nb), name2(nb)
580 double precision bl(nb), bu(nb)
581 double precision xn(nb)
582
583 * ------------------------------------------------------------------
584 * m4load inputs a load file, which may contain a full or partial
585 * list of row and column names and their states and values.
586 * Valid keys are BS, LL, UL, SB.
587 * ------------------------------------------------------------------
588
589 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
590 common /m1file/ iread,iprint,isumm
591 common /m2file/ iback,idump,iload,imps,inewb,insrt,
592 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
593 common /m3mps3/ aijtol,bstruc(2),mlst,mer,
594 $ aijmin,aijmax,na0,line,ier(20)
595 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
596
597 intrinsic abs
598 integer id(5)
599 character*4 key
600 character*4 lbs , lll , lul ,
601 $ lsb , lend
602 data lbs /' BS '/, lll /' LL '/, lul /' UL '/,
603 $ lsb /' SB '/, lend/'ENDA'/
604
605 bplus = 0.9*plinfy
606 if (iprint .gt. 0) write(iprint, 1999) iload
607 if (isumm .gt. 0) write(isumm , 1999) iload
608 read (iload , 1000) id
609 if (iprint .gt. 0) write(iprint, 2000) id
610 l1 = n + 1
611 ignord = 0
612 nbs = 0
613 ns = 0
614 notfnd = 0
615 ncard = 0
616 jmark = 1
617 ndum = n + 100000
618
619 * Read names until ENDATA is found.
620
621 do 300 nloop = 1, ndum
622 read (iload, 1020) key, id1, id2, xj
623 if (key .eq. lend) go to 310
624
625 ncard = nloop
626 call m4name( nb, name1, name2, id1, id2,
627 $ ncard, notfnd, mer, 1, nb, jmark, j )
628 if (j .le. 0) go to 300
629
630 * The name id1-id2 belongs to the j-th variable.
631
632 if (hs(j) .gt. 1) go to 290
633 if (j .eq.jobj) go to 90
634 if (key .eq. lbs) go to 90
635 if (key .eq. lll) go to 100
636 if (key .eq. lul) go to 150
637 if (key .eq. lsb) go to 200
638 go to 290
639
640 * Make basic.
641
642 90 nbs = nbs + 1
643 hs(j) = 3
644 go to 250
645
646 * LO or UP.
647
648 100 hs(j) = 0
649 go to 250
650
651 150 hs(j) = 1
652 go to 250
653
654 * Make superbasic.
655
656 200 ns = ns + 1
657 hs(j) = 2
658
659 * Save x values.
660
661 250 if (abs(xj) .lt. bplus) xn(j) = xj
662 go to 300
663
664 * Card ignored.
665
666 290 ignord = ignord + 1
667 if (ignord .le. mer) then
668 if (iprint .gt. 0) write(iprint, 2010) ncard, key, id1, id2
669 end if
670 300 continue
671
672 310 ignord = ignord + notfnd
673 if (iprint .gt. 0) write(iprint, 2050) ncard, ignord, nbs, ns
674 if (isumm .gt. 0) write(isumm , 2050) ncard, ignord, nbs, ns
675
676 * Make sure the linear objective is basic.
677
678 if (iobj .gt. 0 .and. hs(jobj) .ne. 3) then
679 hs(jobj) = 3
680
681 * Swap obj with last basic variable.
682
683 do 850 j = nb, 1, -1
684 if (hs(j) .eq. 3) go to 860
685 850 continue
686
687 860 hs(j) = 0
688 end if
689
690 if (iload .ne. iread) rewind iload
691 return
692
693 1000 format(14x, 2a4, 2x, 3a4)
694 1005 format(2a4)
695 1020 format(3a4, 12x, e12.5)
696 1999 format(/ ' LOAD file to be input from file', i4)
697 2000 format(/ ' NAME', 10x, 2a4, 2x, 3a4)
698 2010 format(' XXX Line', i6, ' ignored:', 8x, 3a4, 2x, 2a4)
699 2050 format(/ ' No. of lines read ', i6, ' Lines ignored', i6
700 $ / ' No. of basics specified', i6, ' Superbasics ', i6)
701
702 * end of m4load
703 end
704
705 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
706
707 subroutine m4oldb( m, n, nb, ns,
708 $ hs, bl, bu, xn )
709
710 implicit double precision (a-h,o-z)
711 integer*4 hs(nb)
712 double precision bl(nb), bu(nb)
713 double precision xn(nb)
714
715 * ------------------------------------------------------------------
716 * m4oldb inputs a compact basis file from file ioldb.
717 * ------------------------------------------------------------------
718
719 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
720 common /m1file/ iread,iprint,isumm
721 common /m2file/ iback,idump,iload,imps,inewb,insrt,
722 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
723 common /m5log1/ idebug,ierr,lprint
724
725 character*80 id
726
727 bplus = 0.9*plinfy
728 if (iprint .gt. 0) write(iprint, 1999) ioldb
729 if (isumm .gt. 0) write(isumm , 1999) ioldb
730 read (ioldb , 1000) id
731 if (iprint .gt. 0) then
732 write(iprint, 2000) id
733 end if
734
735 read (ioldb , 1005) id(1:52), newm, newn, ns
736 if (iprint .gt. 0) then
737 write(iprint, 2005) id(1:52), newm, newn, ns
738 end if
739
740 if (newm .ne. m .or. newn .ne. n) go to 900
741 read (ioldb , 1010) hs
742
743 * Set values for nonbasic variables.
744
745 do 200 j = 1, nb
746 js = hs(j)
747 if (js .le. 1) then
748 if (js .eq. 0) xj = bl(j)
749 if (js .eq. 1) xj = bu(j)
750 if (abs(xj) .lt. bplus) xn(j) = xj
751 end if
752 200 continue
753
754 * Load superbasics.
755
756 ns = 0
757 ndummy = m + n + 10000
758
759 do 300 idummy = 1, ndummy
760 read(ioldb, 1020, end=310) j, xj
761 if (j .le. 0) go to 310
762 if (j .le. nb) then
763 xn(j) = xj
764 if (hs(j) .eq. 2) ns = ns + 1
765 end if
766 300 continue
767
768 310 if (ns .gt. 0) then
769 if (iprint .gt. 0) write(iprint, 2010) ns
770 if (isumm .gt. 0) write(isumm , 2010) ns
771 end if
772 go to 990
773
774 * Error exits.
775
776 900 call m1page( 1 )
777 if (iprint .gt. 0) write(iprint, 3000)
778 if (isumm .gt. 0) write(isumm , 3000)
779 ierr = 30
780
781 990 if (ioldb .ne. iread) rewind ioldb
782 return
783
784 1000 format(a80)
785 1005 format(a52, 2x, i7, 3x, i7, 4x, i5)
786 1010 format(80i1)
787 1020 format(i8, e24.14)
788 1999 format(/ ' OLD BASIS file to be input from file', i4)
789 2000 format(1x, a80)
790 2005 format(1x, a52, 'M=', i7, ' N=', i7, ' SB=', i5)
791 2010 format(' No. of superbasics loaded', i7)
792 3000 format(' EXIT -- the basis file dimensions do not match',
793 $ ' this problem')
794
795 * end of m4oldb
796 end
797
798 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
799
800 subroutine m4savb( mode, m, mbs, n, nb, nn, nname, nscl, msoln,ns,
801 $ ne, nka, a, ha, ka,
802 $ hs, kb, ascale, bl, bu,
803 $ name1, name2, pi, rc, xn, y, z, nwcore )
804
805 implicit double precision (a-h,o-z)
806 integer*4 ha(ne), hs(nb)
807 integer ka(nka), kb(mbs), name1(nname), name2(nname)
808 double precision a(ne), ascale(nscl), bl(nb), bu(nb)
809 double precision pi(m), rc(nb), xn(nb), y(m), z(nwcore)
810
811 * ------------------------------------------------------------------
812 * m4savb saves basis files and/or prints the solution.
813 *
814 * If mode = 1, the problem is first unscaled, then from 0 to 4 files
815 * are saved (PUNCH file, DUMP file, SOLUTION file, REPORT file,
816 * in that order).
817 * A new BASIS file, if any, will already have been saved by m5solv.
818 * A call with mode = 1 must precede a call with mode = 2.
819 *
820 * If mode = 2, the solution is printed under the control of msoln
821 * (which is set by the Solution keyword in the SPECS file).
822 *
823 * 18 Nov 1991: Scaled pinorm saved for use in m4soln.
824 * 25 Nov 1991: rc added as parameter to return reduced costs.
825 * 31 Jan 1991: Call m4rc to get the reduced costs.
826 * 18 Dec 1992: Maximum primal and dual infeasibilities computed
827 * and printed here.
828 * ------------------------------------------------------------------
829
830 logical alone, AMPL, GAMS, MINT, page1, page2
831 common /m1env / alone, AMPL, GAMS, MINT, page1, page2
832 common /m1file/ iread,iprint,isumm
833 common /m2file/ iback,idump,iload,imps,inewb,insrt,
834 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
835 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
836 common /m3scal/ sclobj,scltol,lscale
837 common /m5inf / prinf, duinf, jprinf, jduinf
838 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
839 common /m5log1/ idebug,ierr,lprint
840 common /m5step/ featol, tolx0,tolinc,kdegen,ndegen,
841 $ itnfix, nfix(2)
842 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
843 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
844 common /m7loc / lgobj ,lgobj2
845 common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg
846 common /m8len / njac ,nncon ,nncon0,nnjac
847 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
848 $ lblslk,lbuslk,lxlam ,lrhs ,
849 $ lgcon ,lgcon2,lxdif ,lxold
850 common /m8save/ vimax ,virel ,maxvi ,majits,minits,nssave
851
852 intrinsic max
853 character*12 istate
854 logical feasbl, prnt
855 save pnorm1, pnorm2
856 parameter ( one = 1.0d+0 )
857
858 feasbl = ninf .eq. 0
859 ms = m + ns
860 k = ierr + 1
861 call m4stat( k, istate )
862
863 if (mode .eq. 1) then
864 * ---------------------------------------------------------------
865 * mode = 1.
866 * Compute rc and unscale everything.
867 * Then save basis files.
868 * ---------------------------------------------------------------
869
870 * Compute reduced costs rc(*) for all columns and rows.
871 * Find the maximum primal and dual infeasibilities.
872
873 call m4rc ( feasbl, featol, minimz,
874 $ m, n, nb, nnobj, nnobj0,
875 $ ne, nka, a, ha, ka,
876 $ hs, bl, bu, z(lgobj), pi, rc, xn )
877 call m4infs( nb, jobj, bl, bu, rc, xn )
878
879 * Unscale a, bl, bu, pi, xn, rc, fcon, gobj and xnorm, pinorm.
880 * (m4soln uses scaled pinorm, so save it.)
881
882 pnorm1 = pinorm
883 xnorm1 = xnorm
884 jpinf1 = jprinf
885 jdinf1 = jduinf
886 prinf1 = prinf
887 duinf1 = duinf
888
889 if (lscale .gt. 0) then
890 call m2scla( 2, m, n, nb, ne, nka,
891 $ ha, ka, a, ascale, bl, bu, pi, xn )
892
893 call dddiv ( nb, ascale, 1, rc, 1 )
894
895 if (lscale .eq. 2) then
896 if (nncon .gt. 0)
897 $ call ddscl ( nncon, ascale(n+1), 1, z(lfcon), 1 )
898
899 if (nnobj .gt. 0)
900 $ call dddiv ( nnobj, ascale , 1, z(lgobj), 1 )
901 end if
902
903 * Call m5setp to redefine pinorm. y is not used.
904
905 xnorm = dnorm1( nb, xn, 1 )
906 call m5setp( 3, m, y, pi, z, nwcore )
907 call m4infs( nb, jobj, bl, bu, rc, xn )
908 end if
909 pnorm2 = pinorm
910
911 * ---------------------------------------------------------------
912 * Print various scaled and unscaled norms.
913 * ---------------------------------------------------------------
914 if (lscale .gt. 0) then
915 if (iprint .gt. 0) write(iprint, 1010) xnorm1, pnorm1
916 if (isumm .gt. 0) write(isumm , 1010) xnorm1, pnorm1
917 end if
918 if (iprint .gt. 0) write(iprint, 1020) xnorm , pinorm
919 if (isumm .gt. 0) write(isumm , 1020) xnorm , pinorm
920 if (lscale .gt. 0) then
921 if (iprint .gt. 0) write(iprint, 1030) jpinf1, prinf1,
922 $ jdinf1, duinf1
923 if (isumm .gt. 0) write(isumm , 1030) jpinf1, prinf1,
924 $ jdinf1, duinf1
925 end if
926 if (iprint .gt. 0) write(iprint, 1040) jprinf, prinf,
927 $ jduinf, duinf
928 if (isumm .gt. 0) write(isumm , 1040) jprinf, prinf,
929 $ jduinf, duinf
930
931 * Change the sign of pi and rc if feasible and maximizing.
932
933 if (ninf .eq. 0 .and. minimz .lt. 0) then
934 call dscal ( m , -one, pi, 1 )
935 call dscal ( nb, -one, rc, 1 )
936 end if
937
938 * Compute nonlinear constraint violations.
939
940 if (nncon .gt. 0) then
941 call m8viol( n, nb, nncon,
942 $ ne, nka, a, ha, ka,
943 $ bl, bu, z(lfcon), xn, y, z, nwcore )
944 if (iprint .gt. 0) write(iprint, 1080) vimax, virel
945 if (isumm .gt. 0) write(isumm , 1080) vimax, virel
946 end if
947
948 * ---------------------------------------------------------------
949 * Output PUNCH, DUMP, SOLUTION and/or REPORT files.
950 * ---------------------------------------------------------------
951 if (ipnch .gt. 0)
952 $ call m4pnch( ipnch, n, nb, hs, bl, bu, xn,
953 $ name1, name2 )
954
955 if (idump .gt. 0)
956 $ call m4dump( idump, n, nb, hs, bl, bu, xn,
957 $ name1, name2 )
958
959 pinorm = pnorm1
960 if (isoln .gt. 0)
961 $ call m4soln( .true., m, n, nb, nname, nscl,
962 $ nn, nnobj, nnobj0, ns,
963 $ ne, nka, a, ha, ka,
964 $ hs, ascale, bl, bu,
965 $ z(lgobj), pi, rc, xn, y,
966 $ name1, name2, istate, z, nwcore )
967
968 if (ireprt .gt. 0)
969 $ call m4rept( .true., m, n, nb, nname, nscl,
970 $ nn, nnobj, nnobj0, ns,
971 $ ne, nka, a, ha, ka,
972 $ hs, ascale, bl, bu,
973 $ z(lgobj), pi, rc, xn, y,
974 $ name1, name2, istate, z, nwcore )
975 pinorm = pnorm2
976 else
977 * ---------------------------------------------------------------
978 * mode = 2. Print solution if requested.
979 *
980 * msoln = 0 means no
981 * = 1 means if optimal, infeasible or unbounded
982 * = 2 means yes
983 * = 3 means if error condition
984 * ---------------------------------------------------------------
985 prnt = iprint .gt. 0 .and. msoln .gt. 0
986 if ((msoln .eq. 1 .and. ierr .gt. 2) .or.
987 $ (msoln .eq. 3 .and. ierr .le. 2)) prnt = .false.
988 if ( prnt ) then
989 pinorm = pnorm1
990 call m4soln( .false., m, n, nb, nname, nscl,
991 $ nn, nnobj, nnobj0, ns,
992 $ ne, nka, a, ha, ka,
993 $ hs, ascale, bl, bu,
994 $ z(lgobj), pi, rc, xn, y,
995 $ name1, name2, istate, z, nwcore )
996 pinorm = pnorm2
997 if (isumm .gt. 0) write(isumm, 1200) iprint
998
999 else if (.not. (GAMS .or. AMPL)) then
1000 if (isumm .gt. 0) write(isumm, 1300)
1001 end if
1002 end if
1003
1004 return
1005
1006 1010 format( ' Norm of x (scaled) ', 1p, e16.3,
1007 $ 2x, ' Norm of pi (scaled) ', e16.3)
1008 1020 format( ' Norm of x ', 1p, e27.3,
1009 $ 2x, ' Norm of pi', e27.3)
1010 1030 format( ' Primal inf (scaled)', i7, 1p, e11.1,
1011 $ 2x, ' Dual inf (scaled)', i7, e11.1)
1012 1040 format( ' Primal infeas ', i7, 1p, e11.1,
1013 $ 2x, ' Dual infeas ', i7, e11.1)
1014 1080 format( ' Constraint violation ', 1p, e16.3,
1015 $ 2x, ' Normalized ', e16.3)
1016 1100 format(2a4)
1017 1200 format(/ ' Solution printed on file', i4)
1018 1300 format(/ ' Solution not printed')
1019
1020 * end of m4savb
1021 end
1022
1023 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1024
1025 subroutine m4dump( idump, n, nb, hs, bl, bu, xn, name1, name2 )
1026
1027 implicit double precision (a-h,o-z)
1028 integer*4 hs(nb)
1029 integer name1(nb), name2(nb)
1030 double precision bl(nb), bu(nb)
1031 double precision xn(nb)
1032
1033 * ------------------------------------------------------------------
1034 * m4dump outputs basis names in a format compatible with m4load.
1035 * This file is normally easier to modify than a punch file.
1036 * ------------------------------------------------------------------
1037
1038 common /m1file/ iread,iprint,isumm
1039 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1040
1041 character*4 key(4)
1042 data key /' LL ', ' UL ', ' SB ', ' BS '/
1043
1044 write(idump, 2000) name
1045
1046 do 500 j = 1, nb
1047 k = hs(j) + 1
1048 write(idump, 2100) key(k), name1(j), name2(j), xn(j)
1049 500 continue
1050
1051 write(idump , 2200)
1052 if (iprint .gt. 0) write(iprint, 3000) idump
1053 if (isumm .gt. 0) write(isumm , 3000) idump
1054 if (idump .ne. iprint) rewind idump
1055 return
1056
1057 2000 format('NAME', 10x, 2a4, 2x, ' DUMP/LOAD')
1058 2100 format(3a4, 12x, 1p, e12.5)
1059 2200 format('ENDATA')
1060 3000 format(/ ' DUMP file saved on file', i4)
1061
1062 * end of m4dump
1063 end
1064
1065 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1066
1067 subroutine m4newb( mode, inewb, m, n, nb, nn, ns, ms, nscl, fsub,
1068 $ kb, hs, ascale, bl, bu, x, xn, istate )
1069
1070 implicit double precision (a-h,o-z)
1071 integer*4 hs(nb)
1072 integer kb(ms)
1073 double precision ascale(nscl), bl(nb), bu(nb)
1074 double precision x(ms), xn(nb)
1075 character*12 istate
1076
1077 * ------------------------------------------------------------------
1078 * m4newb saves a compact basis on file inewb. Called from m5solv.
1079 * If mode = 1, the save is a periodic one due to the save frequency.
1080 * If mode = 2, m5solv has just finished the current problem.
1081 * (Hence, mode 2 calls occur at the end of every cycle.)
1082 * ------------------------------------------------------------------
1083
1084 common /m1file/ iread,iprint,isumm
1085 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1086 common /m3scal/ sclobj,scltol,lscale
1087 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1088 logical prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
1089 common /m5log4/ prnt0 ,prnt1 ,summ0 ,summ1 ,newhed
1090 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1091 common /m8len / njac ,nncon ,nncon0,nnjac
1092
1093 logical scaled
1094
1095 scaled = lscale .gt. 0
1096 obj = sinf
1097 if (ninf .eq. 0) obj = minimz * fsub
1098
1099 * Output header cards and the state vector.
1100
1101 write(inewb, 1000) name, itn, istate, ninf, obj
1102 write(inewb, 1005) mobj, mrhs, mrng, mbnd, m, n, ns
1103 write(inewb, 1010) hs
1104
1105 * Output the superbasic variables.
1106
1107 do 50 k = m + 1, ms
1108 j = kb(k)
1109 xj = x(k)
1110 if (scaled) xj = xj * ascale(j)
1111 write(inewb, 1020) j, xj
1112 50 continue
1113
1114 * If there are nonlinear constraints,
1115 * output the values of all other (non-sb) nonlinear variables.
1116
1117 do 100 j = 1, nnjac
1118 if (hs(j) .ne. 2) then
1119 xj = xn(j)
1120 if (scaled) xj = xj * ascale(j)
1121 write(inewb, 1020) j, xj
1122 end if
1123 100 continue
1124
1125 * Output nonbasic variables that are not at a bound.
1126 * Ignore ones that the EXPAND anti-cycling procedure
1127 * has put slightly outside their bound.
1128
1129 do 250 j = nnjac + 1, nb
1130 if (hs(j) .le. 1 ) then
1131 xj = xn(j)
1132 if (xj .le. bl(j)) go to 250
1133 if (xj .ge. bu(j)) go to 250
1134 if (scaled) xj = xj * ascale(j)
1135 write(inewb, 1020) j, xj
1136 end if
1137 250 continue
1138
1139 * Terminate the list with a zero.
1140
1141 j = 0
1142 write(inewb, 1020) j
1143 if (inewb .ne. iprint) rewind inewb
1144 if (iprint .gt. 0) write(iprint, 1030) inewb, itn
1145 if (isumm .gt. 0) write(isumm , 1030) inewb, itn
1146 return
1147
1148 1000 format(2a4, ' ITN', i8, 4x, a12, ' NINF', i7,
1149 $ ' OBJ', 1p, e21.12)
1150 1005 format('OBJ=', 2a4, ' RHS=', 2a4, ' RNG=', 2a4, ' BND=', 2a4,
1151 $ ' M=', i7, ' N=', i7, ' SB=', i5)
1152 1010 format(80i1)
1153 1020 format(i8, 1p, e24.14)
1154 1030 format(/ ' NEW BASIS file saved on file', i4, ' itn =', i7)
1155
1156 * end of m4newb
1157 end
1158
1159 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1160
1161 subroutine m4pnch( ipnch, n, nb, hs, bl, bu, xn, name1, name2 )
1162
1163 implicit double precision (a-h,o-z)
1164 integer*4 hs(nb)
1165 integer name1(nb), name2(nb)
1166 double precision bl(nb), bu(nb)
1167 double precision xn(nb)
1168
1169 * ------------------------------------------------------------------
1170 * m4pnch outputs a PUNCH file (list of basis names, states and
1171 * values) in a format that is compatible with MPS/360.
1172 * ------------------------------------------------------------------
1173
1174 common /m1file/ iread,iprint,isumm
1175 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1176
1177 parameter ( zero = 0.0d+0 )
1178 character*4 key(5), ibl
1179 data key /' LL ', ' UL ', ' SB ', ' XL ', ' XU '/
1180 data ibl /' '/
1181
1182 write(ipnch, 2000) name
1183 irow = n
1184
1185 do 500 j = 1, n
1186 id1 = name1(j)
1187 id2 = name2(j)
1188 k = hs(j)
1189
1190 if (k .eq. 3) then
1191
1192 * Basics -- find the next row that isn't basic.
1193
1194 300 irow = irow + 1
1195 if (irow .le. nb) then
1196 k = hs(irow)
1197 if (k .eq. 3) go to 300
1198
1199 if (k .eq. 2) k = 0
1200 write(ipnch, 2100) key(k+4), id1, id2,
1201 $ name1(irow), name2(irow), xn(j)
1202 end if
1203 else
1204
1205 * Skip nonbasic variables with zero lower bounds.
1206
1207 if (k .le. 1) then
1208 if (bl(j) .eq. zero .and. xn(j) .eq. zero) go to 500
1209 end if
1210 write(ipnch, 2100) key(k+1), id1, id2, ibl, ibl, xn(j)
1211 end if
1212 500 continue
1213
1214 * Output superbasic slacks.
1215
1216 do 700 j = n + 1, nb
1217 if (hs(j) .eq. 2)
1218 $ write(ipnch, 2100) key(3), name1(j), name2(j), ibl, ibl, xn(j)
1219 700 continue
1220
1221 write(ipnch , 2200)
1222 if (iprint .gt. 0) write(iprint, 3000) ipnch
1223 if (isumm .gt. 0) write(isumm , 3000) ipnch
1224 if (ipnch .ne. iprint) rewind ipnch
1225 return
1226
1227 2000 format('NAME', 10x, 2a4, 2x, 'PUNCH/INSERT')
1228 2100 format(3a4, 2x, 2a4, 2x, 1p, e12.5)
1229 2200 format('ENDATA')
1230 3000 format(/ ' PUNCH file saved on file', i4)
1231
1232 * end of m4pnch
1233 end
1234
1235 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1236
1237 subroutine m4rc ( feasbl, featol, minimz,
1238 $ m, n, nb, nnobj, nnobj0,
1239 $ ne, nka, a, ha, ka,
1240 $ hs, bl, bu, gobj, pi, rc, xn )
1241
1242 implicit double precision (a-h,o-z)
1243 logical feasbl
1244 integer*4 ha(ne), hs(nb)
1245 integer ka(nka)
1246 double precision a(ne), bl(nb), bu(nb)
1247 double precision gobj(nnobj0), pi(m), rc(nb), xn(nb)
1248
1249 * ------------------------------------------------------------------
1250 * m4rc computes reduced costs rc(*) for all columns of ( A I ).
1251 * If xn is feasible, the true nonlinear objective gradient gobj(*)
1252 * is used (not the gradient of the augmented Lagrangian).
1253 * Otherwise, the Phase-1 objective is included.
1254 *
1255 * m4rc is called by m4savb BEFORE unscaling.
1256 * External values of hs(*) are used (0, 1, 2, 3),
1257 * but internal would be ok too since we only test for > 1.
1258 *
1259 * 31 Jan 1992: First version.
1260 * ------------------------------------------------------------------
1261
1262 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1263
1264 do 300 j = 1, n
1265 dj = zero
1266 do 200 k = ka(j), ka(j+1) - 1
1267 i = ha(k)
1268 dj = dj + pi(i) * a(k)
1269 200 continue
1270 rc(j) = - dj
1271 300 continue
1272
1273 do 320 i = 1, m
1274 rc(n+i) = - pi(i)
1275 320 continue
1276
1277 if (feasbl .and. nnobj .gt. 0) then
1278
1279 * Include the nonlinear objective gradient.
1280
1281 sgnobj = minimz
1282 call daxpy ( nnobj, sgnobj, gobj, 1, rc, 1 )
1283 else
1284
1285 * Include the Phase 1 objective.
1286 * Only basics and superbasics can be infeasible.
1287
1288 do 500 j = 1, nb
1289 if (hs(j) .gt. 1) then
1290 d1 = bl(j) - xn(j)
1291 d2 = xn(j) - bu(j)
1292 if (d1 .gt. featol) rc(j) = rc(j) - one
1293 if (d2 .gt. featol) rc(j) = rc(j) + one
1294 end if
1295 500 continue
1296 end if
1297
1298 * end of m4rc
1299 end
1300
1301 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1302
1303 subroutine m4infs( nb, jobj, bl, bu, rc, xn )
1304
1305 implicit double precision (a-h,o-z)
1306 double precision bl(nb), bu(nb), rc(nb), xn(nb)
1307
1308 * ------------------------------------------------------------------
1309 * m4infs computes the maximum primal and dual infeasibilities,
1310 * using bl, bu, rc and xn.
1311 *
1312 * m4infs is called by m4savb before and after unscaling.
1313 *
1314 * 18 Dec 1992: First version.
1315 * prinf , duinf return the max primal and dual infeas.
1316 * jprinf, jduinf return the corresponding column nos.
1317 * ------------------------------------------------------------------
1318
1319 common /m5inf / prinf, duinf, jprinf, jduinf
1320 parameter ( zero = 0.0d+0 )
1321
1322 jprinf = 0
1323 jduinf = 0
1324 prinf = zero
1325 duinf = zero
1326
1327 do 500 j = 1, nb
1328 d1 = bl(j) - xn(j)
1329 d2 = xn(j) - bu(j)
1330 if (prinf .lt. d1) then
1331 prinf = d1
1332 jprinf = j
1333 end if
1334 if (prinf .lt. d2) then
1335 prinf = d2
1336 jprinf = j
1337 end if
1338
1339 if ( j .eq. jobj ) go to 500
1340 if (bl(j) .eq. bu(j)) go to 500
1341
1342 dj = rc(j)
1343 if (d1 .ge. zero) then
1344 dj = - dj
1345 else if (d2 .le. zero) then
1346 * dj = + dj
1347 else
1348 dj = abs( dj )
1349 end if
1350
1351 if (duinf .lt. dj) then
1352 duinf = dj
1353 jduinf = j
1354 end if
1355 500 continue
1356
1357 * end of m4infs
1358 end
1359
1360 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1361
1362 subroutine m4rept( ondisk, m, n, nb, nname, nscl,
1363 $ nn, nnobj, nnobj0, ns,
1364 $ ne, nka, a, ha, ka,
1365 $ hs, ascale, bl, bu, gobj, pi, rc, xn, y,
1366 $ name1, name2, istate, z, nwcore )
1367
1368 implicit double precision (a-h,o-z)
1369 logical ondisk
1370 integer*4 ha(ne), hs(nb)
1371 integer ka(nka), name1(nname), name2(nname)
1372 character*12 istate
1373 double precision a(ne), ascale(nscl), bl(nb), bu(nb)
1374 double precision gobj(nnobj0), pi(m), rc(nb), xn(nb), y(m),
1375 $ z(nwcore)
1376
1377 * ------------------------------------------------------------------
1378 * m4rept has the same parameter list as m4soln, the routine that
1379 * prints the solution. It will be called if the SPECS file
1380 * specifies REPORT file n for some positive value of n.
1381 *
1382 * pi contains the unscaled dual solution.
1383 * xn contains the unscaled primal solution. There are n + m = nb
1384 * values (n structural variables and m slacks, in that order).
1385 * y contains the true slack values for nonlinear constraints
1386 * in its first nncon components (computed by m8viol).
1387 *
1388 * This version of m4rept does nothing. Added for PILOT, Oct 1985.
1389 * 31 Oct 1991: Name changed from "report" to "m4rept".
1390 * Parameters altered to allow for MPS or generic names.
1391 * ------------------------------------------------------------------
1392
1393 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1394 common /m1file/ iread,iprint,isumm
1395 common /m2file/ iback,idump,iload,imps,inewb,insrt,
1396 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
1397 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1398 common /m3scal/ sclobj,scltol,lscale
1399 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1400 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1401 common /cycle2/ objtru,suminf,numinf
1402 * ------------------------------------------------------------------
1403
1404 if (iprint .gt. 0) write(iprint, 1000)
1405 if (isumm .gt. 0) write(isumm , 1000)
1406 return
1407
1408 1000 format(/ ' XXX Report file requested. m4rept does nothing.')
1409
1410 * end of m4rept
1411 end
1412
1413 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1414
1415 subroutine m4soln( ondisk, m, n, nb, nname, nscl,
1416 $ nn, nnobj, nnobj0, ns,
1417 $ ne, nka, a, ha, ka,
1418 $ hs, ascale, bl, bu, gobj, pi, rc, xn, y,
1419 $ name1, name2, istate, z, nwcore )
1420
1421 implicit double precision (a-h,o-z)
1422 logical ondisk
1423 integer*4 ha(ne), hs(nb)
1424 integer ka(nka), name1(nname), name2(nname)
1425 character*12 istate
1426 double precision a(ne), ascale(nscl), bl(nb), bu(nb)
1427 double precision gobj(nnobj0), pi(m), rc(nb), xn(nb), y(m),
1428 $ z(nwcore)
1429
1430 * ------------------------------------------------------------------
1431 * m4soln is the standard output routine for printing the solution.
1432 *
1433 * On entry,
1434 * pi contains the dual solution.
1435 * xn contains the primal solution. There are n + m = nb values
1436 * (n structural variables and m slacks, in that order).
1437 * rc contains reduced costs for all variables:
1438 * rc(1:n) = gobj + c - A'pi for structurals,
1439 * rc(n+1:nb) = -pi for slacks.
1440 * pi and rc corresond to the Phase-1 objective
1441 * if the solution is infeasible.
1442 * y contains the true slack values for nonlinear constraints
1443 * in its first nncon components (computed by m8viol).
1444 *
1445 * All quantities a, bl, bu, pi, rc, xn, y, fcon, gobj are unscaled
1446 * and adjusted in sign if maximizing. (fcon is not used here.)
1447 *
1448 * If ondisk is true, the solution is output to the solution file.
1449 * Otherwise, it is output to the printer.
1450 *
1451 * 31 Jan 1991: rc is now an input parameter. It is set in m4savb.
1452 * ------------------------------------------------------------------
1453
1454 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1455 common /m1file/ iread,iprint,isumm
1456 common /m2file/ iback,idump,iload,imps,inewb,insrt,
1457 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
1458 common /m3mps4/ name(2),mobj(2),mrhs(2),mrng(2),mbnd(2),minmax
1459 common /m3scal/ sclobj,scltol,lscale
1460 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1461 common /m5lp1 / itn,itnlim,nphs,kmodlu,kmodpi
1462 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1463 common /m7tols/ xtol(2),ftol(2),gtol(2),pinorm,rgnorm,tolrg
1464 common /m8len / njac ,nncon ,nncon0,nnjac
1465 common /cycle2/ objtru,suminf,numinf
1466
1467 intrinsic abs
1468 logical feasbl, infsbl, maximz, scaled
1469 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1470 * ------------------------------------------------------------------
1471
1472 bplus = 0.1*plinfy
1473 scale = one
1474 feasbl = ninf .eq. 0
1475 infsbl = .not. feasbl
1476 maximz = minimz .lt. 0
1477 scaled = lscale .gt. 0
1478 lpr = iprint
1479 if (ondisk) lpr = isoln
1480
1481 call m1page( 1 )
1482 if (infsbl) write(lpr, 1000) name, ninf, sinf
1483 if (feasbl) write(lpr, 1002) name, objtru
1484 write(lpr, 1004) istate, itn, ns
1485 write(lpr, 1005) mobj, minmax, mrhs, mrng, mbnd
1486 write(lpr, 1010)
1487 ** tolfea = 0.1 * tolx
1488 ** tolopt = 0.1 * toldj(3) * pinorm
1489 ** 05 Oct 1991: Might as well flag according to the tolerances
1490 * actually used.
1491 tolfea = tolx
1492 tolopt = toldj(3) * pinorm
1493
1494 * ------------------------------------------------------------------
1495 * Output the ROWS section.
1496 * ------------------------------------------------------------------
1497 do 300 iloop = 1, m
1498 i = iloop
1499 j = n + i
1500 if (scaled) scale = ascale(j)
1501 js = hs(j)
1502 b1 = bl(j)
1503 b2 = bu(j)
1504 xj = xn(j)
1505 py = pi(i)
1506 dj = rc(j)
1507
1508 * Define row and slack activities.
1509
1510 if (i .le. nncon) xj = y(i)
1511 row = - xj
1512 d1 = b1 - xj
1513 d2 = xj - b2
1514 slk = - d1
1515 if (abs( d1 ) .gt. abs( d2 )) slk = d2
1516 if (abs( slk ) .ge. bplus ) slk = - row
1517 d1 = d1 / scale
1518 d2 = d2 / scale
1519 djtest = dj * scale
1520 if (feasbl) then
1521 if ( maximz ) djtest = - djtest
1522 end if
1523
1524 * Change slacks into rows.
1525
1526 if (js .le. 1) js = 1 - js
1527 b1 = - b2
1528 b2 = - bl(j)
1529
1530 call m4id ( j, m, n, nb, nname, name1, name2, id1, id2 )
1531 call m4solp( ondisk, bplus, tolfea, tolopt,
1532 $ js, d1, d2, djtest,
1533 $ j , id1, id2, row, slk, b1, b2, py, i )
1534 300 continue
1535
1536 * ------------------------------------------------------------------
1537 * Output the COLUMNS section.
1538 * ------------------------------------------------------------------
1539 call m1page( 1 )
1540 write(lpr, 1020)
1541
1542 do 400 jloop = 1, n
1543 j = jloop
1544 if (scaled) scale = ascale(j)
1545 js = hs(j)
1546 b1 = bl(j)
1547 b2 = bu(j)
1548 xj = xn(j)
1549 cj = zero
1550 dj = rc(j)
1551
1552 do 320 k = ka(j), ka(j+1) - 1
1553 ir = ha(k)
1554 if (ir .eq. iobj) cj = a(k)
1555 320 continue
1556
1557 d1 = (b1 - xj) / scale
1558 d2 = (xj - b2) / scale
1559 djtest = - dj * scale
1560 if (feasbl) then
1561 if (j .le. nnobj) cj = cj + gobj(j)
1562 if ( maximz ) djtest = - djtest
1563 end if
1564
1565 call m4id ( j, m, n, nb, nname, name1, name2, id1, id2 )
1566 call m4solp( ondisk, bplus, tolfea, tolopt,
1567 $ js, d1, d2, djtest,
1568 $ j , id1, id2, xj, cj, b1, b2, dj, m+j )
1569 400 continue
1570
1571 if (ondisk) then
1572 if (isoln .ne. iprint) rewind isoln
1573 if (iprint .gt. 0) write(iprint, 1400) isoln
1574 if (isumm .gt. 0) write(isumm , 1400) isoln
1575 end if
1576 return
1577
1578 1000 format(' NAME', 11x, 2a4, 13x,
1579 $ ' INFEASIBILITIES', i7, 1p, e16.4)
1580 1002 format(' NAME', 11x, 2a4, 13x,
1581 $ ' OBJECTIVE VALUE', 1p, e23.10)
1582 1004 format(/ ' STATUS', 9x, a12, 9x,
1583 $ ' ITERATION', i7, ' SUPERBASICS', i7)
1584 1005 format(/
1585 $ ' OBJECTIVE', 6x, 2a4, ' (', a3, ')' /
1586 $ ' RHS ', 6x, 2a4 /
1587 $ ' RANGES ', 6x, 2a4 /
1588 $ ' BOUNDS ', 6x, 2a4)
1589 1010 format(/ ' SECTION 1 - ROWS' //
1590 $ ' NUMBER ...ROW.. STATE ...ACTIVITY... SLACK ACTIVITY',
1591 $ ' ..LOWER LIMIT. ..UPPER LIMIT. .DUAL ACTIVITY ..I' /)
1592 1020 format( ' SECTION 2 - COLUMNS' //
1593 $ ' NUMBER .COLUMN. STATE ...ACTIVITY... .OBJ GRADIENT.',
1594 $ ' ..LOWER LIMIT. ..UPPER LIMIT. REDUCED GRADNT M+J' /)
1595 1400 format(/ ' SOLUTION file saved on file', i4)
1596
1597 * end of m4soln
1598 end
1599
1600 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1601
1602 subroutine m4solp( ondisk, bplus, tolfea, tolopt,
1603 $ js, d1, d2, djtest,
1604 $ j , id1, id2, xj, cj, b1, b2, dj, k )
1605
1606 implicit double precision (a-h,o-z)
1607 logical ondisk
1608
1609 * ------------------------------------------------------------------
1610 * m4solp prints one line of the Solution file.
1611 *
1612 * The following conditions are marked by key:
1613 *
1614 * D degenerate basic or superbasic variable.
1615 * I infeasible basic or superbasic variable.
1616 * A alternative optimum (degenerate nonbasic dual).
1617 * N nonoptimal superbasic or nonbasic (infeasible dual).
1618 *
1619 * Tests for these conditions are performed on scaled quantities
1620 * d1, d2, djtest,
1621 * since the correct indication is then more likely to be given.
1622 * On badly scaled problems, the unscaled solution may then appear
1623 * to be flagged incorrectly, but this is just an illusion.
1624 * ------------------------------------------------------------------
1625
1626 common /m1file/ iread,iprint,isumm
1627 common /m2file/ iback,idump,iload,imps,inewb,insrt,
1628 $ ioldb,ipnch,iprob,iscr,isoln,ispecs,ireprt
1629
1630 intrinsic abs
1631 character*4 jstat, jstate(6)
1632 character*1 key, lblank, laltop, ldegen, linfea, lnotop
1633 data jstate /' LL ', ' UL ', 'SBS ',
1634 $ ' BS' , ' EQ' , ' FR '/
1635 data lblank /' '/, laltop /'A'/, ldegen /'D'/,
1636 $ linfea /'I'/, lnotop /'N'/
1637 * ------------------------------------------------------------------
1638
1639 key = lblank
1640 if (js .le. 1) then
1641
1642 * Set key for nonbasic variables.
1643
1644 if (b1 .eq. b2) js = 4
1645 if (- d1 .gt. tolfea .and. - d2 .gt. tolfea) js = 5
1646 if (js .eq. 1 ) djtest = - djtest
1647 if (js .ge. 4 ) djtest = abs( djtest )
1648 if ( abs( djtest ) .le. tolopt) key = laltop
1649 if (js .ne. 4 .and. djtest .gt. tolopt) key = lnotop
1650 else
1651
1652 * Set key for basic and superbasic variables.
1653
1654 if (abs( d1 ) .le. tolfea .or.
1655 $ abs( d2 ) .le. tolfea) key = ldegen
1656 if ( js .eq. 2 .and.
1657 $ abs( djtest ) .gt. tolopt) key = lnotop
1658 if ( d1 .gt. tolfea .or.
1659 $ d2 .gt. tolfea) key = linfea
1660 end if
1661
1662 * Select format for printing.
1663
1664 jstat = jstate(js + 1)
1665 if (ondisk) then
1666 write(isoln, 1000) j,id1,id2,key,jstat,xj,cj,b1,b2,dj,k
1667 else
1668 if (b2 .lt. bplus) then
1669 if (b1 .gt. - bplus) then
1670 write(iprint, 1200) j,id1,id2,key,jstat,xj,cj,b1,b2,dj,k
1671 else
1672 write(iprint, 1300) j,id1,id2,key,jstat,xj,cj, b2,dj,k
1673 end if
1674 else
1675 if (b1 .gt. - bplus) then
1676 write(iprint, 1400) j,id1,id2,key,jstat,xj,cj,b1, dj,k
1677 else
1678 write(iprint, 1500) j,id1,id2,key,jstat,xj,cj, dj,k
1679 end if
1680 end if
1681 end if
1682 return
1683
1684 1000 format(i8, 2x, 2a4, 1x, a1, 1x, a3, 1p, e16.6, 4e16.6, i7)
1685 1200 format(i8, 2x, 2a4, 1x, a1, 1x, a3, 5f16.5, i7)
1686 1300 format(i8, 2x, 2a4, 1x, a1, 1x, a3, 2f16.5,
1687 $ ' NONE ', f16.5, f16.5, i7)
1688 1400 format(i8, 2x, 2a4, 1x, a1, 1x, a3, 2f16.5,
1689 $ f16.5, ' NONE ', f16.5, i7)
1690 1500 format(i8, 2x, 2a4, 1x, a1, 1x, a3, 2f16.5,
1691 $ ' NONE ', ' NONE ', f16.5, i7)
1692
1693 * end of m4solp
1694 end
1695
1696 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1697
1698 subroutine m4stat( k, istate )
1699
1700 integer k
1701 character*12 istate
1702
1703 * ------------------------------------------------------------------
1704 * m4stat loads istate with words describing the current state.
1705 * ------------------------------------------------------------------
1706
1707 intrinsic min
1708 character*12 c(0:5)
1709 data c /'PROCEEDING ',
1710 $ 'OPTIMAL SOLN',
1711 $ 'INFEASIBLE ',
1712 $ 'UNBOUNDED ',
1713 $ 'EXCESS ITNS ',
1714 $ 'ERROR CONDN '/
1715
1716 j = min( k, 5 )
1717 istate = c(j)
1718
1719 * end of m4stat
1720 end

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