1 | /* $Id: CoinPresolveForcing.cpp 1448 2011-06-19 15:34:41Z stefan $ */ |
2 | // Copyright (C) 2002, International Business Machines |
3 | // Corporation and others. All Rights Reserved. |
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
5 | |
6 | #include <stdio.h> |
7 | #include <math.h> |
8 | |
9 | #include "CoinPresolveMatrix.hpp" |
10 | #include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW |
11 | #include "CoinPresolveFixed.hpp" |
12 | #include "CoinPresolveSubst.hpp" |
13 | #include "CoinHelperFunctions.hpp" |
14 | #include "CoinPresolveUseless.hpp" |
15 | #include "CoinPresolveForcing.hpp" |
16 | #include "CoinMessage.hpp" |
17 | #include "CoinFinite.hpp" |
18 | |
19 | #if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY |
20 | #include "CoinPresolvePsdebug.hpp" |
21 | #endif |
22 | |
23 | /* |
24 | This just doesn't seem efficient, particularly when used to calculate row |
25 | bounds. Lots of extra work. |
26 | */ |
27 | void implied_bounds (const double *els, |
28 | const double *clo, const double *cup, |
29 | const int *hcol, |
30 | CoinBigIndex krs, CoinBigIndex kre, |
31 | double *maxupp, double *maxdownp, |
32 | int jcol, |
33 | double rlo, double rup, |
34 | double *iclb, double *icub) |
35 | { |
36 | if (rlo<=-PRESOLVE_INF&&rup>=PRESOLVE_INF) { |
37 | *iclb = -PRESOLVE_INF; |
38 | *icub = PRESOLVE_INF; |
39 | return; |
40 | } |
41 | bool posinf = false; |
42 | bool neginf = false; |
43 | double maxup = 0.0; |
44 | double maxdown = 0.0; |
45 | |
46 | int jcolk = -1; |
47 | |
48 | // compute sum of all bounds except for jcol |
49 | CoinBigIndex kk; |
50 | for (kk=krs; kk<kre; kk++) { |
51 | if (hcol[kk] == jcol) |
52 | jcolk = kk; |
53 | |
54 | // swap jcol with hcol[kre-1]; |
55 | // that is, consider jcol last |
56 | // this assumes that jcol occurs in this row |
57 | CoinBigIndex k = (hcol[kk] == jcol |
58 | ? kre-1 |
59 | : kk == kre-1 |
60 | ? jcolk |
61 | : kk); |
62 | |
63 | PRESOLVEASSERT(k != -1); // i.e. jcol had better be in the row |
64 | |
65 | int col = hcol[k]; |
66 | double coeff = els[k]; |
67 | double lb = clo[col]; |
68 | double ub = cup[col]; |
69 | |
70 | // quick! compute the implied col bounds before maxup/maxdown are changed |
71 | if (kk == kre-1) { |
72 | PRESOLVEASSERT(fabs(coeff) > ZTOLDP); |
73 | |
74 | double ilb = (rlo - maxup) / coeff; |
75 | bool finite_ilb = (-PRESOLVE_INF < rlo && !posinf && PRESOLVEFINITE(maxup)); |
76 | |
77 | double iub = (rup - maxdown) / coeff; |
78 | bool finite_iub = ( rup < PRESOLVE_INF && !neginf && PRESOLVEFINITE(maxdown)); |
79 | |
80 | if (coeff > 0.0) { |
81 | *iclb = (finite_ilb ? ilb : -PRESOLVE_INF); |
82 | *icub = (finite_iub ? iub : PRESOLVE_INF); |
83 | } else { |
84 | *iclb = (finite_iub ? iub : -PRESOLVE_INF); |
85 | *icub = (finite_ilb ? ilb : PRESOLVE_INF); |
86 | } |
87 | } |
88 | |
89 | if (coeff > 0.0) { |
90 | if (PRESOLVE_INF <= ub) { |
91 | posinf = true; |
92 | if (neginf) |
93 | break; // pointless |
94 | } else |
95 | maxup += ub * coeff; |
96 | |
97 | if (lb <= -PRESOLVE_INF) { |
98 | neginf = true; |
99 | if (posinf) |
100 | break; // pointless |
101 | } else |
102 | maxdown += lb * coeff; |
103 | } else { |
104 | if (PRESOLVE_INF <= ub) { |
105 | neginf = true; |
106 | if (posinf) |
107 | break; // pointless |
108 | } else |
109 | maxdown += ub * coeff; |
110 | |
111 | if (lb <= -PRESOLVE_INF) { |
112 | posinf = true; |
113 | if (neginf) |
114 | break; // pointless |
115 | } else |
116 | maxup += lb * coeff; |
117 | } |
118 | } |
119 | |
120 | // If we broke from the loop, then the bounds are infinite. |
121 | // However, since we put the column whose implied bounds we want |
122 | // to know at the end, and it doesn't matter if its own bounds |
123 | // are infinite, don't worry about breaking at the last iteration. |
124 | if (kk<kre-1) { |
125 | *iclb = -PRESOLVE_INF; |
126 | *icub = PRESOLVE_INF; |
127 | } else |
128 | PRESOLVEASSERT(jcolk != -1); |
129 | |
130 | // store row bounds |
131 | *maxupp = (posinf) ? PRESOLVE_INF : maxup; |
132 | *maxdownp = (neginf) ? -PRESOLVE_INF : maxdown; |
133 | } |
134 | |
135 | static void implied_row_bounds(const double *els, |
136 | const double *clo, const double *cup, |
137 | const int *hcol, |
138 | CoinBigIndex krs, CoinBigIndex kre, |
139 | double *maxupp, double *maxdownp) |
140 | { |
141 | int jcol = hcol[krs]; |
142 | bool posinf = false; |
143 | bool neginf = false; |
144 | double maxup = 0.0; |
145 | double maxdown = 0.0; |
146 | |
147 | int jcolk = -1; |
148 | |
149 | // compute sum of all bounds except for jcol |
150 | CoinBigIndex kk; |
151 | for (kk=krs; kk<kre; kk++) { |
152 | if (hcol[kk] == jcol) |
153 | jcolk = kk; |
154 | |
155 | // swap jcol with hcol[kre-1]; |
156 | // that is, consider jcol last |
157 | // this assumes that jcol occurs in this row |
158 | CoinBigIndex k = (hcol[kk] == jcol |
159 | ? kre-1 |
160 | : kk == kre-1 |
161 | ? jcolk |
162 | : kk); |
163 | |
164 | PRESOLVEASSERT(k != -1); // i.e. jcol had better be in the row |
165 | |
166 | int col = hcol[k]; |
167 | double coeff = els[k]; |
168 | double lb = clo[col]; |
169 | double ub = cup[col]; |
170 | |
171 | if (coeff > 0.0) { |
172 | if (PRESOLVE_INF <= ub) { |
173 | posinf = true; |
174 | if (neginf) |
175 | break; // pointless |
176 | } else |
177 | maxup += ub * coeff; |
178 | |
179 | if (lb <= -PRESOLVE_INF) { |
180 | neginf = true; |
181 | if (posinf) |
182 | break; // pointless |
183 | } else |
184 | maxdown += lb * coeff; |
185 | } else { |
186 | if (PRESOLVE_INF <= ub) { |
187 | neginf = true; |
188 | if (posinf) |
189 | break; // pointless |
190 | } else |
191 | maxdown += ub * coeff; |
192 | |
193 | if (lb <= -PRESOLVE_INF) { |
194 | posinf = true; |
195 | if (neginf) |
196 | break; // pointless |
197 | } else |
198 | maxup += lb * coeff; |
199 | } |
200 | } |
201 | // store row bounds |
202 | *maxupp = (posinf) ? PRESOLVE_INF : maxup; |
203 | *maxdownp = (neginf) ? -PRESOLVE_INF : maxdown; |
204 | } |
205 | |
206 | const char *forcing_constraint_action::name() const |
207 | { |
208 | return ("forcing_constraint_action" ); |
209 | } |
210 | // defed out to avoid compiler warning |
211 | #if 0 |
212 | static bool some_col_was_fixed(const int *hcol, CoinBigIndex krs, CoinBigIndex kre, |
213 | const double *clo, |
214 | const double *cup) |
215 | { |
216 | CoinBigIndex k; |
217 | for (k=krs; k<kre; k++) { |
218 | int jcol = hcol[k]; |
219 | if (clo[jcol] == cup[jcol]) |
220 | break; |
221 | } |
222 | return (k<kre); |
223 | } |
224 | #endif |
225 | |
226 | // |
227 | // It may be the case that the variable bounds are such that no matter what |
228 | // feasible value they take, the constraint cannot be violated; |
229 | // in this case, we obviously don't need to take it into account, and |
230 | // we just drop it as a USELESS constraint. |
231 | // |
232 | // On the other hand, it may be that the only way to satisfy a constraint |
233 | // is to jam all the constraint variables to one of their bounds; |
234 | // in this case, these variables are essentially fixed, which |
235 | // we do with a FORCING constraint. |
236 | // For now, this just tightens the bounds; subsequently the fixed variables |
237 | // will be removed, then the row will be dropped. |
238 | // |
239 | // Since both operations use similar information (the implied row bounds), |
240 | // we combine them into one presolve routine. |
241 | // |
242 | // I claim that these checks could be performed in parallel, |
243 | // that is, the tests could be carried out for all rows in parallel, |
244 | // and then the rows deleted and columns tightened afterward. |
245 | // Obviously, this is true for useless rows. |
246 | // The potential problem is forcing constraints, but I think |
247 | // that is ok. |
248 | // By doing it in parallel rather than sequentially, |
249 | // we may miss transformations due to variables that were fixed |
250 | // by forcing constraints, though. |
251 | // |
252 | // Note that both of these operations will cause problems |
253 | // if the variables in question really need to exceed their bounds in order |
254 | // to make the problem feasible. |
255 | |
256 | const CoinPresolveAction |
257 | *forcing_constraint_action::presolve(CoinPresolveMatrix *prob, |
258 | const CoinPresolveAction *next) |
259 | { |
260 | double startTime = 0.0; |
261 | int startEmptyRows=0; |
262 | int startEmptyColumns = 0; |
263 | if (prob->tuning_) { |
264 | startTime = CoinCpuTime(); |
265 | startEmptyRows = prob->countEmptyRows(); |
266 | startEmptyColumns = prob->countEmptyCols(); |
267 | } |
268 | double *clo = prob->clo_; |
269 | double *cup = prob->cup_; |
270 | double *csol = prob->sol_ ; |
271 | |
272 | const double *rowels = prob->rowels_; |
273 | const int *hcol = prob->hcol_; |
274 | const CoinBigIndex *mrstrt = prob->mrstrt_; |
275 | const int *hinrow = prob->hinrow_; |
276 | const int nrows = prob->nrows_; |
277 | |
278 | const double *rlo = prob->rlo_; |
279 | const double *rup = prob->rup_; |
280 | |
281 | // const char *integerType = prob->integerType_; |
282 | |
283 | const double tol = ZTOLDP; |
284 | const double inftol = prob->feasibilityTolerance_; |
285 | const int ncols = prob->ncols_; |
286 | |
287 | int *fixed_cols = new int[ncols]; |
288 | int nfixed_cols = 0; |
289 | |
290 | action *actions = new action [nrows]; |
291 | int nactions = 0; |
292 | |
293 | int *useless_rows = new int[nrows]; |
294 | int nuseless_rows = 0; |
295 | |
296 | int numberLook = prob->numberRowsToDo_; |
297 | int iLook; |
298 | int * look = prob->rowsToDo_; |
299 | bool fixInfeasibility = (prob->presolveOptions_&16384)!=0; |
300 | |
301 | # if PRESOLVE_DEBUG |
302 | std::cout << "Entering forcing_constraint_action::presolve." << std::endl ; |
303 | presolve_check_sol(prob) ; |
304 | # endif |
305 | /* |
306 | Open a loop to scan the constraints of interest. There must be variables |
307 | left in the row. |
308 | */ |
309 | for (iLook=0;iLook<numberLook;iLook++) { |
310 | int irow = look[iLook]; |
311 | if (hinrow[irow] > 0) { |
312 | CoinBigIndex krs = mrstrt[irow]; |
313 | CoinBigIndex kre = krs + hinrow[irow]; |
314 | /* |
315 | Calculate upper and lower bounds on the row activity based on upper and lower |
316 | bounds on the variables. If these are finite and incompatible with the given |
317 | row bounds, we have infeasibility. |
318 | */ |
319 | double maxup, maxdown; |
320 | implied_row_bounds(rowels, clo, cup, hcol, krs, kre, |
321 | &maxup, &maxdown); |
322 | |
323 | if (maxup < PRESOLVE_INF && maxup + inftol < rlo[irow]&&!fixInfeasibility) { |
324 | /* max row activity below the row lower bound */ |
325 | prob->status_|= 1; |
326 | prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, |
327 | prob->messages()) |
328 | <<irow |
329 | <<rlo[irow] |
330 | <<rup[irow] |
331 | <<CoinMessageEol; |
332 | break; |
333 | } else if (-PRESOLVE_INF < maxdown && rup[irow] < maxdown - inftol&&!fixInfeasibility) { |
334 | /* min row activity above the row upper bound */ |
335 | prob->status_|= 1; |
336 | prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, |
337 | prob->messages()) |
338 | <<irow |
339 | <<rlo[irow] |
340 | <<rup[irow] |
341 | <<CoinMessageEol; |
342 | break; |
343 | } |
344 | // ADD TOLERANCE TO THESE TESTS |
345 | else if ((rlo[irow] <= -PRESOLVE_INF || |
346 | (-PRESOLVE_INF < maxdown && rlo[irow] <= maxdown)) && |
347 | (rup[irow] >= PRESOLVE_INF || |
348 | (maxup < PRESOLVE_INF && rup[irow] >= maxup))) { |
349 | |
350 | /* |
351 | Original comment: I'm not sure that these transforms don't intefere with |
352 | each other. We can get it next time. |
353 | |
354 | Well, I'll argue that bounds are never really loosened (at worst, they're |
355 | transferred onto some other variable, or inferred to be unnecessary. |
356 | Once useless, always useless. Leaving this hook in place allows for a sort |
357 | of infinite loop where this routine keeps queuing the same constraints over |
358 | and over. -- lh, 040901 -- |
359 | |
360 | if (some_col_was_fixed(hcol, krs, kre, clo, cup)) { |
361 | prob->addRow(irow); |
362 | continue; |
363 | } |
364 | */ |
365 | |
366 | // this constraint must always be satisfied - drop it |
367 | useless_rows[nuseless_rows++] = irow; |
368 | |
369 | } else if ((maxup < PRESOLVE_INF && fabs(rlo[irow] - maxup) < tol) || |
370 | (-PRESOLVE_INF < maxdown && fabs(rup[irow] - maxdown) < tol)) { |
371 | |
372 | // the lower bound can just be reached, or |
373 | // the upper bound can just be reached; |
374 | // called a "forcing constraint" in the paper (p. 226) |
375 | const int lbound_tight = (maxup < PRESOLVE_INF && |
376 | fabs(rlo[irow] - maxup) < tol); |
377 | |
378 | /* |
379 | Original comment and rebuttal as above. |
380 | if (some_col_was_fixed(hcol, krs, kre, clo, cup)) { |
381 | // make sure on next time |
382 | prob->addRow(irow); |
383 | continue; |
384 | } |
385 | */ |
386 | // out of space - this probably never happens (but this routine will |
387 | // often put duplicates in the fixed column list) |
388 | if (nfixed_cols + (kre-krs) >= ncols) |
389 | break; |
390 | |
391 | double *bounds = new double[hinrow[irow]]; |
392 | int *rowcols = new int[hinrow[irow]]; |
393 | int lk = krs; // load fix-to-down in front |
394 | int uk = kre; // load fix-to-up in back |
395 | CoinBigIndex k; |
396 | for ( k=krs; k<kre; k++) { |
397 | int jcol = hcol[k]; |
398 | prob->addCol(jcol); |
399 | double coeff = rowels[k]; |
400 | |
401 | PRESOLVEASSERT(fabs(coeff) > ZTOLDP); |
402 | |
403 | // one of the two contributed to maxup - set the other to that |
404 | if (lbound_tight == (coeff > 0.0)) { |
405 | --uk; |
406 | bounds[uk-krs] = clo[jcol]; |
407 | rowcols[uk-krs] = jcol; |
408 | if (csol != 0) { |
409 | csol[jcol] = cup[jcol] ; |
410 | } |
411 | clo[jcol] = cup[jcol]; |
412 | } else { |
413 | bounds[lk-krs] = cup[jcol]; |
414 | rowcols[lk-krs] = jcol; |
415 | ++lk; |
416 | if (csol != 0) { |
417 | csol[jcol] = clo[jcol] ; |
418 | } |
419 | cup[jcol] = clo[jcol]; |
420 | } |
421 | |
422 | fixed_cols[nfixed_cols++] = jcol; |
423 | } |
424 | PRESOLVEASSERT(uk == lk); |
425 | |
426 | action *f = &actions[nactions]; |
427 | nactions++; |
428 | PRESOLVE_DETAIL_PRINT(printf("pre_forcing %dR E\n" ,irow)); |
429 | |
430 | f->row = irow; |
431 | f->nlo = lk-krs; |
432 | f->nup = kre-uk; |
433 | f->rowcols = rowcols; |
434 | f->bounds = bounds; |
435 | } |
436 | } |
437 | } |
438 | |
439 | |
440 | if (nactions) { |
441 | #if PRESOLVE_SUMMARY |
442 | printf("NFORCED: %d\n" , nactions); |
443 | #endif |
444 | next = new forcing_constraint_action(nactions, |
445 | CoinCopyOfArray(actions,nactions), next); |
446 | } |
447 | deleteAction(actions,action*); |
448 | if (nuseless_rows) { |
449 | next = useless_constraint_action::presolve(prob, |
450 | useless_rows, nuseless_rows, |
451 | next); |
452 | } |
453 | delete[]useless_rows; |
454 | /* |
455 | We need to remove duplicates here, or we get into trouble in |
456 | remove_fixed_action::postsolve when we try to reinstate a column multiple |
457 | times. |
458 | */ |
459 | if (nfixed_cols) |
460 | { if (nfixed_cols > 1) |
461 | { std::sort(fixed_cols,fixed_cols+nfixed_cols) ; |
462 | int *end = std::unique(fixed_cols,fixed_cols+nfixed_cols) ; |
463 | nfixed_cols = static_cast<int>(end-fixed_cols) ; } |
464 | next = remove_fixed_action::presolve(prob,fixed_cols,nfixed_cols,next) ; } |
465 | delete[]fixed_cols ; |
466 | |
467 | if (prob->tuning_) { |
468 | double thisTime=CoinCpuTime(); |
469 | int droppedRows = prob->countEmptyRows() - startEmptyRows ; |
470 | int droppedColumns = prob->countEmptyCols() - startEmptyColumns; |
471 | printf("CoinPresolveForcing(32) - %d rows, %d columns dropped in time %g, total %g\n" , |
472 | droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_); |
473 | } |
474 | |
475 | # if PRESOLVE_DEBUG |
476 | presolve_check_sol(prob) ; |
477 | std::cout << "Leaving forcing_constraint_action::presolve." << std::endl ; |
478 | # endif |
479 | |
480 | return (next); |
481 | } |
482 | |
483 | void forcing_constraint_action::postsolve(CoinPostsolveMatrix *prob) const |
484 | { |
485 | const action *const actions = actions_; |
486 | const int nactions = nactions_; |
487 | |
488 | const double *colels = prob->colels_; |
489 | const int *hrow = prob->hrow_; |
490 | const CoinBigIndex *mcstrt = prob->mcstrt_; |
491 | const int *hincol = prob->hincol_; |
492 | const int *link = prob->link_; |
493 | |
494 | // CoinBigIndex free_list = prob->free_list_; |
495 | |
496 | double *clo = prob->clo_; |
497 | double *cup = prob->cup_; |
498 | double *rlo = prob->rlo_; |
499 | double *rup = prob->rup_; |
500 | |
501 | const double *sol = prob->sol_; |
502 | double *rcosts = prob->rcosts_; |
503 | |
504 | double *acts = prob->acts_; |
505 | double *rowduals = prob->rowduals_; |
506 | |
507 | const double ztoldj = prob->ztoldj_; |
508 | const double ztolzb = prob->ztolzb_; |
509 | |
510 | for (const action *f = &actions[nactions-1]; actions<=f; f--) { |
511 | |
512 | const int irow = f->row; |
513 | const int nlo = f->nlo; |
514 | const int nup = f->nup; |
515 | const int ninrow = nlo + nup; |
516 | const int *rowcols = f->rowcols; |
517 | const double *bounds= f->bounds; |
518 | int k; |
519 | /* |
520 | Original comment: When we restore bounds here, we need to allow for the |
521 | possibility that the restored bound is infinite. This implies a check |
522 | for viable status. |
523 | |
524 | Hmmm ... I'm going to argue that in fact we have no choice: the status |
525 | of the variable must reflect the value it was fixed at, else we lose |
526 | feasibility. We don't care what the other bound does. -- lh, 040903 -- |
527 | */ |
528 | for (k=0; k<nlo; k++) { |
529 | int jcol = rowcols[k]; |
530 | cup[jcol] = bounds[k]; |
531 | prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::atLowerBound) ; |
532 | /* |
533 | PRESOLVEASSERT(prob->getColumnStatus(jcol)!=CoinPrePostsolveMatrix::basic); |
534 | if (cup[jcol] >= PRESOLVE_INF) |
535 | { CoinPrePostsolveMatrix::Status statj = prob->getColumnStatus(jcol) ; |
536 | if (statj == CoinPrePostsolveMatrix::atUpperBound) |
537 | { if (clo[jcol] > -PRESOLVE_INF) |
538 | { statj = CoinPrePostsolveMatrix::atLowerBound ; } |
539 | else |
540 | { statj = CoinPrePostsolveMatrix::isFree ; } |
541 | prob->setColumnStatus(jcol,statj) ; } } |
542 | */ |
543 | } |
544 | |
545 | for (k=nlo; k<ninrow; k++) { |
546 | int jcol = rowcols[k]; |
547 | clo[jcol] = bounds[k]; |
548 | prob->setColumnStatus(jcol,CoinPrePostsolveMatrix::atUpperBound) ; |
549 | /* |
550 | PRESOLVEASSERT(prob->getColumnStatus(jcol)!=CoinPrePostsolveMatrix::basic); |
551 | if (clo[jcol] <= -PRESOLVE_INF) |
552 | { CoinPrePostsolveMatrix::Status statj = prob->getColumnStatus(jcol) ; |
553 | if (statj == CoinPrePostsolveMatrix::atLowerBound) |
554 | { if (cup[jcol] < PRESOLVE_INF) |
555 | { statj = CoinPrePostsolveMatrix::atUpperBound ; } |
556 | else |
557 | { statj = CoinPrePostsolveMatrix::isFree ; } |
558 | prob->setColumnStatus(jcol,statj) ; } } |
559 | */ |
560 | } |
561 | |
562 | PRESOLVEASSERT(prob->getRowStatus(irow)==CoinPrePostsolveMatrix::basic); |
563 | PRESOLVEASSERT(rowduals[irow] == 0.0); |
564 | |
565 | // this is a lazy implementation. |
566 | // we tightened the col bounds, then let them be eliminated |
567 | // by repeated uses of FIX_VARIABLE and a final DROP_ROW. |
568 | // Therefore, by this point the row has been marked basic, |
569 | // the rowdual of this row is 0.0, |
570 | // and the reduced costs for the cols may or may not be ok |
571 | // for the relaxed column bounds. |
572 | // |
573 | // find the one most out of whack and fix it. |
574 | int whacked = -1; |
575 | double whack = 0.0; |
576 | for (k=0; k<ninrow; k++) { |
577 | int jcol = rowcols[k]; |
578 | CoinBigIndex kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link); |
579 | |
580 | // choose rowdual to cancel out reduced cost |
581 | double whack0 = rcosts[jcol] / colels[kk]; |
582 | |
583 | if (((rcosts[jcol] > ztoldj && !(fabs(sol[jcol] - clo[jcol]) <= ztolzb)) || |
584 | (rcosts[jcol] < -ztoldj && !(fabs(sol[jcol] - cup[jcol]) <= ztolzb))) && |
585 | fabs(whack0) > fabs(whack)) { |
586 | whacked = jcol; |
587 | whack = whack0; |
588 | } |
589 | } |
590 | |
591 | if (whacked != -1) { |
592 | prob->setColumnStatus(whacked,CoinPrePostsolveMatrix::basic); |
593 | if (acts[irow]-rlo[irow]<rup[irow]-acts[irow]) |
594 | prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound); |
595 | else |
596 | prob->setRowStatus(irow,CoinPrePostsolveMatrix::atUpperBound); |
597 | rowduals[irow] = whack; |
598 | |
599 | for (k=0; k<ninrow; k++) { |
600 | int jcol = rowcols[k]; |
601 | CoinBigIndex kk = presolve_find_row2(irow, mcstrt[jcol], hincol[jcol], hrow, link); |
602 | |
603 | rcosts[jcol] -= (rowduals[irow] * colels[kk]); |
604 | } |
605 | } |
606 | } |
607 | |
608 | # if PRESOLVE_CONSISTENCY |
609 | presolve_check_threads(prob) ; |
610 | # endif |
611 | |
612 | } |
613 | |
614 | |
615 | |
616 | #if 0 // (A) |
617 | // Determine the maximum and minimum values the constraint sums |
618 | // may take, given the bounds on the variables. |
619 | // If there are infinite terms, record where the first one is, |
620 | // and whether there is more than one. |
621 | // It is possible to compute implied bounds for the (one) variable |
622 | // with no bound. |
623 | static void implied_bounds1(CoinPresolveMatrix * prob, const double *rowels, |
624 | const int *mrstrt, |
625 | const int *hrow, |
626 | const int *hinrow, |
627 | const double *clo, const double *cup, |
628 | const int *hcol, |
629 | int ncols, |
630 | const double *rlo, const double *rup, |
631 | const char *integerType, |
632 | int nrows, |
633 | double *ilbound, double *iubound) |
634 | { |
635 | const double tol = prob->feasibilityTolerance_; |
636 | |
637 | for (int irow=0; irow<nrows; irow++) { |
638 | CoinBigIndex krs = mrstrt[irow]; |
639 | CoinBigIndex kre = krs + hinrow[irow]; |
640 | |
641 | double irlo = rlo[irow]; |
642 | double irup = rup[irow]; |
643 | |
644 | // These are used to set column bounds below. |
645 | // If there are no (positive) infinite terms, |
646 | // the loop will range from krs to kre; |
647 | // if there is just one, it will range over that one variable; |
648 | // otherwise, it will be empty. |
649 | int ub_inf_index = -1; |
650 | int lb_inf_index = -1; |
651 | |
652 | double maxup = 0.0; |
653 | double maxdown = 0.0; |
654 | CoinBigIndex k; |
655 | for (k=krs; k<kre; k++) { |
656 | int jcol = hcol[k]; |
657 | double coeff = rowels[k]; |
658 | double lb = clo[jcol]; |
659 | double ub = cup[jcol]; |
660 | |
661 | // HAVE TO DEAL WITH BOUNDS OF INTEGER VARIABLES |
662 | if (coeff > 0.0) { |
663 | if (PRESOLVE_INF <= ub) { |
664 | if (ub_inf_index == -1) { |
665 | ub_inf_index = k; |
666 | } else { |
667 | ub_inf_index = -2; |
668 | if (lb_inf_index == -2) |
669 | break; // pointless |
670 | } |
671 | } else |
672 | maxup += ub * coeff; |
673 | |
674 | if (lb <= -PRESOLVE_INF) { |
675 | if (lb_inf_index == -1) { |
676 | lb_inf_index = k; |
677 | } else { |
678 | lb_inf_index = -2; |
679 | if (ub_inf_index == -2) |
680 | break; // pointless |
681 | } |
682 | } else |
683 | maxdown += lb * coeff; |
684 | } |
685 | else { |
686 | if (PRESOLVE_INF <= ub) { |
687 | if (lb_inf_index == -1) { |
688 | lb_inf_index = k; |
689 | } else { |
690 | lb_inf_index = -2; |
691 | if (ub_inf_index == -2) |
692 | break; // pointless |
693 | } |
694 | } else |
695 | maxdown += ub * coeff; |
696 | |
697 | if (lb <= -PRESOLVE_INF) { |
698 | if (ub_inf_index == -1) { |
699 | ub_inf_index = k; |
700 | } else { |
701 | ub_inf_index = -2; |
702 | if (lb_inf_index == -2) |
703 | break; // pointless |
704 | } |
705 | } else |
706 | maxup += lb * coeff; |
707 | } |
708 | } |
709 | |
710 | // ub_inf says whether the sum of the "other" ub terms is infinite |
711 | // in the loop below. |
712 | // In the case where we only saw one infinite term, the loop |
713 | // will only cover that case, in which case the other terms |
714 | // are *not* infinite. |
715 | // With two or more terms, it is infinite. |
716 | // If we only saw one infinite term, then |
717 | if (ub_inf_index == -2) |
718 | maxup = PRESOLVE_INF; |
719 | |
720 | if (lb_inf_index == -2) |
721 | maxdown = -PRESOLVE_INF; |
722 | |
723 | const bool maxup_finite = PRESOLVEFINITE(maxup); |
724 | const bool maxdown_finite = PRESOLVEFINITE(maxdown); |
725 | |
726 | if (ub_inf_index == -1 && maxup_finite && maxup + tol < rlo[irow]&&!fixInfeasibility) { |
727 | /* infeasible */ |
728 | prob->status_|= 1; |
729 | prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, |
730 | prob->messages()) |
731 | <<irow |
732 | <<rlo[irow] |
733 | <<rup[irow] |
734 | <<CoinMessageEol; |
735 | break; |
736 | } else if (lb_inf_index == -1 && maxdown_finite && rup[irow] < maxdown - tol&&!fixInfeasibility) { |
737 | /* infeasible */ |
738 | prob->status_|= 1; |
739 | prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, |
740 | prob->messages()) |
741 | <<irow |
742 | <<rlo[irow] |
743 | <<rup[irow] |
744 | <<CoinMessageEol; |
745 | break; |
746 | } |
747 | |
748 | for (k = krs; k<kre; ++k) { |
749 | int jcol = hcol[k]; |
750 | double coeff = rowels[k]; |
751 | |
752 | // SHOULD GET RID OF THIS |
753 | if (fabs(coeff) > ZTOLDP2 && |
754 | !integerType[jcol]) { |
755 | double maxup1 = (ub_inf_index == -1 || ub_inf_index == k |
756 | ? maxup |
757 | : PRESOLVE_INF); |
758 | bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k |
759 | ? maxup_finite |
760 | : false); |
761 | double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k |
762 | ? maxdown |
763 | : PRESOLVE_INF); |
764 | bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k |
765 | ? maxdown_finite |
766 | : false); |
767 | |
768 | double ilb = (irlo - maxup1) / coeff; |
769 | bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1); |
770 | |
771 | double iub = (irup - maxdown1) / coeff; |
772 | bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1); |
773 | |
774 | double ilb1 = (coeff > 0.0 |
775 | ? (finite_ilb ? ilb : -PRESOLVE_INF) |
776 | : (finite_iub ? iub : -PRESOLVE_INF)); |
777 | |
778 | if (ilbound[jcol] < ilb1) { |
779 | ilbound[jcol] = ilb1; |
780 | //if (jcol == 278001) |
781 | //printf("JCOL LB %g\n", ilb1); |
782 | } |
783 | } |
784 | } |
785 | |
786 | for (k = krs; k<kre; ++k) { |
787 | int jcol = hcol[k]; |
788 | double coeff = rowels[k]; |
789 | |
790 | // SHOULD GET RID OF THIS |
791 | if (fabs(coeff) > ZTOLDP2 && |
792 | !integerType[jcol]) { |
793 | double maxup1 = (ub_inf_index == -1 || ub_inf_index == k |
794 | ? maxup |
795 | : PRESOLVE_INF); |
796 | bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k |
797 | ? maxup_finite |
798 | : false); |
799 | double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k |
800 | ? maxdown |
801 | : PRESOLVE_INF); |
802 | bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k |
803 | ? maxdown_finite |
804 | : false); |
805 | |
806 | |
807 | double ilb = (irlo - maxup1) / coeff; |
808 | bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1); |
809 | |
810 | double iub = (irup - maxdown1) / coeff; |
811 | bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1); |
812 | |
813 | double iub1 = (coeff > 0.0 |
814 | ? (finite_iub ? iub : PRESOLVE_INF) |
815 | : (finite_ilb ? ilb : PRESOLVE_INF)); |
816 | |
817 | if (iub1 < iubound[jcol]) { |
818 | iubound[jcol] = iub1; |
819 | //if (jcol == 278001) |
820 | //printf("JCOL UB %g\n", iub1); |
821 | } |
822 | } |
823 | } |
824 | } |
825 | } |
826 | |
827 | #if 0 // (B) |
828 | postsolve for implied_bound |
829 | { |
830 | double lo0 = pa->clo; |
831 | double up0 = pa->cup; |
832 | int irow = pa->irow; |
833 | int jcol = pa->icol; |
834 | int *rowcols = pa->rowcols; |
835 | int ninrow = pa->ninrow; |
836 | |
837 | clo[jcol] = lo0; |
838 | cup[jcol] = up0; |
839 | |
840 | if ((colstat[jcol] & PRESOLVE_XBASIC) == 0 && |
841 | fabs(lo0 - sol[jcol]) > ztolzb && |
842 | fabs(up0 - sol[jcol]) > ztolzb) { |
843 | |
844 | // this non-basic variable is now away from its bound |
845 | // it is ok just to force it to be basic |
846 | // informally: if this variable is at its implied bound, |
847 | // then the other variables must be at their bounds, |
848 | // which means the bounds will stop them even if the aren't basic. |
849 | if (rowstat[irow] & PRESOLVE_XBASIC) |
850 | rowstat[irow] = 0; |
851 | else { |
852 | int k; |
853 | for (k=0; k<ninrow; k++) { |
854 | int col = rowcols[k]; |
855 | if (cdone[col] && |
856 | (colstat[col] & PRESOLVE_XBASIC) && |
857 | ((fabs(clo[col] - sol[col]) <= ztolzb && rcosts[col] >= -ztoldj) || |
858 | (fabs(cup[col] - sol[col]) <= ztolzb && rcosts[col] <= ztoldj))) |
859 | break; |
860 | } |
861 | if (k<ninrow) { |
862 | int col = rowcols[k]; |
863 | // steal this basic variable |
864 | #if PRESOLVE_DEBUG |
865 | printf("PIVOTING ON COL: %d %d -> %d\n" , irow, col, jcol); |
866 | #endif |
867 | colstat[col] = 0; |
868 | |
869 | // since all vars were at their bounds, the slack must be 0 |
870 | PRESOLVEASSERT(fabs(acts[irow]) < ZTOLDP); |
871 | rowstat[irow] = PRESOLVE_XBASIC; |
872 | } |
873 | else { |
874 | // should never happen? |
875 | abort(); |
876 | } |
877 | |
878 | // get rid of any remaining basic structurals, since their rcosts |
879 | // are going to become non-zero in a second. |
880 | abort(); |
881 | ///////////////////zero_pivot(); |
882 | } |
883 | |
884 | double rdual_adjust; |
885 | { |
886 | CoinBigIndex kk = presolve_find_row(irow, mcstrt[jcol], mcstrt[jcol] + hincol[jcol], hrow); |
887 | // adjust rowdual to cancel out reduced cost |
888 | // should probably search for col with largest factor |
889 | rdual_adjust = (rcosts[jcol] / colels[kk]); |
890 | rowduals[irow] += rdual_adjust; |
891 | colstat[jcol] = PRESOLVE_XBASIC; |
892 | } |
893 | |
894 | for (k=0; k<ninrow; k++) { |
895 | int jcol = rowcols[k]; |
896 | CoinBigIndex kk = presolve_find_row(irow, mcstrt[jcol], mcstrt[jcol] + hincol[jcol], hrow); |
897 | |
898 | rcosts[jcol] -= (rdual_adjust * colels[kk]); |
899 | } |
900 | |
901 | { |
902 | int k; |
903 | int badbasic = -1; |
904 | |
905 | // we may have just screwed up the rcost of another basic variable |
906 | for (k=0; k<ninrow; k++) { |
907 | int col = rowcols[k]; |
908 | if (col != jcol && |
909 | cdone[col] && |
910 | (colstat[col] & PRESOLVE_XBASIC) && |
911 | !(fabs(rcosts[col]) < ztoldj)) |
912 | if (badbasic == -1) |
913 | badbasic = k; |
914 | else |
915 | abort(); // two!! what to do??? |
916 | } |
917 | |
918 | if (badbasic != -1) { |
919 | int col = rowcols[badbasic]; |
920 | |
921 | if (fabs(acts[irow]) < ZTOLDP) { |
922 | #if PRESOLVE_DEBUG |
923 | printf("PIVOTING COL TO SLACK!: %d %d\n" , irow, col); |
924 | #endif |
925 | colstat[col] = 0; |
926 | rowstat[irow] = PRESOLVE_XBASIC; |
927 | } |
928 | else |
929 | abort(); |
930 | } |
931 | } |
932 | } |
933 | } |
934 | #endif // #if 0 // (B) |
935 | #endif // #if 0 // (A) |
936 | |
937 | forcing_constraint_action::~forcing_constraint_action() |
938 | { |
939 | int i; |
940 | for (i=0;i<nactions_;i++) { |
941 | //delete [] actions_[i].rowcols; MS Visual C++ V6 can not compile |
942 | //delete [] actions_[i].bounds; MS Visual C++ V6 can not compile |
943 | deleteAction(actions_[i].rowcols,int *); |
944 | deleteAction(actions_[i].bounds,double *); |
945 | } |
946 | // delete [] actions_; MS Visual C++ V6 can not compile |
947 | deleteAction(actions_,action *); |
948 | } |
949 | |