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

Contents of /trunk/minos54/mi20amat.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: 55975 byte(s)
Setting up web subdirectory in repository
1 ************************************************************************
2 *
3 * File mi20amat fortran.
4 *
5 * m2core m2amat m2aprd m2apr1 m2apr5
6 * m2crsh m2scal m2scla m2unpk matcol
7 *
8 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
9
10 subroutine m2core( mode, mincor )
11
12 implicit double precision (a-h,o-z)
13
14 * ------------------------------------------------------------------
15 * m2core allocates all array storage for MINOS,
16 * using various dimensions stored in common blocks as follows:
17 * (m2len ) mrows, mcols, melms
18 * (m3len ) nscl
19 * (m5len ) maxr, maxs, nn
20 * (m7len ) nnobj
21 * (m8len ) njac, nncon, nnjac
22 *
23 * if mode = 1, the call is from minos1 to estimate the storage
24 * needed for all stages of the problem.
25 * if mode = 2, the call is from m3inpt to find the minimum storage
26 * needed for mps input.
27 * if mode = 3, all dimensions should be known exactly.
28 * The call is from m3inpt to find the minimum
29 * storage to solve the problem and output the solution.
30 * The basis factorization routines will say later
31 * if they have sufficient storage.
32 * if mode = 4, the call is from minoss and all dimensions are known
33 * exactly. a, ha, ka etc need not be allocated.
34 * ------------------------------------------------------------------
35
36 common /m1word/ nwordr,nwordi,nwordh
37 common /m2len / mrows,mcols,melms
38 common /m2mapa/ ne ,nka ,la ,lha ,lka
39 common /m2mapz/ maxw ,maxz
40 common /m3len / m ,n ,nb ,nscl
41 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
42 $ lhrtyp,lhs ,lkb
43 common /m3mps1/ lname1,lname2,lkeynm,nname
44 common /m3scal/ sclobj,scltol,lscale
45 common /m5len / maxr ,maxs ,mbs ,nn ,nn0 ,nr ,nx
46 common /m5loc / lpi ,lpi2 ,lw ,lw2 ,
47 $ lx ,lx2 ,ly ,ly2 ,
48 $ lgsub ,lgsub2,lgrd ,lgrd2 ,
49 $ lr ,lrg ,lrg2 ,lxn
50 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
51 common /m7loc / lgobj ,lgobj2
52 common /m7cg2 / lcg1,lcg2,lcg3,lcg4,modtcg,nitncg,nsubsp
53 common /m8len / njac ,nncon ,nncon0,nnjac
54 common /m8loc / lfcon ,lfcon2,lfdif ,lfdif2,lfold ,
55 $ lblslk,lbuslk,lxlam ,lrhs ,
56 $ lgcon ,lgcon2,lxdif ,lxold
57 common /m8al1 / penpar,rowtol,ncom,nden,nlag,nmajor,nminor
58 * ------------------------------------------------------------------
59
60 m = mrows
61 n = mcols
62 ne = melms
63 mbs = m + maxs
64 nka = n + 1
65 nb = n + m
66 nscl = nb
67 if (lscale .eq. 0) nscl = 1
68 nname = nb
69 if (mode .eq. 4) nname = 1
70
71 ncoli = n/nwordi + 1
72 nrowi = m/nwordi + 1
73 lenh = max( 3*nrowi, 100 )
74
75 nn0 = max( nn , 1 )
76 nncon0 = max( nncon, 1 )
77 nnobj0 = max( nnobj, 1 )
78 maxr = max( maxr , 1 )
79 nr = maxr*(maxr + 1)/2 + (maxs - maxr)
80 nx = max( mbs, nn )
81 nx2 = 0
82 if (lscale .ge. 2) nx2 = nn
83
84 if (mode .le. 3) then
85
86 * MINOS is getting its data from an MPS file.
87 * Allocate arrays that are arguments of minoss and misolv.
88 * These are for the data,
89 * A, ha, ka, bl, bu, name1, name2,
90 * and for the solution
91 * hs, xn, pi, rc.
92
93 la = maxw + 1
94 lha = la + ne
95 lka = lha + 1 + (ne /nwordh)
96 lbl = lka + 1 + (nka/nwordi)
97 lbu = lbl + nb
98 lname1 = lbu + nb
99 lname2 = lname1 + 1 + (nname/nwordi)
100 lhs = lname2 + 1 + (nname/nwordi)
101 lxn = lhs + 1 + (nb /nwordh)
102 lpi = lxn + nb
103 lrc = lpi + m
104 minsub = lrc + nb
105 else
106
107 * The subroutine version can use all of z except up to maxw.
108
109 minsub = maxw + 1
110 end if
111
112 * Allocate arrays needed during and after MPS input.
113
114 lhrtyp = minsub
115 lkb = lhrtyp + 1 + (mbs/nwordh)
116 minz = lkb + 1 + (mbs/nwordi)
117
118 if (mode .le. 2) then
119
120 * Allocate additional arrays needed for MPS input.
121 * Currently just keynam -- the hash table.
122
123 lkeynm = minz
124 minmps = lkeynm + lenh
125 end if
126
127 * LP and reduced-gradient algorithm.
128 ***** Beware -- length 0000 is used below for arrays that are
129 ***** not needed in this version of MINOS.
130
131 lascal = minz
132 lpi2 = lascal + nscl
133 lw = lpi2 + 0000
134 lw2 = lw + nx
135 lx = lw2 + 0000
136 lx2 = lx + nx
137 ly = lx2 + nx2
138 ly2 = ly + nx
139 lgsub = ly2 + nx
140 lgsub2 = lgsub + nn
141 lgrd = lgsub2 + 0000
142 lr = lgrd + mbs
143 lrg = lr + nr
144 lrg2 = lrg + maxs
145 minz = lrg2 + maxs
146
147 * Nonlinear objective.
148
149 lgobj = minz
150 lgobj2 = lgobj + nnobj
151 minz = lgobj2 + nnobj
152
153 * Nonlinear constraints.
154 * If Jacobian = Dense we can set njac correctly.
155 * Otherwise we have to guess if the MPS file hasn't been
156 * input yet (mode = 1 or 2).
157
158 if (nden .eq. 1) then
159 njac = nncon * nnjac
160 else if (mode .le. 2) then
161 njac = 5 * nnjac
162 end if
163 njac = max( njac, 1 )
164 lfcon = minz
165 lfcon2 = lfcon + nncon
166 lfdif = lfcon2 + nncon
167 lfdif2 = lfdif + nncon
168 lfold = lfdif2 + 0000
169 lblslk = lfold + nncon
170 lbuslk = lblslk + nncon
171 lxlam = lbuslk + nncon
172 lrhs = lxlam + nncon
173 lgcon = lrhs + nncon
174 lgcon2 = lgcon + njac
175 lxdif = lgcon2 + njac
176 lxold = lxdif + nn
177 minz = lxold + nn
178
179 * Work arrays that can be overwritten by the names (see below)
180 * in between cycles.
181
182 lbbl = minz
183 lbbu = lbbl + mbs
184 lgrd2 = lbbu + mbs
185 minz = lgrd2 + mbs
186
187 * Allocate arrays for the basis factorization routines.
188 * minz points to the beginning of the LU factorization.
189 * This is the beginning of free space between cycles
190 * if the LU factors are allowed to be overwritten.
191
192 call m2bmap( mode, m, n, ne, minz, maxz, nguess )
193
194 if (mode .eq. 1) mincor = max( minmps, nguess, minz )
195 if (mode .eq. 2) mincor = minmps
196 if (mode .ge. 3) mincor = minz
197
198 * end of m2core
199 end
200
201 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
202
203 subroutine m2amat( mode, m, n, nb,
204 $ ne, nka, a, ha, ka,
205 $ bl, bu, hrtype )
206
207 implicit double precision (a-h,o-z)
208 integer*4 ha(ne), hrtype(m)
209 integer ka(nka)
210 double precision a(ne), bl(nb), bu(nb)
211
212 * ------------------------------------------------------------------
213 * If mode = 1 or 2, m2amat defines hrtype, the set of row types.
214 * If mode = 1, m2amat also prints the matrix statistics.
215 *
216 * The vector of row-types is as follows:
217 * hrtype(i) = 0 for e rows (equalities),
218 * hrtype(i) = 2 for n rows (objective or free rows),
219 * hrtype(i) = 1 otherwise.
220 * They are used in m2scal and m2crsh.
221 * ------------------------------------------------------------------
222
223 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
224 common /m1file/ iread,iprint,isumm
225 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
226 common /m8len / njac ,nncon ,nncon0,nnjac
227
228 intrinsic abs, max, min
229 parameter ( zero = 0.0d+0 )
230
231 bplus = 0.9*plinfy
232 bminus = - bplus
233
234 * Construct the vector of row-types.
235
236 nfixed = 0
237 nfree = 0
238 nnorml = 0
239 do 100 i = 1, m
240 j = n + i
241 b1 = bl(j)
242 b2 = bu(j)
243 itype = 1
244 if (b1 .eq. b2 ) itype = 0
245 if (b1 .le. bminus .and. b2 .ge. bplus) itype = 2
246 hrtype(i) = itype
247 if (itype .eq. 0) nfixed = nfixed + 1
248 if (itype .eq. 2) nfree = nfree + 1
249 if (itype .ne. 1) go to 100
250 if (b1 .le. bminus .or. b2 .ge. bplus) nnorml = nnorml + 1
251 100 continue
252
253 if (mode .eq. 2) return
254 if (iprint .le. 0) return
255
256 nbnded = m - nfixed - nfree - nnorml
257 write(iprint, 2200)
258 write(iprint, 2300) ' Rows ', m, nnorml, nfree, nfixed, nbnded
259
260 nfixed = 0
261 nfree = 0
262 nnorml = 0
263
264 do 200 j = 1, n
265 b1 = bl(j)
266 b2 = bu(j)
267 if (b1 .eq. b2) then
268 nfixed = nfixed + 1
269 else
270 if (b1 .eq. zero ) then
271 if (b2 .ge. bplus) then
272 nnorml = nnorml + 1
273 end if
274 else if (b1 .le. bminus) then
275 if (b2 .eq. zero ) then
276 nnorml = nnorml + 1
277 else if (b2 .ge. bplus) then
278 nfree = nfree + 1
279 end if
280 end if
281 end if
282 200 continue
283
284 nbnded = n - nfixed - nfree - nnorml
285 write(iprint, 2300) ' Columns', n, nnorml, nfree, nfixed, nbnded
286
287 * Find the biggest and smallest elements in a, excluding free rows
288 * and fixed columns.
289 * Also find the largest objective coefficients in non-fixed columns.
290
291 aijmax = zero
292 aijmin = bplus
293 cmax = zero
294 cmin = bplus
295 ncost = 0
296
297 do 300 j = 1, n
298 if (bl(j) .lt. bu(j)) then
299 do 250 k = ka(j), ka(j+1) - 1
300 i = ha(k)
301 if (hrtype(i) .eq. 2) then
302 if (i .eq. iobj) then
303 aij = abs( a(k) )
304 if (aij .gt. zero) then
305 ncost = ncost + 1
306 cmax = max( cmax, aij )
307 cmin = min( cmin, aij )
308 end if
309 end if
310 else
311 aij = abs( a(k) )
312 aijmax = max( aijmax, aij )
313 aijmin = min( aijmin, aij )
314 end if
315 250 continue
316 end if
317 300 continue
318
319 if (aijmin .eq. bplus) aijmin = zero
320 adnsty = 100.0*ne / (m*n)
321 write(iprint, 2400) ne, adnsty, aijmax, aijmin, ncost
322 if (ncost .gt. 0) write(iprint, 2410) cmax, cmin
323
324 * Print a few things to be gathered as statistics
325 * from a bunch of test runs.
326
327 write(iprint, 2500) nncon, m-nncon, n
328 return
329
330 2200 format(///
331 $ ' Matrix Statistics' /
332 $ ' -----------------' /
333 $ 15x, 'Total', 6x, 'Normal', 8x, 'Free', 7x, 'Fixed',
334 $ 5x, 'Bounded')
335 2300 format(a, 5i12)
336 2400 format(/ ' No. of matrix elements', i21, 5x, 'Density', f12.3
337 $ / ' Biggest ', 1p, e35.4, ' (excluding fixed columns,'
338 $ / ' Smallest', e35.4, ' free rows, and RHS)'
339 $ // ' No. of objective coefficients', i14)
340 2410 format( ' Biggest ', 1p, e35.4, ' (excluding fixed columns)'
341 $ / ' Smallest', e35.4)
342 2500 format(//' Number of Nonlinear constraints', i12
343 $ / ' Number of Linear constraints', i12
344 $ / ' Number of Variables ', i12)
345
346 * end of m2amat
347 end
348
349 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
350
351 subroutine m2aprd( mode, x, lenx, y, leny,
352 $ ne, nka, a, ha, ka,
353 $ z, nwcore )
354
355 implicit double precision (a-h,o-z)
356 integer*4 ha(ne)
357 integer ka(nka)
358 double precision a(ne)
359 double precision x(lenx), y(leny), z(nwcore)
360
361 * ------------------------------------------------------------------
362 * m2aprd computes various products involving the matrix A.
363 * For mode 1 to 4, x is the basic and/or superbasic variables.
364 * For mode 5 to 8, x is xn, the variables in natural order.
365 *
366 * mode function lenx leny called by
367 * ---- -------- ---- ---- ---------
368 * 1 y = y - B*x m m not used
369 * 1 y = y - (B S)*x ms m not used
370 * 2 y = y - S*x ns m m7rgit
371 * 3 y = y - B(t)*x m m not used
372 * 3 y = y - (B S)(t)*x m ms not used
373 * 4 y = y - S(t)*x m ms m7chzq, m7rg
374 *
375 * 5 y = y - A*xn n m m5setx, m5solv
376 * 6 not used
377 * 7 y = y - (nonlinear A)*xn nnjac nncon m6fun1, m8setj
378 * 8 y = y - (linear A)*xn n nncon m8viol
379 * mode 7 and 8 deal with nonlinear rows only and ignore slacks.
380 *
381 * 14 Oct 1991: a, ha, ka added as parameters to help with minoss.
382 * ------------------------------------------------------------------
383
384 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
385 common /m3len / m ,n ,nb ,nscl
386 common /m3loc / lascal,lbl ,lbu ,lbbl ,lbbu ,
387 $ lhrtyp,lhs ,lkb
388 common /m5len / maxr ,maxs ,mbs ,nn ,nn0 ,nr ,nx
389 common /m8len / njac ,nncon ,nncon0,nnjac
390
391 tolz = eps0
392
393 if (mode .lt. 5) then
394 call m2apr1( mode, m, mbs, n, tolz,
395 $ ne, nka, a, ha, ka, z(lkb),
396 $ x, lenx, y, leny )
397 else
398 call m2apr5( mode, m, n, nb, nncon, nnjac, tolz,
399 $ ne, nka, a, ha, ka,
400 $ x, lenx, y, leny )
401 end if
402
403 * end of m2aprd
404 end
405
406 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
407
408 subroutine m2apr1( mode, m, mbs, n, tolz,
409 $ ne, nka, a, ha, ka, kb, x, lenx, y, leny )
410
411 implicit double precision (a-h,o-z)
412 integer*4 ha(ne)
413 integer ka(nka), kb(mbs)
414 double precision a(ne)
415 double precision x(lenx), y(leny)
416
417 * ------------------------------------------------------------------
418 * m2apr1 computes various matrix-vector products involving
419 * B and S, the basic and superbasic columns of A.
420 *
421 * mode=1 y = y - B*x or y - (B S)*x (depending on lenx)
422 * mode=2 y = y - S*x
423 * mode=3 y = y - B(transpose)*x or y - (B S)(transpose)*x
424 * (mode 3 is not used as yet)
425 * mode=4 y = y - S(transpose)*x
426 * ------------------------------------------------------------------
427
428 kbase = 0
429 if (mode .le. 2) then
430 * -------------
431 * mode 1 and 2.
432 * -------------
433 if (mode .eq. 2) kbase = m
434
435 do 200 k = 1, lenx
436 xj = x(k)
437 if (abs(xj) .le. tolz) go to 200
438 j = kb(kbase + k)
439 if (j .le. n) then
440 do 150 i = ka(j), ka(j+1) - 1
441 ir = ha(i)
442 y(ir) = y(ir) - a(i)*xj
443 150 continue
444 else
445 ir = j - n
446 y(ir) = y(ir) - xj
447 end if
448 200 continue
449 else
450 * -------------
451 * mode 3 and 4.
452 * -------------
453 if (mode .eq. 4) kbase = m
454
455 do 400 k = 1, leny
456 t = y(k)
457 j = kb(kbase + k)
458 if (j .le. n) then
459 do 350 i = ka(j), ka(j+1) - 1
460 ir = ha(i)
461 t = t - x(ir)*a(i)
462 350 continue
463 y(k) = t
464 else
465 ir = j - n
466 y(k) = t - x(ir)
467 end if
468 400 continue
469 end if
470
471 * end of m2apr1
472 end
473
474 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
475
476 subroutine m2apr5( mode, m, n, nb, nncon, nnjac, tolz,
477 $ ne, nka, a, ha, ka, xn, lenx, y, leny )
478
479 implicit double precision (a-h,o-z)
480 integer*4 ha(ne)
481 integer ka(nka)
482 double precision a(ne)
483 double precision xn(lenx), y(leny)
484
485 * ------------------------------------------------------------------
486 * m2apr5 computes products involving parts of A and xn.
487 * ------------------------------------------------------------------
488
489 intrinsic abs
490
491 if (mode .eq. 5) then
492 * ----------------------------------
493 * mode = 5. Set y = y - A*xn.
494 * ----------------------------------
495 do 500 j = 1, n
496 xj = xn(j)
497 if (abs( xj ) .gt. tolz) then
498 do 450 k = ka(j), ka(j+1) - 1
499 ir = ha(k)
500 y(ir) = y(ir) - a(k)*xj
501 450 continue
502 end if
503 500 continue
504
505 else
506 * ------------------------------------------------------------
507 * If mode = 7, set y = y - (Jacobian)*xn.
508 * If mode = 8, set y = y - (linear A)*xn.
509 * Only the first nncon rows are involved, and no slacks.
510 * ------------------------------------------------------------
511 if (mode .eq. 7) then
512 j1 = 1
513 j2 = nnjac
514 else
515 j1 = nnjac + 1
516 j2 = n
517 end if
518
519 do 800 j = j1, j2
520 xj = xn(j)
521 if (abs( xj ) .gt. tolz) then
522 do 750 k = ka(j), ka(j+1) - 1
523 ir = ha(k)
524 if (ir .le. nncon) y(ir) = y(ir) - a(k)*xj
525 750 continue
526 end if
527 800 continue
528 end if
529
530 * end of m2apr5
531 end
532
533 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
534
535 subroutine m2crsh( lcrash, m , n , nb , nn ,
536 $ ne , nka , a , ha , ka ,
537 $ hpiv , hs , hrtype, bl , bu ,
538 $ xn , z , nwcore )
539
540 implicit double precision (a-h,o-z)
541 integer*4 ha(ne), hpiv(m), hs(nb), hrtype(m)
542 integer ka(nka)
543 double precision a(ne), bl(nb), bu(nb), xn(nb), z(nwcore )
544
545 * ------------------------------------------------------------------
546 * m2crsh looks for a basis in the columns of ( A I ).
547 *
548 * ON ENTRY
549 *
550 * iparm(1) = Crash option (icrash) has been used by m4getb
551 * to set lcrash.
552 * dparm(5) = Crash tolerance (tcrash). Default = 0.1
553 *
554 * lcrash specifies the action to be taken by Crash.
555 * 0,1,2,3 The call is from m4getb.
556 * 4 The call is from m5solv.
557 * 5 The call is from m8setj.
558 *
559 * 0 The all-slack basis is set up.
560 *
561 * 1 A triangular crash is applied to the columns of A.
562 * hs(1:n) is used to help select columns.
563 * tcrash is used to ignore small entries in each column.
564 * Depending on the size of tcrash, the resulting basis
565 * will be nearly (but not strictly) lower triangular.
566 *
567 * 2 As for 1, but nonlinear rows are ignored.
568 *
569 * 3 As for 2, but linear LG rows are also ignored.
570 *
571 * 4 Linear LG rows are now included.
572 * All hs(1:nb) and xn(n+i) are defined.
573 * Slack values of xn(n+i) are used to select LG rows.
574 *
575 * 5 Nonlinear rows are now included.
576 *
577 * hrtype(*) should be defined as described in m2amat:
578 * hrtype(i) = 0 for E rows (equalities)
579 * hrtype(i) = 1 for L or G rows (inequalities)
580 * hrtype(i) = 2 for N rows (objective or free rows)
581 *
582 * xn If lcrash <= 4, xn(1:n) is used to initialize
583 * slacks as xn(n+1:nb) = - A*xn.
584 * Used to select slacks from LG rows to be in B (basis).
585 * If lcrash = 5, xn(n+1:n+nncon) contains slack values
586 * evaluated from xn(1:n) and fcon(*).
587 * Used to select slacks from nonlinear rows to be in B.
588 *
589 * hs If lcrash = 1, 2 or 3, hs(1:n) is used.
590 * If lcrash = 4 or 5, hs(1:nb) is used.
591 * If hs(j) = 0, 1 or 3, column j is eligible for B,
592 * with 3 being "preferred".
593 * If hs(j) = 2, 4 or 5, column j is ignored.
594 *
595 *
596 * Crash has several stages.
597 *
598 * Stage 1: Insert any slacks (N, L or G rows, hrtype = 1 or 2).
599 *
600 * Stage 2: Do triangular crash on any free columns (wide bounds)
601 *
602 * Stage 3: Do triangular crash on "preferred" columns (hs(j) < 0).
603 * For the linear crash, this includes variables set
604 * between their bounds in the MPS file via FR INITIAL.
605 * For the nonlinear crash, it includes nonbasics
606 * between their bounds.
607 * (That is, "pegged" variables in both cases.)
608 *
609 * Stage 4: Grab unit columns.
610 *
611 * Stage 5: Grab double columns.
612 *
613 * Stage 6: Do triangular crash on all columns.
614 *
615 * Slacks are then used to pad the basis.
616 *
617 *
618 * ON EXIT
619 *
620 * hs is set to denote an initial (B S N) partition.
621 * hs(j) = 3 denotes variables for the initial basis.
622 * If hs(j) = 2 still, variable j will be superbasic.
623 * If hs(j) = 4 or 5 still, it will be changed to 0 or 1
624 * by m4chek and variable j will be nonbasic.
625 *
626 * xn If lcrash <= 4, slacks xn(n+1:nb) are initialized.
627 *
628 * ------------------------------------------------------------------
629 * Nov 1986: Essentially the same as in 1976.
630 * Crash tolerance added.
631 * Attention paid to various hs values.
632 *
633 * 12 Nov 1988: After free rows and columns have been processed
634 * (stage 1 and 2), slacks on L or G rows are inserted
635 * if their rows have not yet been assigned.
636 *
637 * 28 Dec 1988: Refined as follows.
638 * Stage 1 inserts free and preferred rows (slacks).
639 * Stage 2 performs a triangular crash on free or
640 * preferred columns, ignoring free rows.
641 * Unpivoted L or G slacks are then inserted.
642 * Stage 3 performs a triangular crash on the other
643 * columns, ignoring rows whose slack is basic.
644 * (Such rows form the top part of U. The
645 * remaining rows form the beginning of L.)
646 *
647 * 30 Apr 1989: Stage 1 now also looks for singleton columns
648 * (ignoring free and preferred rows).
649 * 05 May 1989: Stage 2 doesn't insert slacks if crash option < 0.
650 *
651 * 06 Dec 1989: Stage 2, 3, 4 modified. Columns of length 2 are
652 * now treated specially.
653 *
654 * 20 Dec 1989: Stage 2 thru 5 modified. Free columns done before
655 * unit and double columns.
656 *
657 * 19 May 1992: xn now used to help initialize slacks.
658 * Stage 1 thru 7 redefined as above.
659 *
660 * 01 Jun 1992: abs used to define closeness of slacks to bounds.
661 * Unfortunately, xn(1:n) seldom has meaningful values.
662 *
663 * 02 Jun 1992: Poor performance on the larger problems.
664 * Reverted to simple approach: all slacks grabbed.
665 *
666 * 04 Jun 1992: Compromise -- Crash 3 now has 3 phases:
667 * (a) E rows.
668 * (b) LG rows.
669 * (c) Nonlinear rows.
670 * xn(1:n) should then define the slack values better
671 * for (b) and (c).
672 * ------------------------------------------------------------------
673
674 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
675 common /m1file/ iread,iprint,isumm
676 common /m2parm/ dparm(30),iparm(30)
677 common /m5log1/ idebug,ierr,lprint
678 common /m7len / fobj ,fobj2 ,nnobj ,nnobj0
679 common /m8len / njac ,nncon ,nncon0,nnjac
680
681 equivalence ( iparm(1), icrash )
682 equivalence ( dparm(5), tcrash )
683
684 intrinsic abs, max
685 double precision apiv(2)
686 integer ipiv(2)
687 parameter ( nstage = 6 )
688 integer num(nstage), stage
689 logical free, prefer, gotslk,
690 $ stage2, stage3, stage4, stage5
691 parameter ( zero = 0.0d+0,
692 $ small = 1.0d-3, big = 1.0d+4 )
693
694 if (iprint .gt. 0) then
695 if (lcrash .le. 3) write(iprint, 1000) icrash
696 if (lcrash .eq. 3 .and. nncon .lt. m) write(iprint, 1030)
697 if (lcrash .eq. 4) write(iprint, 1040)
698 if (lcrash .eq. 5) write(iprint, 1050)
699 end if
700
701 if (lcrash .le. 4) then
702
703 * Sets slacks xn(n+1:nb) = - A*xn.
704 * This is where the slacks are initialized.
705 * They may be altered later (see the end of Crash).
706
707 call dload ( m, zero, xn(n+1), 1 )
708 call m2aprd( 5, xn, n, xn(n+1), m,
709 $ ne, nka, a, ha, ka,
710 $ z, nwcore )
711 end if
712
713 * ------------------------------------------------------------------
714 * For Crash option 0, set hs(j) = 3 for all slacks and quit.
715 * ------------------------------------------------------------------
716 if (lcrash .eq. 0) then
717 call hload ( n, 0, hs , 1 )
718 call hload ( m, 3, hs(n+1), 1 )
719 go to 900
720 end if
721
722 * ------------------------------------------------------------------
723 * Crash option 1, 2 or 3. lcrash = 1, 2, 3, 4, or 5.
724 * tolslk measures closeness of slacks to bounds.
725 * i1,i2 are the first and last rows of A involved in Stage 1.
726 * ------------------------------------------------------------------
727 tolslk = 0.25
728 call iload ( nstage, 0, num, 1 )
729
730 if (lcrash .le. 3) then
731 * ---------------------------------------------------------------
732 * First call. lcrash = 1, 2 or 3.
733 * Initialize hpiv(*) for all rows and hs(*) for all slacks.
734 * ---------------------------------------------------------------
735 i1 = 1
736 if (lcrash .ge. 2) i1 = nncon + 1
737 i2 = m
738 nrows = i2 - i1 + 1
739
740 * Make sure there are no basic columns already (hs(j) = 3).
741 * If there are, make them "preferred".
742
743 do 20 j = 1, n
744 if (hs(j) .eq. 3) hs(j) = -1
745 20 continue
746
747 * Make relevant rows available: hpiv(i) = 1, hs(n+i) = 0.
748
749 if (nrows .gt. 0) then
750 call hload ( nrows, 1, hpiv(i1), 1 )
751 call hload ( nrows, 0, hs(n+i1), 1 )
752 end if
753
754 if (lcrash .eq. 1) then
755 nbasic = 0
756 else
757
758 * lcrash = 2 or 3: Insert nonlinear slacks.
759
760 nbasic = nncon
761 if (nncon .gt. 0) then
762 call hload ( nncon, 3, hpiv , 1 )
763 call hload ( nncon, 3, hs(n+1), 1 )
764 end if
765 end if
766
767 if (lcrash .eq. 3) then
768
769 * Insert linear inequality slacks (including free rows).
770
771 do 25 i = i1, m
772 if (hrtype(i) .ge. 1) then
773 nbasic = nbasic + 1
774 nrows = nrows - 1
775 hpiv(i) = 3
776 hs(n+i) = 3
777 end if
778 25 continue
779 end if
780
781 * We're done if there are no relevant rows.
782
783 if (nrows .eq. 0) go to 800
784
785 else
786 * ---------------------------------------------------------------
787 * Second or third call. lcrash = 4 or 5.
788 * Initialize hpiv(*) for all rows.
789 * hs(*) already defines a basis for the full problem,
790 * but we want to do better by including only some of the slacks.
791 * ---------------------------------------------------------------
792 if (lcrash .eq. 4) then
793 * ------------------------------------------------------------
794 * Crash on linear LG rows.
795 * ------------------------------------------------------------
796 if (nncon .eq. m) go to 900
797 i1 = nncon + 1
798 i2 = m
799
800 * Mark nonlinear rows as pivoted: hpiv(i) = 3.
801
802 nbasic = nncon
803 if (nbasic .gt. 0) then
804 call hload ( nbasic, 3, hpiv, 1 )
805 end if
806
807 * Mark linear E rows as pivoted: hpiv(i) = 3
808 * Make linear LG rows available: hpiv(i) = 1, hs(n+i) = 0.
809
810 do 30 i = i1, m
811 if (hrtype(i) .eq. 0) then
812 nbasic = nbasic + 1
813 hpiv(i) = 3
814 else
815 hpiv(i) = 1
816 hs(n+i) = 0
817 end if
818 30 continue
819
820 * Mark linear LG rows with hpiv(i) = 2
821 * if any basic columns contain a nonzero in row i.
822
823 do 50 j = 1, n
824 if (hs(j) .eq. 3) then
825 do 40 k = ka(j), ka(j+1) - 1
826 i = ha(k)
827 if (hrtype(i) .eq. 1) then
828 if (i .gt. nncon) then
829 if (a(k) .ne. zero) hpiv(i) = 2
830 end if
831 end if
832 40 continue
833 end if
834 50 continue
835
836 else
837 * ------------------------------------------------------------
838 * lcrash = 5. Crash on nonlinear rows.
839 * ------------------------------------------------------------
840 i1 = 1
841 i2 = nncon
842
843 * Mark all linear rows as pivoted: hpiv(i) = 3
844
845 nbasic = m - nncon
846 if (nbasic .gt. 0) then
847 call hload ( nbasic, 3, hpiv(nncon+1), 1 )
848 end if
849
850 * Make nonlinear rows available: hpiv(i) = 1, hs(n+i) = 0.
851
852 call hload ( nncon, 1, hpiv , 1 )
853 call hload ( nncon, 0, hs(n+1), 1 )
854
855 * Mark nonlinear rows with hpiv(i) = 2
856 * if any basic columns contain a nonzero in row i.
857
858 do 70 j = 1, n
859 if (hs(j) .eq. 3) then
860 do 60 k = ka(j), ka(j+1) - 1
861 i = ha(k)
862 if (i .le. nncon) then
863 if (a(k) .ne. zero) hpiv(i) = 2
864 end if
865 60 continue
866 end if
867 70 continue
868 end if
869 end if
870
871 * ------------------------------------------------------------------
872 * Stage 1: Insert relevant slacks (N, L or G rows, hrtype = 1 or 2).
873 * If lcrash = 4 or 5, grab them only if they are more than
874 * tolslk from their bound.
875 * ------------------------------------------------------------------
876 stage = 1
877 gotslk = lcrash .eq. 4 .or. lcrash .eq. 5
878
879 do 100 i = i1, i2
880 j = n + i
881 if (hs(j) .le. 1 .and. hrtype(i) .gt. 0) then
882 if (gotslk) then
883 d1 = xn(j) - bl(j)
884 d2 = bu(j) - xn(j)
885 if (min( d1, d2 ) .le. tolslk) then
886
887 * The slack is close to a bound or infeasible.
888 * Move it exactly onto the bound.
889
890 if (d1 .le. d2) then
891 xn(j) = bl(j)
892 hs(j) = 0
893 else
894 xn(j) = bu(j)
895 hs(j) = 1
896 end if
897 go to 100
898 end if
899 end if
900
901 nbasic = nbasic + 1
902 num(stage) = num(stage) + 1
903 hpiv(i) = 3
904 hs(j) = 3
905 end if
906 100 continue
907
908 if (nbasic .eq. m) go to 700
909
910 * ------------------------------------------------------------------
911 * Apply a triangular crash to various subsets of the columns of A.
912 *
913 * hpiv(i) = 1 if row i is unmarked (initial state).
914 * hpiv(i) = 3 if row i has been given a pivot
915 * in one of a set of triangular columns.
916 * hpiv(i) = 2 if one of the triangular columns contains
917 * a nonzero in row i below the triangle.
918 * ------------------------------------------------------------------
919
920 do 600 stage = 2, nstage
921 stage2 = stage .eq. 2
922 stage3 = stage .eq. 3
923 stage4 = stage .eq. 4
924 stage5 = stage .eq. 5
925
926 * ---------------------------------------------------------------
927 * Main loop for triangular crash.
928 * ---------------------------------------------------------------
929 do 200 j = 1, n
930 js = hs(j)
931 if (js .gt. 1 ) go to 200
932 if (bl(j) .eq. bu(j)) go to 200
933
934 if ( stage2 ) then
935 free = bl(j) .le. - big .and. bu(j) .ge. big
936 if ( .not. free ) go to 200
937
938 else if ( stage3 ) then
939 prefer = js .lt. 0
940 if ( .not. prefer) go to 200
941 end if
942
943 * Find the biggest aij, ignoring free rows.
944
945 k1 = ka(j)
946 k2 = ka(j+1) - 1
947 aimax = zero
948
949 do 130 k = k1, k2
950 i = ha(k)
951 if (hrtype(i) .ne. 2) then
952 ai = a(k)
953 aimax = max( aimax, abs( ai ) )
954 end if
955 130 continue
956
957 * Prevent small pivots if crash tol is too small.
958
959 if (aimax .le. small) go to 200
960
961 * Find the biggest pivots in rows that are still
962 * unpivoted and unmarked. Ignore smallish elements.
963 * nz counts the number of relevant nonzeros.
964
965 aitol = aimax * tcrash
966 nz = 0
967 npiv = 0
968 ipiv(1) = 0
969 ipiv(2) = 0
970 apiv(1) = zero
971 apiv(2) = zero
972
973 do 150 k = k1, k2
974 i = ha(k)
975 if (hs(n+i) .ne. 3) then
976 ai = abs( a(k) )
977 if (ai .gt. aitol) then
978 nz = nz + 1
979 ip = hpiv(i)
980 if (ip .le. 2) then
981 if (apiv(ip) .lt. ai) then
982 apiv(ip) = ai
983 ipiv(ip) = i
984 end if
985 else
986 npiv = npiv + 1
987 end if
988 end if
989 end if
990 150 continue
991
992 * Grab unit or double columns.
993
994 if ( stage4 ) then
995 if (nz .ne. 1) go to 200
996 else if ( stage5 ) then
997 if (nz .ne. 2) go to 200
998 end if
999
1000 * See if the column contained a potential pivot.
1001 * An unmarked row is favored over a marked row.
1002
1003 ip = 1
1004 if (ipiv(1) .eq. 0 .and. npiv .eq. 0) ip = 2
1005 i = ipiv(ip)
1006
1007 if (i .gt. 0) then
1008 nbasic = nbasic + 1
1009 num(stage) = num(stage) + 1
1010 hpiv(i) = 3
1011 hs(j) = 3
1012 if (nbasic .ge. m) go to 700
1013
1014 * Mark off all relevant unmarked rows.
1015
1016 do 180 k = k1, k2
1017 i = ha(k)
1018 if (hs(n+i) .ne. 3) then
1019 ai = abs( a(k) )
1020 if (ai .gt. aitol) then
1021 if (hpiv(i) .eq. 1) hpiv(i) = 2
1022 end if
1023 end if
1024 180 continue
1025 end if
1026 200 continue
1027 600 continue
1028
1029 * ------------------------------------------------------------------
1030 * All stages finished.
1031 * Fill remaining gaps with slacks.
1032 * ------------------------------------------------------------------
1033 700 npad = m - nbasic
1034 if (iprint .gt. 0) write(iprint, 1200) num, npad
1035
1036 if (npad .gt. 0) then
1037 do 720 i = 1, m
1038 if (hpiv(i) .lt. 3) then
1039 nbasic = nbasic + 1
1040 hs(n+i) = 3
1041 if (nbasic .ge. m) go to 800
1042 end if
1043 720 continue
1044 end if
1045
1046 * ------------------------------------------------------------------
1047 * Make sure there aren't lots of nonbasic slacks floating off
1048 * their bounds. They could take lots of iterations to move.
1049 * ------------------------------------------------------------------
1050 800 do 850 i = i1, i2
1051 j = n + i
1052 if (hs(j) .le. 1 .and. hrtype(i) .gt. 0) then
1053 d1 = xn(j) - bl(j)
1054 d2 = bu(j) - xn(j)
1055 if (min( d1, d2 ) .le. tolslk) then
1056
1057 * The slack is close to a bound or infeasible.
1058 * Move it exactly onto the bound.
1059
1060 if (d1 .le. d2) then
1061 xn(j) = bl(j)
1062 hs(j) = 0
1063 else
1064 xn(j) = bu(j)
1065 hs(j) = 1
1066 end if
1067 end if
1068 end if
1069 850 continue
1070
1071 900 return
1072
1073 1000 format(// ' Crash option', i3)
1074 1030 format( ' Crash on linear E rows:')
1075 1040 format( ' Crash on linear LG rows:')
1076 1050 format( ' Crash on nonlinear rows:')
1077 1200 format(
1078 $ ' Slacks', i6, ' Free cols', i6, ' Preferred', i6
1079 $ / ' Unit ', i6, ' Double ', i6, ' Triangle ', i6,
1080 $ ' Pad', i6)
1081
1082 * end of m2crsh
1083 end
1084
1085 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1086
1087 subroutine m2scal( m, n, nb, ne, nka, nn, nncon, nnjac,
1088 $ hrtype, ha, ka, a, ascale, bl, bu, rmin, rmax )
1089
1090 implicit double precision (a-h,o-z)
1091 integer*4 hrtype(m), ha(ne)
1092 integer ka(nka)
1093 double precision a(ne), ascale(nb), bl(nb), bu(nb)
1094 double precision rmin(m), rmax(m)
1095
1096 * ------------------------------------------------------------------
1097 * m2scal computes scale factors ascale from A, bl, bu.
1098 *
1099 * In phase 1, an iterative procedure based on geometric means is
1100 * used to compute scales from a alone. This procedure is derived
1101 * from a routine written by Robert Fourer, 1979. The main steps
1102 * are:
1103 *
1104 * (1) Compute aratio = max(i1,i2,j) a(i1,j) / a(i2,j).
1105 * (2) Divide each row i by
1106 * ( min(j) a(i,j) * max(j) a(i,j) ) ** 1/2.
1107 * (3) Divide each column j by
1108 * ( min(i) a(i,j) * max(i) a(i,j) ) ** 1/2.
1109 * (4) Compute sratio as in (1).
1110 * (5) If sratio .lt. scltol * aratio, repeat from step (1).
1111 *
1112 * Free rows (hrtype=2) and fixed columns (bl=bu) are not used
1113 * at this stage.
1114 *
1115 * In phase 2, the scales for free rows are set to be their largest
1116 * element.
1117 *
1118 * In phase 3, fixed columns are summed in order to compute
1119 * a scale factor sigma that allows for the effective rhs of the
1120 * constraints. All scales are then multiplied by sigma.
1121 *
1122 *
1123 * If lscale = 1, the first nncon rows and the first nn columns will
1124 * retain scales of 1.0 during phases 1-2, and phase 3 will not be
1125 * performed. (This takes effect if the problem is nonlinear but
1126 * the user has specified scale linear variables only.)
1127 * However, all rows contribute to the linear column scales,
1128 * and all columns contribute to the linear row scales.
1129 *
1130 * If lscale = 2, all rows and columns are scaled. To guard against
1131 * misleadingly small Jacobians, if the maximum element in any of
1132 * the first nncon rows or the first nnjac columns is less than
1133 * smallj, the corresponding row or column retains a scale of 1.0.
1134 *
1135 * Initial version March 1981.
1136 * Feb 1983. For convenience, everything is now double-precision.
1137 * Aug 1983. Revised to damp the effect of very small elements.
1138 * On each pass, a new row or column scale will not be
1139 * smaller than sqrt(damp) times the largest (scaled)
1140 * element in that row or column.
1141 * Oct 1986. Scale option 2 and phase 3 implemented.
1142 * Nov 1989. Large scales are now reduced to at most 0.1/tolx (say)
1143 * to help prevent obvious infeasibilities when the
1144 * problem is unscaled.
1145 * ------------------------------------------------------------------
1146
1147 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1148 common /m1file/ iread,iprint,isumm
1149 common /m2parm/ dparm(30),iparm(30)
1150 common /m3scal/ sclobj,scltol,lscale
1151 common /m5log1/ idebug,ierr,lprint
1152 common /m5tols/ toldj(3),tolx,tolpiv,tolrow,rowerr,xnorm
1153
1154 intrinsic abs, min, max, sqrt
1155
1156 double precision ac, ar, amin, amax, aratio, bnd, b1, b2,
1157 $ cmin, cmax, cratio, damp, sigma,
1158 $ small, smallj, sratio
1159
1160 parameter ( zero = 0.0d+0, one = 1.0d+0,
1161 $ damp = 1.0d-4, smallj = 1.0d-2 )
1162
1163 if (iprint .gt. 0 ) write(iprint, 1000)
1164 if (scltol .ge. one) scltol = 0.99
1165 bplus = 0.1*plinfy
1166 aratio = bplus
1167 maxk = 11
1168 if (lscale .eq. 1) then
1169 * Do linear rows and columns only.
1170 istart = nncon + 1
1171 jstart = nn + 1
1172 else
1173 * Do all rows and cols if lscale .ge. 2
1174 istart = 1
1175 jstart = 1
1176 end if
1177
1178 do 50 j = 1, nb
1179 ascale(j) = one
1180 50 continue
1181
1182 if (jstart .gt. n) go to 900
1183
1184 * ------------------------------------------------------------------
1185 * Main loop for phase 1.
1186 * Only the following row-types are used:
1187 * hrtype(i) = 2 for n rows (objective or free rows),
1188 * hrtype(i) = 0 or 1 otherwise.
1189 * ------------------------------------------------------------------
1190
1191 do 400 kk = 1, maxk
1192
1193 * Find the largest column ratio.
1194 * Also set new column scales (except on pass 0).
1195
1196 npass = kk - 1
1197 amin = bplus
1198 amax = zero
1199 small = smallj
1200 sratio = one
1201
1202 do 250 j = jstart, n
1203 if (bl(j) .lt. bu(j)) then
1204 cmin = bplus
1205 cmax = zero
1206
1207 do 230 k = ka(j), ka(j+1) - 1
1208 i = ha(k)
1209 if (hrtype(i) .eq. 2 ) go to 230
1210 ar = abs( a(k) )
1211 if (ar .eq. zero) go to 230
1212 ar = ar / ascale(n+i)
1213 cmin = min( cmin, ar )
1214 cmax = max( cmax, ar )
1215 230 continue
1216
1217 ac = max( cmin, damp*cmax )
1218 ac = sqrt( ac ) * sqrt( cmax )
1219 if (j .gt. nnjac) small = zero
1220 if (cmax .le. small) ac = one
1221 if (npass .gt. 0 ) ascale(j) = ac
1222 amin = min( amin, cmin / ascale(j) )
1223 amax = max( amax, cmax / ascale(j) )
1224 cratio = cmax / cmin
1225 sratio = max( sratio, cratio )
1226 end if
1227 250 continue
1228
1229 if (iprint .gt. 0) then
1230 write(iprint, 1200) npass, amin, amax, sratio
1231 end if
1232 if (npass .ge. 3 .and.
1233 $ sratio .ge. aratio*scltol) go to 420
1234 if (kk .eq. maxk ) go to 420
1235 aratio = sratio
1236
1237
1238 * Set new row scales for the next pass.
1239
1240 if (istart .le. m) then
1241 do 300 i = istart, m
1242 rmin(i) = bplus
1243 rmax(i) = zero
1244 300 continue
1245
1246 do 350 j = 1, n
1247 if (bl(j) .lt. bu(j)) then
1248 ac = ascale(j)
1249
1250 do 330 k = ka(j), ka(j+1) - 1
1251 i = ha(k)
1252 if (i .lt. istart) go to 330
1253 ar = abs( a(k) )
1254 if (ar .eq. zero ) go to 330
1255 ar = ar / ac
1256 rmin(i) = min( rmin(i), ar )
1257 rmax(i) = max( rmax(i), ar )
1258 330 continue
1259 end if
1260 350 continue
1261
1262 do 360 i = istart, m
1263 j = n + i
1264 ar = rmax(i)
1265 if (i .le. nncon .and. ar .le. smallj) then
1266 ascale(j) = one
1267 else
1268 ac = max( rmin(i), damp*ar )
1269 ascale(j) = sqrt( ac ) * sqrt( ar )
1270 end if
1271 360 continue
1272 end if
1273 400 continue
1274
1275 * ------------------------------------------------------------------
1276 * End of main loop.
1277 * ------------------------------------------------------------------
1278
1279 * Invert the column scales, so that structurals and logicals
1280 * can be treated the same way during subsequent unscaling.
1281 * Find the min and max column scales while we're at it.
1282 * Nov 1989: nclose counts how many are "close" to 1.
1283 * For problems that are already well-scaled, it seemed sensible to
1284 * set the "close" ones exactly equal to 1.
1285 * Tried "close" = (0.5,2.0) and (0.9,1.1), but they helped only
1286 * occasionally. Have decided not to interfree.
1287
1288 420 amin = bplus
1289 amax = zero
1290 close1 = 0.5
1291 close2 = 2.0
1292 nclose = 0
1293
1294
1295 do 430 j = 1, n
1296 ac = one / ascale(j)
1297
1298 if (amin .gt. ac) then
1299 amin = ac
1300 jmin = j
1301 end if
1302 if (amax .lt. ac) then
1303 amax = ac
1304 jmax = j
1305 end if
1306
1307 if (ac .gt. close1 .and. ac .lt. close2) then
1308 nclose = nclose + 1
1309 *---- ac = one
1310 end if
1311
1312 ascale(j) = ac
1313 430 continue
1314
1315 * Remember, column scales are upside down.
1316
1317 amax = one / amax
1318 amin = one / amin
1319 k = max( 1, n - jstart + 1 )
1320 percnt = nclose * 100.0 / k
1321 if (iprint .gt. 0) then
1322 write(iprint, 1300)
1323 write(iprint, 1310) 'Col', jmax, amax, 'Col', jmin, amin,
1324 $ nclose, percnt
1325 end if
1326
1327 * ------------------------------------------------------------------
1328 * Phase 2. Deal with empty rows and free rows.
1329 * Find the min and max row scales while we're at it.
1330 * ------------------------------------------------------------------
1331 amin = bplus
1332 amax = zero
1333 imin = 0
1334 imax = 0
1335 nclose = 0
1336
1337 do 440 i = istart, m
1338 j = n + i
1339
1340 if (hrtype(i) .eq. 2) then
1341 ar = rmax(i)
1342 if (ar .eq. zero) ar = one
1343 else
1344 ar = ascale(j)
1345 if (ar .eq. zero) ar = one
1346
1347 if (amin .gt. ar) then
1348 amin = ar
1349 imin = i
1350 end if
1351 if (amax .lt. ar) then
1352 amax = ar
1353 imax = i
1354 end if
1355
1356 if (ar .gt. close1 .and. ar .lt. close2) then
1357 nclose = nclose + 1
1358 *---- ar = one
1359 end if
1360 end if
1361
1362 ascale(j) = ar
1363 440 continue
1364
1365 if (imin .eq. 0) then
1366 amin = zero
1367 amax = zero
1368 end if
1369 k = max( 1, m - istart + 1 )
1370 percnt = nclose * 100.0 / k
1371 if (iprint .gt. 0) then
1372 write(iprint, 1310) 'Row', imin, amin, 'Row', imax, amax,
1373 $ nclose, percnt
1374 end if
1375
1376 * ------------------------------------------------------------------
1377 * Phase 3.
1378 * Compute what is effectively the rhs for the constraints.
1379 * We set rmax = ( A I ) * x for fixed columns and slacks,
1380 * including positive lower bounds and negative upper bounds.
1381 * ------------------------------------------------------------------
1382 if (lscale .ge. 2) then
1383 call dload ( m, zero, rmax, 1 )
1384
1385 do 500 j = 1, nb
1386 bnd = zero
1387 b1 = bl(j)
1388 b2 = bu(j)
1389 if (b1 .eq. b2 ) bnd = b1
1390 if (b1 .gt. zero) bnd = b1
1391 if (b2 .lt. zero) bnd = b2
1392 if (bnd .eq. zero) go to 500
1393
1394 if (j .le. n ) then
1395 do 480 k = ka(j), ka(j+1) - 1
1396 i = ha(k)
1397 rmax(i) = rmax(i) + a(k) * bnd
1398 480 continue
1399 else
1400 * Free slacks never get here,
1401 * so we don't have to skip them.
1402 i = j - n
1403 rmax(i) = rmax(i) + bnd
1404 end if
1405 500 continue
1406
1407 * We don't want nonzeros in free rows to interfere.
1408
1409 do 520 i = 1, m
1410 if (hrtype(i) .eq. 2) rmax(i) = zero
1411 520 continue
1412
1413 * Scale rmax = rmax / (row scales), and use its norm sigma
1414 * to adjust both row and column scales.
1415
1416 ac = dnorm1( m, rmax, 1 )
1417 call dddiv ( m, ascale(n+1), 1, rmax, 1 )
1418 sigma = dnorm1( m, rmax, 1 )
1419 if (iprint .gt. 0) write(iprint, 1400) ac, sigma
1420 sigma = max( sigma, one )
1421 call dscal ( nb, sigma, ascale, 1 )
1422
1423 * Big scales might lead to excessive infeasibility when the
1424 * problem is unscaled. If any are too big, scale them down.
1425
1426 amax = zero
1427 do 540 j = 1, n
1428 amax = max( amax, ascale(j) )
1429 540 continue
1430
1431 do 550 i = 1, m
1432 if (hrtype(i) .ne. 2) then
1433 amax = max( amax, ascale(n + i) )
1434 end if
1435 550 continue
1436
1437 big = 0.1 / tolx
1438 sigma = big / amax
1439 if (sigma .lt. one) then
1440 call dscal ( nb, sigma, ascale, 1 )
1441 if (iprint .gt. 0) write(iprint, 1410) sigma
1442 end if
1443 end if
1444
1445 if (iparm(4) .gt. 0 .and. iprint .gt. 0) then
1446 if (istart .le. m) write(iprint,1500) (i,ascale(n+i),i=istart,m)
1447 if (jstart .le. n) write(iprint,1600) (j,ascale(j ),j=jstart,n)
1448 end if
1449
1450 900 return
1451
1452 1000 format(// ' Scaling' / ' -------'
1453 $ / ' Min elem Max elem Max col ratio')
1454 1200 format( ' After', i3, 1p, e12.2, e12.2, 0p, f20.2)
1455 1300 format(/ 12x, 'Min scale', 23x, 'Max scale', 6x,
1456 $ 'Between 0.5 and 2.0')
1457 1310 format(1x, a, i7, 1p, e10.1, 12x, a, i7, e10.1, i17, 0p, f8.1)
1458 1400 format(/ ' Norm of fixed columns and slacks', 1p, e20.1
1459 $ / ' (before and after row scaling) ', 1p, e20.1)
1460 1410 format( ' Scales are large --- reduced by ', 1p, e20.1)
1461 1500 format(// ' Row scales r(i)',
1462 $ 8x, ' a(i,j) = r(i) * scaled a(i,j) / c(j)'
1463 $ / ' ----------------' // 5(i6, g16.5))
1464 1600 format(// ' Column scales c(j)',
1465 $ 5x, ' x(j) = c(j) * scaled x(j)'
1466 $ / ' -------------------' // 5(i6, g16.5))
1467
1468 * end of m2scal
1469 end
1470
1471 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1472
1473 subroutine m2scla( mode, m, n, nb, ne, nka,
1474 $ ha, ka, a, ascale, bl, bu, pi, xn )
1475
1476 implicit double precision (a-h,o-z)
1477 integer*4 ha(ne)
1478 integer ka(nka)
1479 double precision a(ne), ascale(nb), bl(nb), bu(nb)
1480 double precision pi(m), xn(nb)
1481
1482 * ------------------------------------------------------------------
1483 * m2scla scales or unscales a, bl, bu, pi, xn using ascale.
1484 * ------------------------------------------------------------------
1485
1486 common /m1eps / eps,eps0,eps1,eps2,eps3,eps4,eps5,plinfy
1487 common /m3scal/ sclobj,scltol,lscale
1488 common /m5lobj/ sinf,wtobj,minimz,ninf,iobj,jobj,kobj
1489
1490 parameter ( one = 1.0d+0 )
1491
1492 bplus = 0.1*plinfy
1493 if (mode .eq. 1) then
1494 * ---------------------------------------------------------------
1495 * mode = 1 --- scale a, bl, bu, xn and pi.
1496 * ---------------------------------------------------------------
1497 do 150 j = 1, nb
1498 scale = ascale(j)
1499 if (j .le. n) then
1500 do 110 k = ka(j), ka(j+1) - 1
1501 i = ha(k)
1502 a(k) = a(k) * ( scale / ascale(n+i) )
1503 110 continue
1504 end if
1505 xn(j) = xn(j) / scale
1506 if (bl(j) .gt. -bplus) bl(j) = bl(j) / scale
1507 if (bu(j) .lt. bplus) bu(j) = bu(j) / scale
1508 150 continue
1509
1510 call ddscl ( m, ascale(n+1), 1, pi, 1 )
1511 if (jobj .gt. 0) sclobj = ascale(jobj)
1512 else
1513 * ---------------------------------------------------------------
1514 * mode = 2 --- unscale everything.
1515 * ---------------------------------------------------------------
1516 do 250 j = 1, nb
1517 scale = ascale(j)
1518 if (j .le. n) then
1519 do 210 k = ka(j), ka(j+1) - 1
1520 i = ha(k)
1521 a(k) = a(k) * ( ascale(n+i) / scale )
1522 210 continue
1523 end if
1524 xn(j) = xn(j) * scale
1525 if (bl(j) .gt. -bplus) bl(j) = bl(j) * scale
1526 if (bu(j) .lt. bplus) bu(j) = bu(j) * scale
1527 250 continue
1528
1529 call dddiv ( m, ascale(n+1), 1, pi, 1 )
1530 sclobj = one
1531 end if
1532
1533 * end of m2scla
1534 end
1535
1536 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1537
1538 subroutine m2unpk( jq, m, n, ne, nka, a, ha, ka, y )
1539
1540 implicit double precision (a-h,o-z)
1541 integer*4 ha(ne)
1542 integer ka(nka)
1543 double precision a(ne)
1544 double precision y(m)
1545
1546 * ------------------------------------------------------------------
1547 * m2unpk expands the jq-th column of ( A I ) into y.
1548 * ------------------------------------------------------------------
1549
1550 parameter ( zero = 0.0d+0, one = 1.0d+0 )
1551
1552 call dload ( m, zero, y, 1 )
1553
1554 if (jq .le. n) then
1555 do 10 k = ka(jq), ka(jq+1) - 1
1556 i = ha(k)
1557 y(i) = a(k)
1558 10 continue
1559 else
1560 islack = jq - n
1561 y(islack) = one
1562 end if
1563
1564 * end of m2unpk
1565 end
1566
1567 *+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1568
1569 subroutine matcol( m, n, nb, ne, nka,
1570 $ a, ha, ka, bl, bu, col, ztol )
1571
1572 implicit double precision (a-h,o-z)
1573 integer*4 ha(ne)
1574 integer ka(nka)
1575 double precision a(ne), bl(nb), bu(nb), col(m)
1576
1577 * ------------------------------------------------------------------
1578 * matcol creates a new matrix column from the dense vector col.
1579 * It becomes column number jnew, which is updated accordingly.
1580 * Elements as small as the zero tolerance ztol are not retained.
1581 * Default bounds are given to the new variable.
1582 * ------------------------------------------------------------------
1583
1584 common /m1file/ iread,iprint,isumm
1585 common /m3mps3/ aijtol,bstruc(2),mlst,mer,
1586 $ aijmin,aijmax,na0,line,ier(20)
1587 common /cyclcm/ cnvtol,jnew,materr,maxcy,nephnt,nphant,nprint
1588
1589 intrinsic abs
1590
1591 if (jnew .ge. n) go to 900
1592 jnew = jnew + 1
1593 ia = ka(jnew)
1594
1595 do 100 i = 1, m
1596 if (abs( col(i) ) .le. ztol) go to 100
1597 if (ia .gt. ne) go to 910
1598 a (ia) = col(i)
1599 ha(ia) = i
1600 ia = ia + 1
1601 100 continue
1602
1603 if (ia .eq. ka(jnew)) go to 920
1604 bl(jnew) = bstruc(1)
1605 bu(jnew) = bstruc(2)
1606 ka(jnew + 1) = ia
1607 return
1608
1609 * Too many columns.
1610
1611 900 if (iprint .gt. 0) write(iprint, 1000)
1612 if (isumm .gt. 0) write(isumm , 1000)
1613 go to 990
1614
1615 * Too many elements.
1616
1617 910 if (iprint .gt. 0) write(iprint, 1100)
1618 if (isumm .gt. 0) write(isumm , 1100)
1619 go to 990
1620
1621 * Zero column.
1622
1623 920 if (iprint .gt. 0) write(iprint, 1200)
1624 if (isumm .gt. 0) write(isumm , 1200)
1625
1626 * Error exit.
1627
1628 990 materr = materr + 1
1629 return
1630
1631 1000 format(/ ' XXX MATCOL error. Not enough Phantom columns.')
1632 1100 format(/ ' XXX MATCOL error. Not enough Phantom elements.')
1633 1200 format(/ ' XXX MATCOL error. New column of A was zero.')
1634
1635 * end of matcol
1636 end

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