| 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 |
Parent Directory
| 
Revision Log