1/* $Id: CoinPresolveImpliedFree.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 "CoinPresolveSubst.hpp"
11#include "CoinPresolveIsolated.hpp"
12#include "CoinPresolveFixed.hpp"
13#include "CoinPresolveImpliedFree.hpp"
14#include "CoinPresolveUseless.hpp"
15#include "CoinPresolveForcing.hpp"
16#include "CoinMessage.hpp"
17#include "CoinHelperFunctions.hpp"
18#include "CoinSort.hpp"
19#include "CoinFinite.hpp"
20
21#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
22#include "CoinPresolvePsdebug.hpp"
23#endif
24
25namespace { // begin unnamed file-local namespace
26
27const CoinPresolveAction *testRedundant (CoinPresolveMatrix *prob,
28 const CoinPresolveAction *next,
29 int & numberInfeasible)
30{
31 numberInfeasible=0;
32 int numberColumns = prob->ncols_;
33 double * columnLower = new double[numberColumns];
34 double * columnUpper = new double[numberColumns];
35 CoinMemcpyN(prob->clo_,numberColumns,columnLower);
36 CoinMemcpyN(prob->cup_,numberColumns,columnUpper);
37
38 const double *element = prob->rowels_;
39 const int *column = prob->hcol_;
40 const CoinBigIndex *rowStart = prob->mrstrt_;
41 const int *rowLength = prob->hinrow_;
42 int numberRows = prob->nrows_;
43 const int *hrow = prob->hrow_;
44 const CoinBigIndex *mcstrt = prob->mcstrt_;
45 const int *hincol = prob->hincol_;
46
47 int *useless_rows = prob->usefulRowInt_+numberRows; //new int[numberRows];
48 int nuseless_rows = 0;
49
50 double *rowLower = prob->rlo_;
51 double *rowUpper = prob->rup_;
52
53 double tolerance = prob->feasibilityTolerance_;
54 int numberChanged=1,iPass=0;
55#define USE_SMALL_LARGE
56#ifdef USE_SMALL_LARGE
57 double large = 1.0e15; // treat bounds > this as infinite
58#else
59 double large = 1.0e20; // treat bounds > this as infinite
60#endif
61#ifndef NDEBUG
62 double large2= 1.0e10*large;
63#endif
64 int totalTightened = 0;
65
66 int iRow, iColumn;
67
68 char * markRow = reinterpret_cast<char *>(prob->usefulRowInt_); // wasnew int[numberRows];
69 for (iRow=0;iRow<numberRows;iRow++) {
70 if ((rowLower[iRow]>-large||rowUpper[iRow]<large)&&rowLength[iRow]>0) {
71 markRow[iRow]=-1;
72 } else {
73 markRow[iRow]=1;
74 if (rowLength[iRow]>0) {
75 // Row is redundant
76 useless_rows[nuseless_rows++] = iRow;
77 prob->addRow(iRow);
78 }
79 }
80 }
81#define MAXPASS 10
82 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
83 double relaxedTolerance = 100.0*tolerance;
84
85 // Loop round seeing if we can tighten bounds
86 // Would be faster to have a stack of possible rows
87 // and we put altered rows back on stack
88 int numberCheck=-1;
89 while(numberChanged>numberCheck) {
90
91 numberChanged = 0; // Bounds tightened this pass
92
93 if (iPass==MAXPASS) break;
94 iPass++;
95
96 for (iRow = 0; iRow < numberRows; iRow++) {
97
98 if (markRow[iRow]==-1) {
99
100 // possible row - but mark as useless next time
101 markRow[iRow]=-2;
102 int infiniteUpper = 0;
103 int infiniteLower = 0;
104 double maximumUp = 0.0;
105 double maximumDown = 0.0;
106 double newBound;
107 CoinBigIndex rStart = rowStart[iRow];
108 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
109 CoinBigIndex j;
110 // Compute possible lower and upper ranges
111
112 for (j = rStart; j < rEnd; ++j) {
113 double value=element[j];
114 iColumn = column[j];
115 if (value > 0.0) {
116 if (columnUpper[iColumn] < large)
117 maximumUp += columnUpper[iColumn] * value;
118 else
119 ++infiniteUpper;
120 if (columnLower[iColumn] > -large)
121 maximumDown += columnLower[iColumn] * value;
122 else
123 ++infiniteLower;
124 } else if (value<0.0) {
125 if (columnUpper[iColumn] < large)
126 maximumDown += columnUpper[iColumn] * value;
127 else
128 ++infiniteLower;
129 if (columnLower[iColumn] > -large)
130 maximumUp += columnLower[iColumn] * value;
131 else
132 ++infiniteUpper;
133 }
134 }
135 // Build in a margin of error
136 maximumUp += 1.0e-8*fabs(maximumUp);
137 maximumDown -= 1.0e-8*fabs(maximumDown);
138 double maxUp = maximumUp+infiniteUpper*1.0e31;
139 double maxDown = maximumDown-infiniteLower*1.0e31;
140 if (maxUp <= rowUpper[iRow] + tolerance &&
141 maxDown >= rowLower[iRow] - tolerance) {
142
143 } else {
144 if (maxUp < rowLower[iRow] -relaxedTolerance ||
145 maxDown > rowUpper[iRow]+relaxedTolerance) {
146 if(!fixInfeasibility) {
147 // problem is infeasible - exit at once
148 numberInfeasible++;
149 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
150 prob->messages())
151 <<iRow
152 <<rowLower[iRow]
153 <<rowUpper[iRow]
154 <<CoinMessageEol;
155 break;
156 } else {
157 continue;
158 }
159 }
160 double lower = rowLower[iRow];
161 double upper = rowUpper[iRow];
162 // Clean up
163 if (maximumUp < lower && maximumUp > lower -relaxedTolerance)
164 maximumUp=lower;
165 if (maximumDown > upper && maximumDown < upper +relaxedTolerance)
166 maximumDown=upper;
167 for (j = rStart; j < rEnd; ++j) {
168 double value=element[j];
169 iColumn = column[j];
170 double nowLower = columnLower[iColumn];
171 double nowUpper = columnUpper[iColumn];
172 if (value > 0.0) {
173 // positive value
174 if (lower>-large) {
175 if (!infiniteUpper) {
176 assert(nowUpper < large2);
177 newBound = nowUpper +
178 (lower - maximumUp) / value;
179 // relax if original was large
180 if (fabs(maximumUp)>1.0e8)
181 newBound -= 1.0e-12*fabs(maximumUp);
182 } else if (infiniteUpper==1&&nowUpper>=large) {
183 newBound = (lower -maximumUp) / value;
184 // relax if original was large
185 if (fabs(maximumUp)>1.0e8)
186 newBound -= 1.0e-12*fabs(maximumUp);
187 } else {
188 newBound = -COIN_DBL_MAX;
189 }
190 if (newBound > nowLower + 1.0e-12&&newBound>-large) {
191 // Tighten the lower bound
192 columnLower[iColumn] = newBound;
193 markRow[iRow]=1;
194 numberChanged++;
195 // Mark rows as possible
196 CoinBigIndex kcs = mcstrt[iColumn];
197 CoinBigIndex kce = kcs + hincol[iColumn];
198 CoinBigIndex k;
199 for (k=kcs; k<kce; ++k) {
200 int row = hrow[k];
201 if (markRow[row]==-2) {
202 // on list for next time
203 markRow[row]=-1;
204 }
205 }
206 // check infeasible (relaxed)
207 if (nowUpper - newBound <
208 -relaxedTolerance) {
209 numberInfeasible++;
210 }
211 // adjust
212 double now;
213 if (nowLower<=-large) {
214 now=0.0;
215 infiniteLower--;
216 } else {
217 now = nowLower;
218 }
219 maximumDown += (newBound-now) * value;
220 nowLower = newBound;
221 //#define FREE_DEBUG 2
222#if FREE_DEBUG >1
223 if (fabs((newBound-now)*value)>1.0e8) {
224 // recompute
225 infiniteLower = 0;
226 maximumDown = 0.0;
227 CoinBigIndex j2;
228 // Compute possible lower and upper ranges
229 for (j2 = rStart; j2 < rEnd; ++j2) {
230 double value=element[j2];
231 int iColumn = column[j2];
232 if (value > 0.0) {
233 if (columnLower[iColumn] > -large)
234 maximumDown += columnLower[iColumn] * value;
235 else
236 ++infiniteLower;
237 } else if (value<0.0) {
238 if (columnUpper[iColumn] < large)
239 maximumDown += columnUpper[iColumn] * value;
240 else
241 ++infiniteLower;
242 }
243 }
244 // Build in a margin of error
245 maximumDown -= 1.0e-8*fabs(maximumDown);
246 if (maximumDown > upper && maximumDown < upper +relaxedTolerance)
247 maximumDown=upper;
248 }
249#endif
250#if FREE_DEBUG
251#define DEBUG_TOLERANCE 1.0e-10
252 { // DEBUG
253 int infiniteUpper2 = 0;
254 int infiniteLower2 = 0;
255 double maximumUp2 = 0.0;
256 double maximumDown2 = 0.0;
257 CoinBigIndex j2;
258 // Compute possible lower and upper ranges
259 for (j2 = rStart; j2 < rEnd; ++j2) {
260 double value=element[j2];
261 int iColumn = column[j2];
262 if (value > 0.0) {
263 if (columnUpper[iColumn] < large)
264 maximumUp2 += columnUpper[iColumn] * value;
265 else
266 ++infiniteUpper2;
267 if (columnLower[iColumn] > -large)
268 maximumDown2 += columnLower[iColumn] * value;
269 else
270 ++infiniteLower2;
271 } else if (value<0.0) {
272 if (columnUpper[iColumn] < large)
273 maximumDown2 += columnUpper[iColumn] * value;
274 else
275 ++infiniteLower2;
276 if (columnLower[iColumn] > -large)
277 maximumUp2 += columnLower[iColumn] * value;
278 else
279 ++infiniteUpper2;
280 }
281 }
282 // Build in a margin of error
283 maximumUp2 += 1.0e-8*fabs(maximumUp2);
284 maximumDown2 -= 1.0e-8*fabs(maximumDown2);
285 if (maximumUp2 < lower && maximumUp2 > lower -relaxedTolerance)
286 maximumUp2=lower;
287 if (maximumDown2 > upper && maximumDown2 < upper +relaxedTolerance)
288 maximumDown2=upper;
289 assert (infiniteLower==infiniteLower2);
290 assert (infiniteUpper==infiniteUpper2);
291 assert (fabs(maximumDown-maximumDown2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumDown)));
292 assert (fabs(maximumUp-maximumUp2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumUp)));
293 } // END DEBUG
294#endif
295 }
296 }
297 if (upper <large) {
298 if (!infiniteLower) {
299 assert(nowLower >- large2);
300 newBound = nowLower +
301 (upper - maximumDown) / value;
302 // relax if original was large
303 if (fabs(maximumDown)>1.0e8)
304 newBound += 1.0e-12*fabs(maximumDown);
305 } else if (infiniteLower==1&&nowLower<=-large) {
306 newBound = (upper - maximumDown) / value;
307 // relax if original was large
308 if (fabs(maximumDown)>1.0e8)
309 newBound += 1.0e-12*fabs(maximumDown);
310 } else {
311 newBound = COIN_DBL_MAX;
312 }
313 if (newBound < nowUpper - 1.0e-12&&newBound<large) {
314 // Tighten the upper bound
315 columnUpper[iColumn] = newBound;
316 markRow[iRow]=1;
317 numberChanged++;
318 // Mark rows as possible
319 CoinBigIndex kcs = mcstrt[iColumn];
320 CoinBigIndex kce = kcs + hincol[iColumn];
321 CoinBigIndex k;
322 for (k=kcs; k<kce; ++k) {
323 int row = hrow[k];
324 if (markRow[row]==-2) {
325 // on list for next time
326 markRow[row]=-1;
327 }
328 }
329 // check infeasible (relaxed)
330 if (newBound - nowLower <
331 -relaxedTolerance) {
332 numberInfeasible++;
333 }
334 // adjust
335 double now;
336 if (nowUpper>=large) {
337 now=0.0;
338 infiniteUpper--;
339 } else {
340 now = nowUpper;
341 }
342 maximumUp += (newBound-now) * value;
343 nowUpper = newBound;
344#if FREE_DEBUG >1
345 if (fabs((newBound-now)*value)>1.0e8) {
346 // recompute
347 infiniteUpper = 0;
348 maximumUp = 0.0;
349 CoinBigIndex j2;
350 // Compute possible lower and upper ranges
351 for (j2 = rStart; j2 < rEnd; ++j2) {
352 double value=element[j2];
353 int iColumn = column[j2];
354 if (value > 0.0) {
355 if (columnUpper[iColumn] < large)
356 maximumUp += columnUpper[iColumn] * value;
357 else
358 ++infiniteUpper;
359 } else if (value<0.0) {
360 if (columnLower[iColumn] > -large)
361 maximumUp += columnLower[iColumn] * value;
362 else
363 ++infiniteUpper;
364 }
365 }
366 // Build in a margin of error
367 maximumUp += 1.0e-8*fabs(maximumUp);
368 if (maximumUp < lower && maximumUp > lower -relaxedTolerance)
369 maximumUp=lower;
370 }
371#endif
372#if FREE_DEBUG
373 { // DEBUG
374 int infiniteUpper2 = 0;
375 int infiniteLower2 = 0;
376 double maximumUp2 = 0.0;
377 double maximumDown2 = 0.0;
378 CoinBigIndex j2;
379 // Compute possible lower and upper ranges
380 for (j2 = rStart; j2 < rEnd; ++j2) {
381 double value=element[j2];
382 int iColumn = column[j2];
383 if (value > 0.0) {
384 if (columnUpper[iColumn] < large)
385 maximumUp2 += columnUpper[iColumn] * value;
386 else
387 ++infiniteUpper2;
388 if (columnLower[iColumn] > -large)
389 maximumDown2 += columnLower[iColumn] * value;
390 else
391 ++infiniteLower2;
392 } else if (value<0.0) {
393 if (columnUpper[iColumn] < large)
394 maximumDown2 += columnUpper[iColumn] * value;
395 else
396 ++infiniteLower2;
397 if (columnLower[iColumn] > -large)
398 maximumUp2 += columnLower[iColumn] * value;
399 else
400 ++infiniteUpper2;
401 }
402 }
403 // Build in a margin of error
404 maximumUp2 += 1.0e-8*fabs(maximumUp2);
405 maximumDown2 -= 1.0e-8*fabs(maximumDown2);
406 if (maximumUp2 < lower && maximumUp2 > lower -relaxedTolerance)
407 maximumUp2=lower;
408 if (maximumDown2 > upper && maximumDown2 < upper +relaxedTolerance)
409 maximumDown2=upper;
410 assert (infiniteLower==infiniteLower2);
411 assert (infiniteUpper==infiniteUpper2);
412 assert (fabs(maximumDown-maximumDown2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumDown)));
413 assert (fabs(maximumUp-maximumUp2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumUp)));
414 } // END DEBUG
415#endif
416 }
417 }
418 } else {
419 // negative value
420 if (lower>-large) {
421 if (!infiniteUpper) {
422 assert(nowLower < large2);
423 newBound = nowLower +
424 (lower - maximumUp) / value;
425 // relax if original was large
426 if (fabs(maximumUp)>1.0e8)
427 newBound += 1.0e-12*fabs(maximumUp);
428 } else if (infiniteUpper==1&&nowLower<=-large) {
429 newBound = (lower -maximumUp) / value;
430 // relax if original was large
431 if (fabs(maximumUp)>1.0e8)
432 newBound += 1.0e-12*fabs(maximumUp);
433 } else {
434 newBound = COIN_DBL_MAX;
435 }
436 if (newBound < nowUpper - 1.0e-12&&newBound<large) {
437 // Tighten the upper bound
438 columnUpper[iColumn] = newBound;
439 markRow[iRow]=1;
440 numberChanged++;
441 // Mark rows as possible
442 CoinBigIndex kcs = mcstrt[iColumn];
443 CoinBigIndex kce = kcs + hincol[iColumn];
444 CoinBigIndex k;
445 for (k=kcs; k<kce; ++k) {
446 int row = hrow[k];
447 if (markRow[row]==-2) {
448 // on list for next time
449 markRow[row]=-1;
450 }
451 }
452 // check infeasible (relaxed)
453 if (newBound - nowLower <
454 -relaxedTolerance) {
455 numberInfeasible++;
456 }
457 // adjust
458 double now;
459 if (nowUpper>=large) {
460 now=0.0;
461 infiniteLower--;
462 } else {
463 now = nowUpper;
464 }
465 maximumDown += (newBound-now) * value;
466 nowUpper = newBound;
467#if FREE_DEBUG >1
468 if (fabs((newBound-now)*value)>1.0e8) {
469 // recompute
470 infiniteLower = 0;
471 maximumDown = 0.0;
472 CoinBigIndex j2;
473 // Compute possible lower and upper ranges
474 for (j2 = rStart; j2 < rEnd; ++j2) {
475 double value=element[j2];
476 int iColumn = column[j2];
477 if (value > 0.0) {
478 if (columnLower[iColumn] > -large)
479 maximumDown += columnLower[iColumn] * value;
480 else
481 ++infiniteLower;
482 } else if (value<0.0) {
483 if (columnUpper[iColumn] < large)
484 maximumDown += columnUpper[iColumn] * value;
485 else
486 ++infiniteLower;
487 }
488 }
489 // Build in a margin of error
490 maximumDown -= 1.0e-8*fabs(maximumDown);
491 if (maximumDown > upper && maximumDown < upper +relaxedTolerance)
492 maximumDown=upper;
493 }
494#endif
495#if FREE_DEBUG
496 { // DEBUG
497 int infiniteUpper2 = 0;
498 int infiniteLower2 = 0;
499 double maximumUp2 = 0.0;
500 double maximumDown2 = 0.0;
501 CoinBigIndex j2;
502 // Compute possible lower and upper ranges
503 for (j2 = rStart; j2 < rEnd; ++j2) {
504 double value=element[j2];
505 int iColumn = column[j2];
506 if (value > 0.0) {
507 if (columnUpper[iColumn] < large)
508 maximumUp2 += columnUpper[iColumn] * value;
509 else
510 ++infiniteUpper2;
511 if (columnLower[iColumn] > -large)
512 maximumDown2 += columnLower[iColumn] * value;
513 else
514 ++infiniteLower2;
515 } else if (value<0.0) {
516 if (columnUpper[iColumn] < large)
517 maximumDown2 += columnUpper[iColumn] * value;
518 else
519 ++infiniteLower2;
520 if (columnLower[iColumn] > -large)
521 maximumUp2 += columnLower[iColumn] * value;
522 else
523 ++infiniteUpper2;
524 }
525 }
526 // Build in a margin of error
527 maximumUp2 += 1.0e-8*fabs(maximumUp2);
528 maximumDown2 -= 1.0e-8*fabs(maximumDown2);
529 if (maximumUp2 < lower && maximumUp2 > lower -relaxedTolerance)
530 maximumUp2=lower;
531 if (maximumDown2 > upper && maximumDown2 < upper +relaxedTolerance)
532 maximumDown2=upper;
533 assert (infiniteLower==infiniteLower2);
534 assert (infiniteUpper==infiniteUpper2);
535 assert (fabs(maximumDown-maximumDown2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumDown)));
536 assert (fabs(maximumUp-maximumUp2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumUp)));
537 } // END DEBUG
538#endif
539 }
540 }
541 if (upper <large) {
542 if (!infiniteLower) {
543 assert(nowUpper < large2);
544 newBound = nowUpper +
545 (upper - maximumDown) / value;
546 // relax if original was large
547 if (fabs(maximumDown)>1.0e8)
548 newBound -= 1.0e-12*fabs(maximumDown);
549 } else if (infiniteLower==1&&nowUpper>=large) {
550 newBound = (upper - maximumDown) / value;
551 // relax if original was large
552 if (fabs(maximumDown)>1.0e8)
553 newBound -= 1.0e-12*fabs(maximumDown);
554 } else {
555 newBound = -COIN_DBL_MAX;
556 }
557 if (newBound > nowLower + 1.0e-12&&newBound>-large) {
558 // Tighten the lower bound
559 columnLower[iColumn] = newBound;
560 markRow[iRow]=1;
561 numberChanged++;
562 // Mark rows as possible
563 CoinBigIndex kcs = mcstrt[iColumn];
564 CoinBigIndex kce = kcs + hincol[iColumn];
565 CoinBigIndex k;
566 for (k=kcs; k<kce; ++k) {
567 int row = hrow[k];
568 if (markRow[row]==-2) {
569 // on list for next time
570 markRow[row]=-1;
571 }
572 }
573 // check infeasible (relaxed)
574 if (nowUpper - newBound <
575 -relaxedTolerance) {
576 numberInfeasible++;
577 }
578 // adjust
579 double now;
580 if (nowLower<=-large) {
581 now=0.0;
582 infiniteUpper--;
583 } else {
584 now = nowLower;
585 }
586 maximumUp += (newBound-now) * value;
587 nowLower = newBound;
588#if FREE_DEBUG >1
589 if (fabs((newBound-now)*value)>1.0e8) {
590 // recompute
591 infiniteUpper = 0;
592 maximumUp = 0.0;
593 CoinBigIndex j2;
594 // Compute possible lower and upper ranges
595 for (j2 = rStart; j2 < rEnd; ++j2) {
596 double value=element[j2];
597 int iColumn = column[j2];
598 if (value > 0.0) {
599 if (columnUpper[iColumn] < large)
600 maximumUp += columnUpper[iColumn] * value;
601 else
602 ++infiniteUpper;
603 } else if (value<0.0) {
604 if (columnLower[iColumn] > -large)
605 maximumUp += columnLower[iColumn] * value;
606 else
607 ++infiniteUpper;
608 }
609 }
610 // Build in a margin of error
611 maximumUp += 1.0e-8*fabs(maximumUp);
612 if (maximumUp < lower && maximumUp > lower -relaxedTolerance)
613 maximumUp=lower;
614 }
615#endif
616#if FREE_DEBUG
617 { // DEBUG
618 int infiniteUpper2 = 0;
619 int infiniteLower2 = 0;
620 double maximumUp2 = 0.0;
621 double maximumDown2 = 0.0;
622 CoinBigIndex j2;
623 // Compute possible lower and upper ranges
624 for (j2 = rStart; j2 < rEnd; ++j2) {
625 double value=element[j2];
626 int iColumn = column[j2];
627 if (value > 0.0) {
628 if (columnUpper[iColumn] < large)
629 maximumUp2 += columnUpper[iColumn] * value;
630 else
631 ++infiniteUpper2;
632 if (columnLower[iColumn] > -large)
633 maximumDown2 += columnLower[iColumn] * value;
634 else
635 ++infiniteLower2;
636 } else if (value<0.0) {
637 if (columnUpper[iColumn] < large)
638 maximumDown2 += columnUpper[iColumn] * value;
639 else
640 ++infiniteLower2;
641 if (columnLower[iColumn] > -large)
642 maximumUp2 += columnLower[iColumn] * value;
643 else
644 ++infiniteUpper2;
645 }
646 }
647 // Build in a margin of error
648 maximumUp2 += 1.0e-8*fabs(maximumUp2);
649 maximumDown2 -= 1.0e-8*fabs(maximumDown2);
650 if (maximumUp2 < lower && maximumUp2 > lower -relaxedTolerance)
651 maximumUp2=lower;
652 if (maximumDown2 > upper && maximumDown2 < upper +relaxedTolerance)
653 maximumDown2=upper;
654 assert (infiniteLower==infiniteLower2);
655 assert (infiniteUpper==infiniteUpper2);
656 assert (fabs(maximumDown-maximumDown2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumDown)));
657 assert (fabs(maximumUp-maximumUp2)<DEBUG_TOLERANCE*(1.0e6+fabs(maximumUp)));
658 } // END DEBUG
659#endif
660 }
661 }
662 }
663 }
664 }
665 }
666 }
667 totalTightened += numberChanged;
668 if (iPass==1)
669 numberCheck=CoinMax(10,numberChanged>>5);
670 if (numberInfeasible) break;
671 }
672 if (!numberInfeasible) {
673 //#define NO_FORCING
674#ifdef NO_FORCING
675#define ZZ 1
676#ifdef ZZ
677 double * clo2 = CoinCopyOfArray(prob->clo_,numberColumns);
678 double * cup2 = CoinCopyOfArray(prob->cup_,numberColumns);
679#endif
680 forcing_constraint_action::action * actions = NULL;
681 int numActions=0;
682 int maxActions=0;
683 double * clo = prob->clo_;
684 double * cup = prob->cup_;
685 double *csol = prob->sol_ ;
686 int * fixed = prob->usefulColumnInt_;
687 for (int i=0;i<numberColumns;i++) {
688 if (cup[i]>clo[i]+tolerance)
689 fixed[i]=0;
690 else
691 fixed[i]=-1;
692 }
693#endif
694 for (iRow = 0; iRow < numberRows; iRow++) {
695
696 if (markRow[iRow]<0) {
697
698 // possible row
699 int infiniteUpper = 0;
700 int infiniteLower = 0;
701 double maximumUp = 0.0;
702 double maximumDown = 0.0;
703 CoinBigIndex rStart = rowStart[iRow];
704 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
705 CoinBigIndex j;
706 // Compute possible lower and upper ranges
707
708 for (j = rStart; j < rEnd; ++j) {
709 double value=element[j];
710 iColumn = column[j];
711 if (value > 0.0) {
712 if (columnUpper[iColumn] < large)
713 maximumUp += columnUpper[iColumn] * value;
714 else
715 ++infiniteUpper;
716 if (columnLower[iColumn] > -large)
717 maximumDown += columnLower[iColumn] * value;
718 else
719 ++infiniteLower;
720 } else if (value<0.0) {
721 if (columnUpper[iColumn] < large)
722 maximumDown += columnUpper[iColumn] * value;
723 else
724 ++infiniteLower;
725 if (columnLower[iColumn] > -large)
726 maximumUp += columnLower[iColumn] * value;
727 else
728 ++infiniteUpper;
729 }
730 }
731 // Build in a margin of error
732 maximumUp += 1.0e-8*fabs(maximumUp);
733 maximumDown -= 1.0e-8*fabs(maximumDown);
734 double maxUp = maximumUp+infiniteUpper*1.0e31;
735 double maxDown = maximumDown-infiniteLower*1.0e31;
736 if (maxUp <= rowUpper[iRow] + tolerance &&
737 maxDown >= rowLower[iRow] - tolerance) {
738#ifndef NO_FORCING
739 // Row is redundant
740 useless_rows[nuseless_rows++] = iRow;
741 prob->addRow(iRow);
742#else
743 if (maxUp <= rowUpper[iRow] - tolerance &&
744 maxDown >= rowLower[iRow] + tolerance) {
745 // Row is redundant
746 useless_rows[nuseless_rows++] = iRow;
747 prob->addRow(iRow);
748 } else {
749 CoinBigIndex k;
750 double tol2 = 10.0*tolerance;
751 int nFree=0;
752 bool bad=false;
753 for ( k=rStart; k<rEnd; k++) {
754 int jcol = column[k];
755 if (columnUpper[jcol]-columnLower[jcol]>tol2)
756 nFree++;
757 if (fixed[jcol]>0)
758 bad=true;
759 }
760 if (nFree&&!bad) {
761 // the lower bound can just be reached, or
762 // the upper bound can just be reached;
763 // called a "forcing constraint" in the paper (p. 226)
764 const int lbound_tight = (maxUp < PRESOLVE_INF &&
765 fabs(rowLower[iRow] - maxUp) < tolerance);
766 double *bounds = new double[rowLength[iRow]];
767 int *rowcols = new int[rowLength[iRow]];
768 int lk = rStart; // load fix-to-down in front
769 int uk = rEnd; // load fix-to-up in back
770 for ( k=rStart; k<rEnd; k++) {
771 int jcol = column[k];
772 if (!fixed[jcol])
773 fixed[jcol]=1;
774 prob->addCol(jcol);
775 double coeff = element[k];
776 // one of the two contributed to maxup - set the other to that
777 if (lbound_tight == (coeff > 0.0)) {
778 --uk;
779 bounds[uk-rStart] = clo[jcol];
780 rowcols[uk-rStart] = jcol;
781 if (csol != 0) {
782 csol[jcol] = columnUpper[jcol] ;
783 }
784#ifdef ZZ
785 assert (columnLower[jcol]==clo[jcol]);
786 clo[jcol] = columnUpper[jcol];
787#endif
788 columnLower[jcol] = columnUpper[jcol];
789 } else {
790 bounds[lk-rStart] = cup[jcol];
791 rowcols[lk-rStart] = jcol;
792 ++lk;
793 if (csol != 0) {
794 csol[jcol] = columnLower[jcol] ;
795 }
796#ifdef ZZ
797 assert (columnUpper[jcol]==cup[jcol]);
798 cup[jcol] = columnLower[jcol];
799#endif
800 columnUpper[jcol] = columnLower[jcol];
801 }
802 }
803 PRESOLVEASSERT(uk == lk);
804 if (numActions==maxActions) {
805 maxActions = (11*maxActions)/10+10;
806 forcing_constraint_action::action * temp =
807 new forcing_constraint_action::action [maxActions];
808 for (int i=0;i<numActions;i++)
809 temp[i]=actions[i];
810 delete [] actions;
811 actions = temp;
812 }
813
814 forcing_constraint_action::action * f = &actions[numActions];
815 numActions++;
816 PRESOLVE_DETAIL_PRINT(printf("pre_impliedfree %dR E\n",iRow));
817 f->row = iRow;
818 f->nlo = lk-rStart;
819 f->nup = rEnd-uk;
820 f->rowcols = rowcols;
821 f->bounds = bounds;
822 } else if (!nFree&&!bad) {
823 // Row is redundant
824 useless_rows[nuseless_rows++] = iRow;
825 prob->addRow(iRow);
826 }
827 }
828#endif
829 }
830 }
831 }
832#ifdef NO_FORCING
833 if (numActions) {
834#if PRESOLVE_SUMMARY
835 printf("NFORCED_FROM_IMPLIED: %d\n", numActions);
836#endif
837 next = new forcing_constraint_action(numActions,
838 CoinCopyOfArray(actions,numActions), next);
839 deleteAction(actions,action*);
840 actions=NULL;
841 }
842 assert (!actions);
843#endif
844 if (nuseless_rows) {
845 next = useless_constraint_action::presolve(prob,
846 useless_rows, nuseless_rows,
847 next);
848 }
849#ifdef NO_FORCING
850#ifdef ZZ
851 // temp
852 for (int i=0;i<numberColumns;i++) {
853 clo[i]=clo2[i];
854 cup[i]=cup2[i];
855 }
856#endif
857 //delete [] clo2;
858 //delete [] cup2;
859 if (numActions) {
860 //int * fixed = prob->usefulColumnInt_;
861 // See if any columns fixed
862 int nFixed=0;
863 for (int i=0;i<numberColumns;i++) {
864 //if (columnUpper[i]<columnLower[i]+tolerance&&
865 // cup[i]>clo[i]+tolerance) {
866 if (fixed[i]>0) {
867 fixed[nFixed++]=i;
868 assert (columnUpper[i]>columnLower[i]-10.0*tolerance);
869 }
870 }
871 //if (nFixed)
872 //printf("could fix %d\n",nFixed);
873 //next = remove_fixed_action::presolve(prob,fixed,nFixed,next) ;
874 }
875#endif
876 if ((prob->presolveOptions_&16)!=0) {
877 // may not unroll
878 const unsigned char *integerType = prob->integerType_;
879 double *csol = prob->sol_ ;
880 double * clo = prob->clo_;
881 double * cup = prob->cup_;
882 int * fixed = prob->usefulColumnInt_;
883 int nFixed=0;
884 int nChanged=0;
885 for (int i=0;i<numberColumns;i++) {
886 if (clo[i]==cup[i])
887 continue;
888 double lower = columnLower[i];
889 double upper = columnUpper[i];
890 if (integerType[i]) {
891 //if (floor(upper+1.0e-4)<upper)
892 upper=floor(upper+1.0e-4);
893 //if (ceil(lower-1.0e-4)>lower)
894 lower=ceil(lower-1.0e-4);
895 }
896 if (upper-lower<1.0e-8) {
897 if (upper-lower<-tolerance)
898 numberInfeasible++;
899 if (CoinMin(fabs(upper),fabs(lower))<=1.0e-7)
900 upper=0.0;
901 fixed[nFixed++]=i;
902 //printf("fixing %d to %g\n",i,upper);
903 prob->addCol(i);
904 cup[i]=upper;
905 clo[i]=upper;
906 if (csol != 0)
907 csol[i] = upper;
908 } else {
909 if (integerType[i]) {
910 if (upper<cup[i]) {
911 cup[i]=upper;
912 nChanged++;
913 prob->addCol(i);
914 }
915 if (lower>clo[i]) {
916 clo[i]=lower;
917 nChanged++;
918 prob->addCol(i);
919 }
920 }
921 }
922 }
923#ifdef CLP_INVESTIGATE
924 if (nFixed||nChanged)
925 printf("%d fixed in impliedfree, %d changed\n",nFixed,nChanged);
926#endif
927 if (nFixed)
928 next = remove_fixed_action::presolve(prob,fixed,nFixed,next) ;
929 }
930 }
931
932 delete [] columnLower;
933 delete [] columnUpper;
934 return next;
935}
936
937} // end unnamed file-local namespace
938
939// If there is a row with a singleton column such that no matter what
940// the values of the other variables are, the constraint forces the singleton
941// column to have a feasible value, then we can drop the column and row,
942// since we just compute the value of the column from the row in postsolve.
943// This seems vaguely similar to the case of a useless constraint, but it
944// is different. For example, if the singleton column is already free,
945// then this operation will eliminate it, but the constraint is not useless
946// (assuming the constraint is not trivial), since the variables do not imply an
947// upper or lower bound.
948//
949// If the column is not singleton, we can still do something similar if the
950// constraint is an equality constraint. In that case, we substitute away
951// the variable in the other constraints it appears in. This introduces
952// new coefficients, but the total number of coefficients never increases
953// if the column has only two constraints, and may not increase much even
954// if there are more.
955//
956// There is nothing to prevent us from substituting away a variable
957// in an equality from the other constraints it appears in, but since
958// that causes fill-in, it wouldn't make sense unless we could then
959// drop the equality itself. We can't do that if the bounds on the
960// variable in equation aren't implied by the equality.
961// Another way of thinking of this is that there is nothing special
962// about an equality; just like one can't always drop an inequality constraint
963// with a column singleton, one can't always drop an equality.
964//
965// It is possible for two singleton columns to be in the same row.
966// In that case, the other one will become empty. If its bounds and
967// costs aren't just right, this signals an unbounded problem.
968// We don't need to check that specially here.
969//
970// invariant: loosely packed
971const CoinPresolveAction *implied_free_action::presolve(CoinPresolveMatrix *prob,
972 const CoinPresolveAction *next,
973 int & fill_level)
974{
975 double startTime = 0.0;
976 int startEmptyRows=0;
977 int startEmptyColumns = 0;
978 if (prob->tuning_) {
979 startTime = CoinCpuTime();
980 startEmptyRows = prob->countEmptyRows();
981 startEmptyColumns = prob->countEmptyCols();
982 }
983 double *colels = prob->colels_;
984 int *hrow = prob->hrow_;
985 const CoinBigIndex *mcstrt = prob->mcstrt_;
986 int *hincol = prob->hincol_;
987 const int ncols = prob->ncols_;
988
989 const double *clo = prob->clo_;
990 const double *cup = prob->cup_;
991
992 const double *rowels = prob->rowels_;
993 const int *hcol = prob->hcol_;
994 const CoinBigIndex *mrstrt = prob->mrstrt_;
995 int *hinrow = prob->hinrow_;
996 int nrows = prob->nrows_;
997
998 /*const*/ double *rlo = prob->rlo_;
999 /*const*/ double *rup = prob->rup_;
1000
1001 double *cost = prob->cost_;
1002
1003 presolvehlink *rlink = prob->rlink_;
1004 presolvehlink *clink = prob->clink_;
1005
1006 const unsigned char *integerType = prob->integerType_;
1007 bool stopSomeStuff = (prob->presolveOptions()&4)!=0;
1008
1009 const double tol = prob->feasibilityTolerance_;
1010#if 1
1011 // This needs to be made faster
1012 int numberInfeasible;
1013 //printf("Imp pass %d\n",prob->pass_);
1014#ifdef COIN_LIGHTWEIGHT_PRESOLVE
1015 if (prob->pass_==1) {
1016#endif
1017 next = testRedundant(prob,next,numberInfeasible);
1018 if ((prob->presolveOptions_&16384)!=0)
1019 numberInfeasible=0;
1020 if (numberInfeasible) {
1021 // infeasible
1022 prob->status_|= 1;
1023 return (next);
1024 }
1025#ifdef COIN_LIGHTWEIGHT_PRESOLVE
1026 }
1027#endif
1028 if (prob->pass_>15&&(prob->presolveOptions_&0x10000)!=0) {
1029 fill_level=2;
1030 return next;
1031 }
1032#endif
1033 // int nbounds = 0;
1034
1035 action *actions = new action [ncols];
1036# ifdef ZEROFAULT
1037 CoinZeroN(reinterpret_cast<char *>(actions),ncols*sizeof(action)) ;
1038# endif
1039 int nactions = 0;
1040 bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
1041
1042 int *implied_free = prob->usefulColumnInt_; //new int[ncols];
1043 int * whichFree = implied_free+ncols;
1044 int numberFree=0;
1045 int i;
1046
1047 // memory for min max
1048 int * infiniteDown = new int [nrows];
1049 int * infiniteUp = new int [nrows];
1050 double * maxDown = new double[nrows];
1051 double * maxUp = new double[nrows];
1052
1053 // mark as not computed
1054 // -1 => not computed, -2 give up (singleton), -3 give up (other)
1055 for (i=0;i<nrows;i++) {
1056 if (hinrow[i]>1)
1057 infiniteUp[i]=-1;
1058 else
1059 infiniteUp[i]=-2;
1060 }
1061#ifdef USE_SMALL_LARGE
1062 double large=1.0e10;
1063#else
1064 double large=1.0e20;
1065#endif
1066
1067 int numberLook = prob->numberColsToDo_;
1068 int iLook;
1069 int * look = prob->colsToDo_;
1070 //int * look2 = NULL;
1071 // if gone from 2 to 3 look at all
1072 // why does this loop not check for prohibited columns? -- lh, 040818 --
1073 // changed 040923
1074 if (fill_level<0) {
1075 look = prob->usefulColumnInt_+ncols; //new int[ncols];
1076 //look=look2;
1077 if (!prob->anyProhibited()) {
1078 for (iLook=0;iLook<ncols;iLook++)
1079 look[iLook]=iLook;
1080 numberLook=ncols;
1081 } else {
1082 // some prohibited
1083 numberLook=0;
1084 for (iLook=0;iLook<ncols;iLook++)
1085 if (!prob->colProhibited(iLook))
1086 look[numberLook++]=iLook;
1087 }
1088 }
1089 int maxLook = CoinAbs(fill_level);
1090 for (iLook=0;iLook<numberLook;iLook++) {
1091 int j=look[iLook];
1092 if (hincol[j] <= maxLook&&hincol[j]) {
1093 CoinBigIndex kcs = mcstrt[j];
1094 CoinBigIndex kce = kcs + hincol[j];
1095 bool singletonColumn = (hincol[j]==1);
1096 bool possible = false;
1097 bool singleton = false;
1098 CoinBigIndex k;
1099 double largestElement=0.0;
1100 for (k=kcs; k<kce; ++k) {
1101 int row = hrow[k];
1102 double coeffj = colels[k];
1103
1104 // if its row is an equality constraint...
1105 if (hinrow[row] > 1 ) {
1106 if ( fabs(rlo[row] - rup[row]) < tol &&
1107 fabs(coeffj) > ZTOLDP2) {
1108 possible=true;
1109 }
1110 largestElement = CoinMax(largestElement,fabs(coeffj));
1111 } else {
1112 singleton=true;
1113 }
1114 }
1115 if (possible&&!singleton) {
1116 double low=-COIN_DBL_MAX;
1117 double high=COIN_DBL_MAX;
1118 // get bound implied by all rows
1119 for (k=kcs; k<kce; ++k) {
1120 int row = hrow[k];
1121 double coeffj = colels[k];
1122 if (fabs(coeffj) > ZTOLDP2) {
1123 if (infiniteUp[row]==-1) {
1124 // compute
1125 CoinBigIndex krs = mrstrt[row];
1126 CoinBigIndex kre = krs + hinrow[row];
1127 int infiniteUpper = 0;
1128 int infiniteLower = 0;
1129 double maximumUp = 0.0;
1130 double maximumDown = 0.0;
1131 CoinBigIndex kk;
1132 // Compute possible lower and upper ranges
1133 for (kk = krs; kk < kre; ++kk) {
1134 double value=rowels[kk];
1135 int iColumn = hcol[kk];
1136 if (value > 0.0) {
1137 if (cup[iColumn] < large)
1138 maximumUp += cup[iColumn] * value;
1139 else
1140 ++infiniteUpper;
1141 if (clo[iColumn] > -large)
1142 maximumDown += clo[iColumn] * value;
1143 else
1144 ++infiniteLower;
1145 } else if (value<0.0) {
1146 if (cup[iColumn] < large)
1147 maximumDown += cup[iColumn] * value;
1148 else
1149 ++infiniteLower;
1150 if (clo[iColumn] > -large)
1151 maximumUp += clo[iColumn] * value;
1152 else
1153 ++infiniteUpper;
1154 }
1155 }
1156 double maxUpx = maximumUp+infiniteUpper*1.0e31;
1157 double maxDownx = maximumDown-infiniteLower*1.0e31;
1158 if (maxUpx <= rup[row] + tol &&
1159 maxDownx >= rlo[row] - tol) {
1160
1161 // Row is redundant
1162 infiniteUp[row]=-3;
1163
1164 } else if (maxUpx < rlo[row] -tol &&!fixInfeasibility) {
1165 /* there is an upper bound and it can't be reached */
1166 prob->status_|= 1;
1167 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
1168 prob->messages())
1169 <<row
1170 <<rlo[row]
1171 <<rup[row]
1172 <<CoinMessageEol;
1173 infiniteUp[row]=-3;
1174 break;
1175 } else if ( maxDownx > rup[row]+tol&&!fixInfeasibility) {
1176 /* there is a lower bound and it can't be reached */
1177 prob->status_|= 1;
1178 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
1179 prob->messages())
1180 <<row
1181 <<rlo[row]
1182 <<rup[row]
1183 <<CoinMessageEol;
1184 infiniteUp[row]=-3;
1185 break;
1186 } else {
1187 infiniteUp[row]=infiniteUpper;
1188 infiniteDown[row]=infiniteLower;
1189 maxUp[row]=maximumUp;
1190 maxDown[row]=maximumDown;
1191 }
1192 }
1193 if (infiniteUp[row]>=0) {
1194 double lower = rlo[row];
1195 double upper = rup[row];
1196 double value=coeffj;
1197 double nowLower = clo[j];
1198 double nowUpper = cup[j];
1199 double newBound;
1200 int infiniteUpper=infiniteUp[row];
1201 int infiniteLower=infiniteDown[row];
1202 double maximumUp = maxUp[row];
1203 double maximumDown = maxDown[row];
1204 if (value > 0.0) {
1205 // positive value
1206 if (lower>-large) {
1207 if (!infiniteUpper) {
1208 assert(nowUpper < large);
1209 newBound = nowUpper +
1210 (lower - maximumUp) / value;
1211 // relax if original was large
1212 if (fabs(maximumUp)>1.0e8&&!singletonColumn)
1213 newBound -= 1.0e-12*fabs(maximumUp);
1214 } else if (infiniteUpper==1&&nowUpper>large) {
1215 newBound = (lower -maximumUp) / value;
1216 // relax if original was large
1217 if (fabs(maximumUp)>1.0e8&&!singletonColumn)
1218 newBound -= 1.0e-12*fabs(maximumUp);
1219 } else {
1220 newBound = -COIN_DBL_MAX;
1221 }
1222 if (newBound<=-large)
1223 newBound = -COIN_DBL_MAX;
1224 if (newBound > nowLower + 1.0e-12) {
1225 // Tighten the lower bound
1226 // adjust
1227 double now;
1228 if (nowLower<-large) {
1229 now=0.0;
1230 infiniteLower--;
1231 } else {
1232 now = nowLower;
1233 }
1234 maximumDown += (newBound-now) * value;
1235 nowLower = newBound;
1236 }
1237 low=CoinMax(low,newBound);
1238 }
1239 if (upper <large) {
1240 if (!infiniteLower) {
1241 assert(nowLower >- large);
1242 newBound = nowLower +
1243 (upper - maximumDown) / value;
1244 // relax if original was large
1245 if (fabs(maximumDown)>1.0e8&&!singletonColumn)
1246 newBound += 1.0e-12*fabs(maximumDown);
1247 } else if (infiniteLower==1&&nowLower<-large) {
1248 newBound = (upper - maximumDown) / value;
1249 // relax if original was large
1250 if (fabs(maximumDown)>1.0e8&&!singletonColumn)
1251 newBound += 1.0e-12*fabs(maximumDown);
1252 } else {
1253 newBound = COIN_DBL_MAX;
1254 }
1255 if (newBound>=large)
1256 newBound = COIN_DBL_MAX;
1257 if (newBound < nowUpper - 1.0e-12) {
1258 // Tighten the upper bound
1259 // adjust
1260 double now;
1261 if (nowUpper>large) {
1262 now=0.0;
1263 infiniteUpper--;
1264 } else {
1265 now = nowUpper;
1266 }
1267 maximumUp += (newBound-now) * value;
1268 nowUpper = newBound;
1269 }
1270 high=CoinMin(high,newBound);
1271 }
1272 } else {
1273 // negative value
1274 if (lower>-large) {
1275 if (!infiniteUpper) {
1276 assert(nowLower >- large);
1277 newBound = nowLower +
1278 (lower - maximumUp) / value;
1279 // relax if original was large
1280 if (fabs(maximumUp)>1.0e8&&!singletonColumn)
1281 newBound += 1.0e-12*fabs(maximumUp);
1282 } else if (infiniteUpper==1&&nowLower<-large) {
1283 newBound = (lower -maximumUp) / value;
1284 // relax if original was large
1285 if (fabs(maximumUp)>1.0e8&&!singletonColumn)
1286 newBound += 1.0e-12*fabs(maximumUp);
1287 } else {
1288 newBound = COIN_DBL_MAX;
1289 }
1290 if (newBound>=large)
1291 newBound = COIN_DBL_MAX;
1292 if (newBound < nowUpper - 1.0e-12) {
1293 // Tighten the upper bound
1294 // adjust
1295 double now;
1296 if (nowUpper>large) {
1297 now=0.0;
1298 infiniteLower--;
1299 } else {
1300 now = nowUpper;
1301 }
1302 maximumDown += (newBound-now) * value;
1303 nowUpper = newBound;
1304 }
1305 high=CoinMin(high,newBound);
1306 }
1307 if (upper <large) {
1308 if (!infiniteLower) {
1309 assert(nowUpper < large);
1310 newBound = nowUpper +
1311 (upper - maximumDown) / value;
1312 // relax if original was large
1313 if (fabs(maximumDown)>1.0e8&&!singletonColumn)
1314 newBound -= 1.0e-12*fabs(maximumDown);
1315 } else if (infiniteLower==1&&nowUpper>large) {
1316 newBound = (upper - maximumDown) / value;
1317 // relax if original was large
1318 if (fabs(maximumDown)>1.0e8&&!singletonColumn)
1319 newBound -= 1.0e-12*fabs(maximumDown);
1320 } else {
1321 newBound = -COIN_DBL_MAX;
1322 }
1323 if (newBound<=-large)
1324 newBound = -COIN_DBL_MAX;
1325 if (newBound > nowLower + 1.0e-12) {
1326 // Tighten the lower bound
1327 // adjust
1328 double now;
1329 if (nowLower<-large) {
1330 now=0.0;
1331 infiniteUpper--;
1332 } else {
1333 now = nowLower;
1334 }
1335 maximumUp += (newBound-now) * value;
1336 nowLower = newBound;
1337 }
1338 low = CoinMax(low,newBound);
1339 }
1340 }
1341 } else if (infiniteUp[row]==-3) {
1342 // give up
1343 high=COIN_DBL_MAX;
1344 low=-COIN_DBL_MAX;
1345 break;
1346 }
1347 }
1348 }
1349 if (clo[j] <= low && high <= cup[j]) {
1350
1351 // both column bounds implied by the constraints of the problem
1352 // get row
1353 // If more than one equality is present, how do I know the one I
1354 // select here will be the one that actually implied tighter
1355 // bounds? Seems like I should care. -- lh, 040818 --
1356 largestElement *= 0.1;
1357 int krow=-1;
1358 int ninrow=ncols+1;
1359 double thisValue=0.0;
1360 for (k=kcs; k<kce; ++k) {
1361 int row = hrow[k];
1362 double coeffj = colels[k];
1363 if ( fabs(rlo[row] - rup[row]) < tol &&
1364 fabs(coeffj) > largestElement) {
1365 if (hinrow[row]<ninrow) {
1366 ninrow=hinrow[row];
1367 krow=row;
1368 thisValue=coeffj;
1369 }
1370 }
1371 }
1372 if (krow>=0) {
1373 bool goodRow=true;
1374 if (integerType[j]) {
1375 // can only accept if good looking row
1376 double scaleFactor = 1.0/thisValue;
1377 double rhs = rlo[krow]*scaleFactor;
1378 if (fabs(rhs-floor(rhs+0.5))<tol) {
1379 CoinBigIndex krs = mrstrt[krow];
1380 CoinBigIndex kre = krs + hinrow[krow];
1381 CoinBigIndex kk;
1382 bool allOnes=true;
1383 for (kk = krs; kk < kre; ++kk) {
1384 double value=rowels[kk]*scaleFactor;
1385 if (fabs(value)!=1.0)
1386 allOnes=false;
1387 int iColumn = hcol[kk];
1388 if (!integerType[iColumn]||fabs(value-floor(value+0.5))>tol) {
1389 goodRow=false;
1390 break;
1391 }
1392 }
1393 if (rlo[krow]==1.0&&hinrow[krow]>=5&&stopSomeStuff&&allOnes)
1394 goodRow=false; // may spoil SOS
1395 } else {
1396 goodRow=false;
1397 }
1398 }
1399 if (goodRow) {
1400 implied_free[numberFree] = krow;
1401 whichFree[numberFree++] = j;
1402 // And say row no good for further use
1403 infiniteUp[krow]=-3;
1404 //printf("column %d implied free by row %d hincol %d hinrow %d\n",
1405 // j,krow,hincol[j],hinrow[krow]);
1406 }
1407 }
1408 }
1409 }
1410 }
1411 }
1412
1413 delete [] infiniteDown;
1414 delete [] infiniteUp;
1415 delete [] maxDown;
1416 delete [] maxUp;
1417
1418 int isolated_row = -1;
1419
1420 // first pick off the easy ones
1421 // note that this will only deal with columns that were originally
1422 // singleton; it will not deal with doubleton columns that become
1423 // singletons as a result of dropping rows.
1424 for (iLook=0;iLook<numberFree;iLook++) {
1425 int j=whichFree[iLook];
1426 if (hincol[j] == 1) {
1427 CoinBigIndex kcs = mcstrt[j];
1428 int row = hrow[kcs];
1429 double coeffj = colels[kcs];
1430
1431 CoinBigIndex krs = mrstrt[row];
1432 CoinBigIndex kre = krs + hinrow[row];
1433
1434
1435 // isolated rows are weird
1436 {
1437 int n = 0;
1438 for (CoinBigIndex k=krs; k<kre; ++k)
1439 n += hincol[hcol[k]];
1440 if (n==hinrow[row]) {
1441 isolated_row = row;
1442 break;
1443 }
1444 }
1445
1446 const bool nonzero_cost = (cost[j] != 0.0&&fabs(rup[row]-rlo[row])<=tol);
1447
1448 double *save_costs = nonzero_cost ? new double[hinrow[row]] : NULL;
1449
1450 {
1451 action *s = &actions[nactions++];
1452
1453 s->row = row;
1454 s->col = j;
1455 PRESOLVE_DETAIL_PRINT(printf("pre_impliedfree2 %dC %dR E\n",j,row));
1456
1457 s->clo = clo[j];
1458 s->cup = cup[j];
1459 s->rlo = rlo[row];
1460 s->rup = rup[row];
1461
1462 s->ninrow = hinrow[row];
1463 s->rowels = presolve_dupmajor(rowels,hcol,hinrow[row],krs) ;
1464 s->costs = save_costs;
1465 }
1466
1467 if (nonzero_cost) {
1468 double rhs = rlo[row];
1469 double costj = cost[j];
1470
1471# if PRESOLVE_DEBUG
1472 printf("FREE COSTS: %g ", costj);
1473# endif
1474 for (CoinBigIndex k=krs; k<kre; k++) {
1475 int jcol = hcol[k];
1476 save_costs[k-krs] = cost[jcol];
1477
1478 if (jcol != j) {
1479 double coeff = rowels[k];
1480
1481# if PRESOLVE_DEBUG
1482 printf("%g %g ", cost[jcol], coeff/coeffj);
1483# endif
1484 /*
1485 * Similar to eliminating doubleton:
1486 * cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
1487 * cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
1488 */
1489 cost[jcol] += costj * (-coeff / coeffj);
1490 }
1491 }
1492# if PRESOLVE_DEBUG
1493 printf("\n");
1494
1495 /* similar to doubleton */
1496 printf("BIAS??????? %g %g %g %g\n",
1497 costj * rhs / coeffj,
1498 costj, rhs, coeffj);
1499# endif
1500 prob->change_bias(costj * rhs / coeffj);
1501 // ??
1502 cost[j] = 0.0;
1503 }
1504
1505 /* remove the row from the columns in the row */
1506 for (CoinBigIndex k=krs; k<kre; k++) {
1507 int jcol=hcol[k];
1508 prob->addCol(jcol);
1509 presolve_delete_from_col(row,jcol,mcstrt,hincol,hrow,colels) ;
1510 if (hincol[jcol] == 0)
1511 { PRESOLVE_REMOVE_LINK(prob->clink_,jcol) ; }
1512 }
1513 PRESOLVE_REMOVE_LINK(rlink, row);
1514 hinrow[row] = 0;
1515
1516 // just to make things squeeky
1517 rlo[row] = 0.0;
1518 rup[row] = 0.0;
1519
1520 PRESOLVE_REMOVE_LINK(clink, j);
1521 hincol[j] = 0;
1522
1523 implied_free[iLook] = -1;
1524 }
1525 }
1526
1527 //delete [] look2;
1528 if (nactions) {
1529# if PRESOLVE_SUMMARY
1530 printf("NIMPLIED FREE: %d\n", nactions);
1531# endif
1532 action *actions1 = new action[nactions];
1533 CoinMemcpyN(actions, nactions, actions1);
1534 next = new implied_free_action(nactions, actions1, next);
1535 }
1536 delete [] actions;
1537
1538 if (isolated_row != -1) {
1539 const CoinPresolveAction *nextX = isolated_constraint_action::presolve(prob,
1540 isolated_row, next);
1541 if (nextX)
1542 next = nextX; // may fail
1543 }
1544 // try more complex ones
1545 if (fill_level) {
1546 next = subst_constraint_action::presolve(prob, implied_free,
1547 whichFree,numberFree,
1548 next,fill_level);
1549 }
1550 //delete[]implied_free;
1551
1552 if (prob->tuning_) {
1553 double thisTime=CoinCpuTime();
1554 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
1555 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
1556 printf("CoinPresolveImpliedFree(64) - %d rows, %d columns dropped in time %g, total %g\n",
1557 droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
1558 }
1559 return (next);
1560}
1561
1562
1563
1564const char *implied_free_action::name() const
1565{
1566 return ("implied_free_action");
1567}
1568
1569void implied_free_action::postsolve(CoinPostsolveMatrix *prob) const
1570{
1571 const action *const actions = actions_;
1572 const int nactions = nactions_;
1573
1574 double *elementByColumn = prob->colels_;
1575 int *hrow = prob->hrow_;
1576 CoinBigIndex *columnStart = prob->mcstrt_;
1577 int *numberInColumn = prob->hincol_;
1578 int *link = prob->link_;
1579
1580 double *clo = prob->clo_;
1581 double *cup = prob->cup_;
1582
1583 double *rlo = prob->rlo_;
1584 double *rup = prob->rup_;
1585
1586 double *sol = prob->sol_;
1587
1588 double *rcosts = prob->rcosts_;
1589 double *dcost = prob->cost_;
1590
1591 double *acts = prob->acts_;
1592 double *rowduals = prob->rowduals_;
1593
1594 const double maxmin = prob->maxmin_;
1595
1596# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1597 char *cdone = prob->cdone_;
1598 char *rdone = prob->rdone_;
1599# endif
1600
1601 CoinBigIndex &free_list = prob->free_list_;
1602
1603 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
1604
1605 int irow = f->row;
1606 int icol = f->col;
1607
1608 int ninrow = f->ninrow;
1609 const double *rowels = f->rowels;
1610 const int *rowcols = reinterpret_cast<const int *>(rowels+ninrow) ;
1611 const double *save_costs = f->costs;
1612
1613 // put back coefficients in the row
1614 // this includes recreating the singleton column
1615 {
1616 for (int k = 0; k<ninrow; k++) {
1617 int jcol = rowcols[k];
1618 double coeff = rowels[k];
1619
1620 if (save_costs) {
1621 rcosts[jcol] += maxmin*(save_costs[k]-dcost[jcol]);
1622 dcost[jcol] = save_costs[k];
1623 }
1624 {
1625 CoinBigIndex kk = free_list;
1626 assert(kk >= 0 && kk < prob->bulk0_) ;
1627 free_list = link[free_list];
1628 link[kk] = columnStart[jcol];
1629 columnStart[jcol] = kk;
1630 elementByColumn[kk] = coeff;
1631 hrow[kk] = irow;
1632 }
1633
1634 if (jcol == icol) {
1635 // initialize the singleton column
1636 numberInColumn[jcol] = 1;
1637 clo[icol] = f->clo;
1638 cup[icol] = f->cup;
1639
1640# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1641 cdone[icol] = IMPLIED_FREE;
1642# endif
1643 } else {
1644 numberInColumn[jcol]++;
1645 }
1646 }
1647# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1648 rdone[irow] = IMPLIED_FREE;
1649# endif
1650# if PRESOLVE_CONSISTENCY
1651 presolve_check_free_list(prob) ;
1652# endif
1653
1654 rlo[irow] = f->rlo;
1655 rup[irow] = f->rup;
1656 }
1657 //deleteAction( f->costs,double*); do on delete
1658 // coeff has now been initialized
1659
1660 // compute solution
1661 {
1662 double act = 0.0;
1663 double coeff = 0.0;
1664
1665 for (int k = 0; k<ninrow; k++)
1666 if (rowcols[k] == icol)
1667 coeff = rowels[k];
1668 else {
1669 int jcol = rowcols[k];
1670 PRESOLVE_STMT(CoinBigIndex kk = presolve_find_row2(irow, columnStart[jcol], numberInColumn[jcol], hrow, link));
1671 act += rowels[k] * sol[jcol];
1672 }
1673
1674 PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
1675 double thisCost = maxmin*dcost[icol];
1676 double loActivity,upActivity;
1677 if (coeff>0) {
1678 loActivity = (rlo[irow]-act)/coeff;
1679 upActivity = (rup[irow]-act)/coeff;
1680 } else {
1681 loActivity = (rup[irow]-act)/coeff;
1682 upActivity = (rlo[irow]-act)/coeff;
1683 }
1684 loActivity = CoinMax(loActivity,clo[icol]);
1685 upActivity = CoinMin(upActivity,cup[icol]);
1686 int where; //0 in basis, -1 at lb, +1 at ub
1687 const double tolCheck = 0.1*prob->ztolzb_;
1688 if (loActivity<clo[icol]+tolCheck/fabs(coeff)&&thisCost>=0.0)
1689 where=-1;
1690 else if (upActivity>cup[icol]-tolCheck/fabs(coeff)&&thisCost<0.0)
1691 where=1;
1692 else
1693 where =0;
1694 // But we may need to put in basis to stay dual feasible
1695 double possibleDual = thisCost/coeff;
1696 if (where) {
1697 double worst= prob->ztoldj_;
1698 for (int k = 0; k<ninrow; k++) {
1699 int jcol = rowcols[k];
1700 if (jcol!=icol) {
1701 CoinPrePostsolveMatrix::Status status = prob->getColumnStatus(jcol);
1702 // can only trust basic
1703 if (status==CoinPrePostsolveMatrix::basic) {
1704 if (fabs(rcosts[jcol])>worst)
1705 worst=fabs(rcosts[jcol]);
1706 } else if (sol[jcol]<clo[jcol]+ZTOLDP) {
1707 if (-rcosts[jcol]>worst)
1708 worst=-rcosts[jcol];
1709 } else if (sol[jcol]>cup[jcol]-ZTOLDP) {
1710 if (rcosts[jcol]>worst)
1711 worst=rcosts[jcol];
1712 }
1713 }
1714 }
1715 if (worst>prob->ztoldj_) {
1716 // see if better if in basis
1717 double worst2 = prob->ztoldj_;
1718 for (int k = 0; k<ninrow; k++) {
1719 int jcol = rowcols[k];
1720 if (jcol!=icol) {
1721 double coeff = rowels[k];
1722 double newDj = rcosts[jcol]-possibleDual*coeff;
1723 CoinPrePostsolveMatrix::Status status = prob->getColumnStatus(jcol);
1724 // can only trust basic
1725 if (status==CoinPrePostsolveMatrix::basic) {
1726 if (fabs(newDj)>worst2)
1727 worst2=fabs(newDj);
1728 } else if (sol[jcol]<clo[jcol]+ZTOLDP) {
1729 if (-newDj>worst2)
1730 worst2=-newDj;
1731 } else if (sol[jcol]>cup[jcol]-ZTOLDP) {
1732 if (newDj>worst2)
1733 worst2=newDj;
1734 }
1735 }
1736 }
1737 if (worst2<worst)
1738 where=0; // put in basis
1739 }
1740 }
1741 if (!where) {
1742 // choose rowdual to make this col basic
1743 rowduals[irow] = possibleDual;
1744 if ((rlo[irow] < rup[irow] && rowduals[irow] < 0.0)
1745 || rlo[irow]< -1.0e20) {
1746 //if (rlo[irow]<-1.0e20&&rowduals[irow]>ZTOLDP)
1747 //printf("IMP %g %g %g\n",rlo[irow],rup[irow],rowduals[irow]);
1748 sol[icol] = (rup[irow] - act) / coeff;
1749 //assert (sol[icol]>=clo[icol]-1.0e-5&&sol[icol]<=cup[icol]+1.0e-5);
1750 acts[irow] = rup[irow];
1751 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atUpperBound);
1752 } else {
1753 sol[icol] = (rlo[irow] - act) / coeff;
1754 //assert (sol[icol]>=clo[icol]-1.0e-5&&sol[icol]<=cup[icol]+1.0e-5);
1755 acts[irow] = rlo[irow];
1756 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
1757 }
1758 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::basic);
1759 for (int k = 0; k<ninrow; k++) {
1760 int jcol = rowcols[k];
1761 double coeff = rowels[k];
1762 rcosts[jcol] -= possibleDual*coeff;
1763 }
1764 rcosts[icol] = 0.0;
1765 } else {
1766 rowduals[irow] = 0.0;
1767 rcosts[icol] = thisCost;
1768 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
1769 if (where<0) {
1770 // to lb
1771 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
1772 sol[icol]=clo[icol];
1773 } else {
1774 // to ub
1775 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
1776 sol[icol]=cup[icol];
1777 }
1778 acts[irow] = act + sol[icol]*coeff;
1779 //assert (acts[irow]>=rlo[irow]-1.0e-5&&acts[irow]<=rup[irow]+1.0e-5);
1780 }
1781# if PRESOLVE_DEBUG
1782 {
1783 double *colels = prob->colels_;
1784 int *hrow = prob->hrow_;
1785 const CoinBigIndex *mcstrt = prob->mcstrt_;
1786 int *hincol = prob->hincol_;
1787 for (int j = 0; j<ninrow; j++) {
1788 int jcol = rowcols[j];
1789 CoinBigIndex k = mcstrt[jcol];
1790 int nx = hincol[jcol];
1791 double dj = dcost[jcol];
1792 for (int i=0; i<nx; ++i) {
1793 int row = hrow[k];
1794 double coeff = colels[k];
1795 k = link[k];
1796 dj -= rowduals[row] * coeff;
1797 //printf("col jcol row %d coeff %g dual %g new dj %g\n",
1798 // row,coeff,rowduals[row],dj);
1799 }
1800 if (fabs(dj-rcosts[jcol])>1.0e-3)
1801 printf("changed\n");
1802 }
1803 }
1804# endif
1805 }
1806 }
1807
1808 return ;
1809}
1810
1811
1812/*
1813 Why do we delete costs during postsolve() execution, but none of the other
1814 components of the action?
1815*/
1816implied_free_action::~implied_free_action()
1817{
1818 int i;
1819 for (i=0;i<nactions_;i++) {
1820 deleteAction(actions_[i].rowels,double *);
1821 deleteAction( actions_[i].costs,double *);
1822 }
1823 deleteAction(actions_,action *);
1824}
1825