1 /* $Id: ClpPresolve.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
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 //#define PRESOLVE_CONSISTENCY 1
7 //#define PRESOLVE_DEBUG 1
8
9 #include <stdio.h>
10
11 #include <cassert>
12 #include <iostream>
13
14 #include "CoinHelperFunctions.hpp"
15 #include "ClpConfig.h"
16 #ifdef CLP_HAS_ABC
17 #include "CoinAbcCommon.hpp"
18 #endif
19
20 #include "CoinPackedMatrix.hpp"
21 #include "ClpPackedMatrix.hpp"
22 #include "ClpSimplex.hpp"
23 #include "ClpSimplexOther.hpp"
24 #ifndef SLIM_CLP
25 #include "ClpQuadraticObjective.hpp"
26 #endif
27
28 #include "ClpPresolve.hpp"
29 #include "CoinPresolveMatrix.hpp"
30
31 #include "CoinPresolveEmpty.hpp"
32 #include "CoinPresolveFixed.hpp"
33 #include "CoinPresolvePsdebug.hpp"
34 #include "CoinPresolveSingleton.hpp"
35 #include "CoinPresolveDoubleton.hpp"
36 #include "CoinPresolveTripleton.hpp"
37 #include "CoinPresolveZeros.hpp"
38 #include "CoinPresolveSubst.hpp"
39 #include "CoinPresolveForcing.hpp"
40 #include "CoinPresolveDual.hpp"
41 #include "CoinPresolveTighten.hpp"
42 #include "CoinPresolveUseless.hpp"
43 #include "CoinPresolveDupcol.hpp"
44 #include "CoinPresolveImpliedFree.hpp"
45 #include "CoinPresolveIsolated.hpp"
46 #include "CoinMessage.hpp"
47
ClpPresolve()48 ClpPresolve::ClpPresolve()
49 : originalModel_(NULL)
50 , presolvedModel_(NULL)
51 , nonLinearValue_(0.0)
52 , originalColumn_(NULL)
53 , originalRow_(NULL)
54 , rowObjective_(NULL)
55 , paction_(0)
56 , ncols_(0)
57 , nrows_(0)
58 , nelems_(0)
59 ,
60 #ifdef CLP_INHERIT_MODE
61 numberPasses_(20)
62 ,
63 #else
64 numberPasses_(5)
65 ,
66 #endif
67 substitution_(3)
68 ,
69 #ifndef CLP_NO_STD
70 saveFile_("")
71 ,
72 #endif
73 presolveActions_(0)
74 {
75 }
76
~ClpPresolve()77 ClpPresolve::~ClpPresolve()
78 {
79 destroyPresolve();
80 }
81 // Gets rid of presolve actions (e.g.when infeasible)
destroyPresolve()82 void ClpPresolve::destroyPresolve()
83 {
84 const CoinPresolveAction *paction = paction_;
85 while (paction) {
86 const CoinPresolveAction *next = paction->next;
87 delete paction;
88 paction = next;
89 }
90 delete[] originalColumn_;
91 delete[] originalRow_;
92 paction_ = NULL;
93 originalColumn_ = NULL;
94 originalRow_ = NULL;
95 delete[] rowObjective_;
96 rowObjective_ = NULL;
97 }
98
99 /* This version of presolve returns a pointer to a new presolved
100 model. NULL if infeasible
101 */
102 ClpSimplex *
presolvedModel(ClpSimplex & si,double feasibilityTolerance,bool keepIntegers,int numberPasses,bool dropNames,bool doRowObjective,const char * prohibitedRows,const char * prohibitedColumns)103 ClpPresolve::presolvedModel(ClpSimplex &si,
104 double feasibilityTolerance,
105 bool keepIntegers,
106 int numberPasses,
107 bool dropNames,
108 bool doRowObjective,
109 const char *prohibitedRows,
110 const char *prohibitedColumns)
111 {
112 // Check matrix
113 int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
114 if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
115 1.0e20, checkType))
116 return NULL;
117 else
118 return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
119 doRowObjective,
120 prohibitedRows,
121 prohibitedColumns);
122 }
123 #ifndef CLP_NO_STD
124 /* This version of presolve updates
125 model and saves original data to file. Returns non-zero if infeasible
126 */
presolvedModelToFile(ClpSimplex & si,std::string fileName,double feasibilityTolerance,bool keepIntegers,int numberPasses,bool dropNames,bool doRowObjective)127 int ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
128 double feasibilityTolerance,
129 bool keepIntegers,
130 int numberPasses,
131 bool dropNames,
132 bool doRowObjective)
133 {
134 // Check matrix
135 if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
136 1.0e20))
137 return 2;
138 saveFile_ = fileName;
139 si.saveModel(saveFile_.c_str());
140 ClpSimplex *model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
141 doRowObjective);
142 if (model == &si) {
143 return 0;
144 } else {
145 si.restoreModel(saveFile_.c_str());
146 remove(saveFile_.c_str());
147 return 1;
148 }
149 }
150 #endif
151 // Return pointer to presolved model
152 ClpSimplex *
model() const153 ClpPresolve::model() const
154 {
155 return presolvedModel_;
156 }
157 // Return pointer to original model
158 ClpSimplex *
originalModel() const159 ClpPresolve::originalModel() const
160 {
161 return originalModel_;
162 }
163 // Return presolve status (0,1,2)
presolveStatus() const164 int ClpPresolve::presolveStatus() const
165 {
166 if (nelems_ >= 0) {
167 // feasible (or not done yet)
168 return 0;
169 } else {
170 int presolveStatus = -static_cast< int >(nelems_);
171 // If both infeasible and unbounded - say infeasible
172 if (presolveStatus > 2)
173 presolveStatus = 1;
174 return presolveStatus;
175 }
176 }
postsolve(bool updateStatus)177 void ClpPresolve::postsolve(bool updateStatus)
178 {
179 // Return at once if no presolved model
180 if (!presolvedModel_)
181 return;
182 // Messages
183 CoinMessages messages = originalModel_->coinMessages();
184 if (!presolvedModel_->isProvenOptimal()) {
185 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
186 messages)
187 << CoinMessageEol;
188 }
189
190 // this is the size of the original problem
191 const int ncols0 = ncols_;
192 const int nrows0 = nrows_;
193 const CoinBigIndex nelems0 = nelems_;
194
195 // this is the reduced problem
196 int ncols = presolvedModel_->getNumCols();
197 int nrows = presolvedModel_->getNumRows();
198
199 double *acts = NULL;
200 double *sol = NULL;
201 unsigned char *rowstat = NULL;
202 unsigned char *colstat = NULL;
203 #ifndef CLP_NO_STD
204 if (saveFile_ == "") {
205 #endif
206 // reality check
207 assert(ncols0 == originalModel_->getNumCols());
208 assert(nrows0 == originalModel_->getNumRows());
209 acts = originalModel_->primalRowSolution();
210 sol = originalModel_->primalColumnSolution();
211 if (updateStatus) {
212 // postsolve does not know about fixed
213 int i;
214 for (i = 0; i < nrows + ncols; i++) {
215 if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
216 presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
217 }
218 unsigned char *status = originalModel_->statusArray();
219 if (!status) {
220 originalModel_->createStatus();
221 status = originalModel_->statusArray();
222 }
223 rowstat = status + ncols0;
224 colstat = status;
225 CoinMemcpyN(presolvedModel_->statusArray(), ncols, colstat);
226 CoinMemcpyN(presolvedModel_->statusArray() + ncols, nrows, rowstat);
227 }
228 #ifndef CLP_NO_STD
229 } else {
230 // from file
231 acts = new double[nrows0];
232 sol = new double[ncols0];
233 CoinZeroN(acts, nrows0);
234 CoinZeroN(sol, ncols0);
235 if (updateStatus) {
236 unsigned char *status = new unsigned char[nrows0 + ncols0];
237 rowstat = status + ncols0;
238 colstat = status;
239 CoinMemcpyN(presolvedModel_->statusArray(), ncols, colstat);
240 CoinMemcpyN(presolvedModel_->statusArray() + ncols, nrows, rowstat);
241 }
242 }
243 #endif
244
245 // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
246 // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
247 // cause duplicate free. In case where saveFile == "", as best I can see
248 // arrays are owned by originalModel_. fix is to
249 // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
250 CoinPostsolveMatrix prob(presolvedModel_,
251 ncols0,
252 nrows0,
253 nelems0,
254 presolvedModel_->getObjSense(),
255 // end prepost
256
257 sol, acts,
258 colstat, rowstat);
259
260 postsolve(prob);
261
262 #ifndef CLP_NO_STD
263 if (saveFile_ != "") {
264 // From file
265 assert(originalModel_ == presolvedModel_);
266 originalModel_->restoreModel(saveFile_.c_str());
267 remove(saveFile_.c_str());
268 CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
269 // delete [] acts;
270 CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
271 // delete [] sol;
272 if (updateStatus) {
273 CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
274 // delete [] colstat;
275 }
276 } else {
277 #endif
278 prob.sol_ = 0;
279 prob.acts_ = 0;
280 prob.colstat_ = 0;
281 #ifndef CLP_NO_STD
282 }
283 #endif
284 // put back duals
285 CoinMemcpyN(prob.rowduals_, nrows_, originalModel_->dualRowSolution());
286 double maxmin = originalModel_->getObjSense();
287 if (maxmin < 0.0) {
288 // swap signs
289 int i;
290 double *pi = originalModel_->dualRowSolution();
291 for (i = 0; i < nrows_; i++)
292 pi[i] = -pi[i];
293 }
294 // Now check solution
295 double offset;
296 CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
297 originalModel_->primalColumnSolution(), offset, true),
298 ncols_, originalModel_->dualColumnSolution());
299 originalModel_->clpMatrix()->transposeTimes(-1.0,
300 originalModel_->dualRowSolution(),
301 originalModel_->dualColumnSolution());
302 memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
303 originalModel_->clpMatrix()->times(1.0,
304 originalModel_->primalColumnSolution(),
305 originalModel_->primalRowSolution());
306 originalModel_->checkSolutionInternal();
307 if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
308 // See if we can fix easily
309 static_cast< ClpSimplexOther * >(originalModel_)->cleanupAfterPostsolve();
310 }
311 // Messages
312 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
313 messages)
314 << originalModel_->objectiveValue()
315 << originalModel_->sumDualInfeasibilities()
316 << originalModel_->numberDualInfeasibilities()
317 << originalModel_->sumPrimalInfeasibilities()
318 << originalModel_->numberPrimalInfeasibilities()
319 << CoinMessageEol;
320
321 //originalModel_->objectiveValue_=objectiveValue_;
322 originalModel_->setNumberIterations(presolvedModel_->numberIterations());
323 if (!presolvedModel_->status()) {
324 if (!originalModel_->numberDualInfeasibilities() && !originalModel_->numberPrimalInfeasibilities()) {
325 originalModel_->setProblemStatus(0);
326 } else {
327 originalModel_->setProblemStatus(-1);
328 // Say not optimal after presolve
329 originalModel_->setSecondaryStatus(7);
330 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
331 messages)
332 << CoinMessageEol;
333 }
334 } else {
335 originalModel_->setProblemStatus(presolvedModel_->status());
336 // but not if close to feasible
337 if (originalModel_->sumPrimalInfeasibilities() < 1.0e-1) {
338 originalModel_->setProblemStatus(-1);
339 // Say not optimal after presolve
340 originalModel_->setSecondaryStatus(7);
341 }
342 }
343 #ifndef CLP_NO_STD
344 if (saveFile_ != "")
345 presolvedModel_ = NULL;
346 #endif
347 }
348
349 // return pointer to original columns
350 const int *
originalColumns() const351 ClpPresolve::originalColumns() const
352 {
353 return originalColumn_;
354 }
355 // return pointer to original rows
356 const int *
originalRows() const357 ClpPresolve::originalRows() const
358 {
359 return originalRow_;
360 }
361 // Set pointer to original model
setOriginalModel(ClpSimplex * model)362 void ClpPresolve::setOriginalModel(ClpSimplex *model)
363 {
364 originalModel_ = model;
365 }
366 #if 0
367 // A lazy way to restrict which transformations are applied
368 // during debugging.
369 static int ATOI(const char *name)
370 {
371 return true;
372 #if PRESOLVE_DEBUG || PRESOLVE_SUMMARY
373 if (getenv(name)) {
374 int val = atoi(getenv(name));
375 printf("%s = %d\n", name, val);
376 return (val);
377 } else {
378 if (strcmp(name, "off"))
379 return (true);
380 else
381 return (false);
382 }
383 #else
384 return (true);
385 #endif
386 }
387 #endif
388 #ifdef PRESOLVE_DEBUG
389 #define PRESOLVE_CHECK_SOL 1
390 #endif
391 //#define PRESOLVE_CHECK_SOL 1
392 #if PRESOLVE_CHECK_SOL
check_sol(CoinPresolveMatrix * prob,double tol)393 void check_sol(CoinPresolveMatrix *prob, double tol)
394 {
395 double *colels = prob->colels_;
396 int *hrow = prob->hrow_;
397 int *mcstrt = prob->mcstrt_;
398 int *hincol = prob->hincol_;
399 int *hinrow = prob->hinrow_;
400 int ncols = prob->ncols_;
401
402 double *csol = prob->sol_;
403 double *acts = prob->acts_;
404 double *clo = prob->clo_;
405 double *cup = prob->cup_;
406 int nrows = prob->nrows_;
407 double *rlo = prob->rlo_;
408 double *rup = prob->rup_;
409
410 int colx;
411
412 double *rsol = new double[nrows];
413 memset(rsol, 0, nrows * sizeof(double));
414
415 for (colx = 0; colx < ncols; ++colx) {
416 if (1) {
417 CoinBigIndex k = mcstrt[colx];
418 int nx = hincol[colx];
419 double solutionValue = csol[colx];
420 for (int i = 0; i < nx; ++i) {
421 int row = hrow[k];
422 double coeff = colels[k];
423 k++;
424 rsol[row] += solutionValue * coeff;
425 }
426 if (csol[colx] < clo[colx] - tol) {
427 printf("low CSOL: %d - %g %g %g\n",
428 colx, clo[colx], csol[colx], cup[colx]);
429 } else if (csol[colx] > cup[colx] + tol) {
430 printf("high CSOL: %d - %g %g %g\n",
431 colx, clo[colx], csol[colx], cup[colx]);
432 }
433 }
434 }
435 int rowx;
436 for (rowx = 0; rowx < nrows; ++rowx) {
437 if (hinrow[rowx]) {
438 if (fabs(rsol[rowx] - acts[rowx]) > tol)
439 printf("inacc RSOL: %d - %g %g (acts_ %g) %g\n",
440 rowx, rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
441 if (rsol[rowx] < rlo[rowx] - tol) {
442 printf("low RSOL: %d - %g %g %g\n",
443 rowx, rlo[rowx], rsol[rowx], rup[rowx]);
444 } else if (rsol[rowx] > rup[rowx] + tol) {
445 printf("high RSOL: %d - %g %g %g\n",
446 rowx, rlo[rowx], rsol[rowx], rup[rowx]);
447 }
448 }
449 }
450 delete[] rsol;
451 }
452 #endif
tightenDoubletons2(CoinPresolveMatrix * prob)453 static int tightenDoubletons2(CoinPresolveMatrix *prob)
454 {
455 // column-major representation
456 const int ncols = prob->ncols_;
457 const CoinBigIndex *const mcstrt = prob->mcstrt_;
458 const int *const hincol = prob->hincol_;
459 const int *const hrow = prob->hrow_;
460 double *colels = prob->colels_;
461 double *cost = prob->cost_;
462
463 // column type, bounds, solution, and status
464 const unsigned char *const integerType = prob->integerType_;
465 double *clo = prob->clo_;
466 double *cup = prob->cup_;
467 // row-major representation
468 //const int nrows = prob->nrows_ ;
469 const CoinBigIndex *const mrstrt = prob->mrstrt_;
470 const int *const hinrow = prob->hinrow_;
471 const int *const hcol = prob->hcol_;
472 double *rowels = prob->rowels_;
473
474 // row bounds
475 double *const rlo = prob->rlo_;
476 double *const rup = prob->rup_;
477
478 // tolerances
479 //const double ekkinf2 = PRESOLVE_SMALL_INF ;
480 //const double ekkinf = ekkinf2*1.0e8 ;
481 //const double ztolcbarj = prob->ztoldj_ ;
482 //const CoinRelFltEq relEq(prob->ztolzb_) ;
483 int numberChanged = 0;
484 double bound[2];
485 double alpha[2] = { 0.0, 0.0 };
486 double offset = 0.0;
487
488 for (int icol = 0; icol < ncols; icol++) {
489 if (hincol[icol] == 2) {
490 CoinBigIndex start = mcstrt[icol];
491 int row0 = hrow[start];
492 if (hinrow[row0] != 2)
493 continue;
494 int row1 = hrow[start + 1];
495 if (hinrow[row1] != 2)
496 continue;
497 double element0 = colels[start];
498 double rowUpper0 = rup[row0];
499 bool swapSigns0 = false;
500 if (rlo[row0] > -1.0e30) {
501 if (rup[row0] > 1.0e30) {
502 swapSigns0 = true;
503 rowUpper0 = -rlo[row0];
504 element0 = -element0;
505 } else {
506 // range or equality
507 continue;
508 }
509 } else if (rup[row0] > 1.0e30) {
510 // free
511 continue;
512 }
513 #if 0
514 // skip here for speed
515 // skip if no cost (should be able to get rid of)
516 if (!cost[icol]) {
517 printf("should be able to get rid of %d with no cost\n",icol);
518 continue;
519 }
520 // skip if negative cost for now
521 if (cost[icol]<0.0) {
522 printf("code for negative cost\n");
523 continue;
524 }
525 #endif
526 double element1 = colels[start + 1];
527 double rowUpper1 = rup[row1];
528 bool swapSigns1 = false;
529 if (rlo[row1] > -1.0e30) {
530 if (rup[row1] > 1.0e30) {
531 swapSigns1 = true;
532 rowUpper1 = -rlo[row1];
533 element1 = -element1;
534 } else {
535 // range or equality
536 continue;
537 }
538 } else if (rup[row1] > 1.0e30) {
539 // free
540 continue;
541 }
542 double lowerX = clo[icol];
543 double upperX = cup[icol];
544 int otherCol = -1;
545 CoinBigIndex startRow = mrstrt[row0];
546 for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
547 int jcol = hcol[j];
548 if (jcol != icol) {
549 alpha[0] = swapSigns0 ? -rowels[j] : rowels[j];
550 otherCol = jcol;
551 }
552 }
553 startRow = mrstrt[row1];
554 bool possible = true;
555 for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
556 int jcol = hcol[j];
557 if (jcol != icol) {
558 if (jcol == otherCol) {
559 alpha[1] = swapSigns1 ? -rowels[j] : rowels[j];
560 } else {
561 possible = false;
562 }
563 }
564 }
565 if (possible) {
566 // skip if no cost (should be able to get rid of)
567 if (!cost[icol]) {
568 PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n", icol));
569 continue;
570 }
571 // skip if negative cost for now
572 if (cost[icol] < 0.0) {
573 PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
574 continue;
575 }
576 bound[0] = clo[otherCol];
577 bound[1] = cup[otherCol];
578 double lowestLowest = COIN_DBL_MAX;
579 double highestLowest = -COIN_DBL_MAX;
580 double lowestHighest = COIN_DBL_MAX;
581 double highestHighest = -COIN_DBL_MAX;
582 int binding0 = 0;
583 int binding1 = 0;
584 for (int k = 0; k < 2; k++) {
585 bool infLow0 = false;
586 bool infLow1 = false;
587 double sum0 = 0.0;
588 double sum1 = 0.0;
589 double value = bound[k];
590 if (fabs(value) < 1.0e30) {
591 sum0 += alpha[0] * value;
592 sum1 += alpha[1] * value;
593 } else {
594 if (alpha[0] > 0.0) {
595 if (value < 0.0)
596 infLow0 = true;
597 } else if (alpha[0] < 0.0) {
598 if (value > 0.0)
599 infLow0 = true;
600 }
601 if (alpha[1] > 0.0) {
602 if (value < 0.0)
603 infLow1 = true;
604 } else if (alpha[1] < 0.0) {
605 if (value > 0.0)
606 infLow1 = true;
607 }
608 }
609 /* Got sums
610 */
611 double thisLowest0 = -COIN_DBL_MAX;
612 double thisHighest0 = COIN_DBL_MAX;
613 if (element0 > 0.0) {
614 // upper bound unless inf&2 !=0
615 if (!infLow0)
616 thisHighest0 = (rowUpper0 - sum0) / element0;
617 } else {
618 // lower bound unless inf&2 !=0
619 if (!infLow0)
620 thisLowest0 = (rowUpper0 - sum0) / element0;
621 }
622 double thisLowest1 = -COIN_DBL_MAX;
623 double thisHighest1 = COIN_DBL_MAX;
624 if (element1 > 0.0) {
625 // upper bound unless inf&2 !=0
626 if (!infLow1)
627 thisHighest1 = (rowUpper1 - sum1) / element1;
628 } else {
629 // lower bound unless inf&2 !=0
630 if (!infLow1)
631 thisLowest1 = (rowUpper1 - sum1) / element1;
632 }
633 if (thisLowest0 > thisLowest1 + 1.0e-12) {
634 if (thisLowest0 > lowerX + 1.0e-12)
635 binding0 |= 1 << k;
636 } else if (thisLowest1 > thisLowest0 + 1.0e-12) {
637 if (thisLowest1 > lowerX + 1.0e-12)
638 binding1 |= 1 << k;
639 thisLowest0 = thisLowest1;
640 }
641 if (thisHighest0 < thisHighest1 - 1.0e-12) {
642 if (thisHighest0 < upperX - 1.0e-12)
643 binding0 |= 1 << k;
644 } else if (thisHighest1 < thisHighest0 - 1.0e-12) {
645 if (thisHighest1 < upperX - 1.0e-12)
646 binding1 |= 1 << k;
647 thisHighest0 = thisHighest1;
648 }
649 lowestLowest = CoinMin(lowestLowest, thisLowest0);
650 highestHighest = CoinMax(highestHighest, thisHighest0);
651 lowestHighest = CoinMin(lowestHighest, thisHighest0);
652 highestLowest = CoinMax(highestLowest, thisLowest0);
653 }
654 // see if any good
655 //#define PRINT_VALUES
656 if (!binding0 || !binding1) {
657 PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n", icol));
658 } else {
659 #ifdef PRINT_VALUES
660 printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
661 icol, lowerX, upperX, lowestLowest, lowestHighest,
662 highestLowest, highestHighest);
663 #endif
664 // if integer adjust
665 if (integerType[icol]) {
666 lowestLowest = ceil(lowestLowest - 1.0e-5);
667 highestLowest = ceil(highestLowest - 1.0e-5);
668 lowestHighest = floor(lowestHighest + 1.0e-5);
669 highestHighest = floor(highestHighest + 1.0e-5);
670 }
671 // if costed may be able to adjust
672 if (cost[icol] >= 0.0) {
673 if (highestLowest < upperX && highestLowest >= lowerX && highestHighest < 1.0e30) {
674 highestHighest = CoinMin(highestHighest, highestLowest);
675 }
676 }
677 if (cost[icol] <= 0.0) {
678 if (lowestHighest > lowerX && lowestHighest <= upperX && lowestHighest > -1.0e30) {
679 lowestLowest = CoinMax(lowestLowest, lowestHighest);
680 }
681 }
682 #if 1
683 if (lowestLowest > lowerX + 1.0e-8) {
684 #ifdef PRINT_VALUES
685 printf("Can increase lower bound on %d from %g to %g\n",
686 icol, lowerX, lowestLowest);
687 #endif
688 lowerX = lowestLowest;
689 }
690 if (highestHighest < upperX - 1.0e-8) {
691 #ifdef PRINT_VALUES
692 printf("Can decrease upper bound on %d from %g to %g\n",
693 icol, upperX, highestHighest);
694 #endif
695 upperX = highestHighest;
696 }
697 #endif
698 // see if we can move costs
699 double xValue;
700 double yValue0;
701 double yValue1;
702 double newLower = COIN_DBL_MAX;
703 double newUpper = -COIN_DBL_MAX;
704 #ifdef PRINT_VALUES
705 double ranges0[2];
706 double ranges1[2];
707 #endif
708 double costEqual;
709 double slope[2];
710 assert(binding0 + binding1 == 3);
711 // get where equal
712 xValue = (rowUpper0 * element1 - rowUpper1 * element0) / (alpha[0] * element1 - alpha[1] * element0);
713 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
714 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
715 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
716 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
717 double xValueEqual = xValue;
718 double yValueEqual = yValue0;
719 costEqual = xValue * cost[otherCol] + yValueEqual * cost[icol];
720 if (binding0 == 1) {
721 #ifdef PRINT_VALUES
722 ranges0[0] = bound[0];
723 ranges0[1] = yValue0;
724 ranges1[0] = yValue0;
725 ranges1[1] = bound[1];
726 #endif
727 // take x 1.0 down
728 double x = xValue - 1.0;
729 double y = (rowUpper0 - x * alpha[0]) / element0;
730 double costTotal = x * cost[otherCol] + y * cost[icol];
731 slope[0] = costEqual - costTotal;
732 // take x 1.0 up
733 x = xValue + 1.0;
734 y = (rowUpper1 - x * alpha[1]) / element0;
735 costTotal = x * cost[otherCol] + y * cost[icol];
736 slope[1] = costTotal - costEqual;
737 } else {
738 #ifdef PRINT_VALUES
739 ranges1[0] = bound[0];
740 ranges1[1] = yValue0;
741 ranges0[0] = yValue0;
742 ranges0[1] = bound[1];
743 #endif
744 // take x 1.0 down
745 double x = xValue - 1.0;
746 double y = (rowUpper1 - x * alpha[1]) / element0;
747 double costTotal = x * cost[otherCol] + y * cost[icol];
748 slope[1] = costEqual - costTotal;
749 // take x 1.0 up
750 x = xValue + 1.0;
751 y = (rowUpper0 - x * alpha[0]) / element0;
752 costTotal = x * cost[otherCol] + y * cost[icol];
753 slope[0] = costTotal - costEqual;
754 }
755 #ifdef PRINT_VALUES
756 printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
757 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
758 printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
759 costEqual, ranges0[0], ranges0[1], slope[0], ranges1[0], ranges1[1], slope[1]);
760 #endif
761 xValue = bound[0];
762 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
763 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
764 #ifdef PRINT_VALUES
765 printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
766 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
767 #endif
768 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
769 // cost>0 so will be at lower
770 //double yValueAtBound0=newLower;
771 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
772 xValue = bound[1];
773 yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
774 yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
775 #ifdef PRINT_VALUES
776 printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
777 otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
778 #endif
779 newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
780 // cost>0 so will be at lower
781 //double yValueAtBound1=newLower;
782 newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
783 lowerX = CoinMax(lowerX, newLower - 1.0e-12 * fabs(newLower));
784 upperX = CoinMin(upperX, newUpper + 1.0e-12 * fabs(newUpper));
785 // Now make duplicate row
786 // keep row 0 so need to adjust costs so same
787 #ifdef PRINT_VALUES
788 printf("Costs for x %g,%g,%g are %g,%g,%g\n",
789 xValueEqual - 1.0, xValueEqual, xValueEqual + 1.0,
790 costEqual - slope[0], costEqual, costEqual + slope[1]);
791 #endif
792 double costOther = cost[otherCol] + slope[1];
793 double costThis = cost[icol] + slope[1] * (element0 / alpha[0]);
794 xValue = xValueEqual;
795 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
796 double thisOffset = costEqual - (costOther * xValue + costThis * yValue0);
797 offset += thisOffset;
798 #ifdef PRINT_VALUES
799 printf("new cost at equal %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
800 #endif
801 xValue = xValueEqual - 1.0;
802 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
803 #ifdef PRINT_VALUES
804 printf("new cost at -1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
805 #endif
806 assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual - slope[0])) < 1.0e-5);
807 xValue = xValueEqual + 1.0;
808 yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
809 #ifdef PRINT_VALUES
810 printf("new cost at +1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
811 #endif
812 assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual + slope[1])) < 1.0e-5);
813 numberChanged++;
814 // continue;
815 cost[otherCol] = costOther;
816 cost[icol] = costThis;
817 clo[icol] = lowerX;
818 cup[icol] = upperX;
819 CoinBigIndex startCol[2];
820 CoinBigIndex endCol[2];
821 startCol[0] = mcstrt[icol];
822 endCol[0] = startCol[0] + 2;
823 startCol[1] = mcstrt[otherCol];
824 endCol[1] = startCol[1] + hincol[otherCol];
825 double values[2] = { 0.0, 0.0 };
826 for (int k = 0; k < 2; k++) {
827 for (CoinBigIndex i = startCol[k]; i < endCol[k]; i++) {
828 if (hrow[i] == row0)
829 values[k] = colels[i];
830 }
831 for (CoinBigIndex i = startCol[k]; i < endCol[k]; i++) {
832 if (hrow[i] == row1)
833 colels[i] = values[k];
834 }
835 }
836 for (CoinBigIndex i = mrstrt[row1]; i < mrstrt[row1] + 2; i++) {
837 if (hcol[i] == icol)
838 rowels[i] = values[0];
839 else
840 rowels[i] = values[1];
841 }
842 }
843 }
844 }
845 }
846 #if ABC_NORMAL_DEBUG > 0
847 if (offset)
848 printf("Cost offset %g\n", offset);
849 #endif
850 return numberChanged;
851 }
852 //#define COIN_PRESOLVE_BUG
853 #ifdef COIN_PRESOLVE_BUG
854 static int counter = 1000000;
855 static int startEmptyRows = 0;
856 static int startEmptyColumns = 0;
break2(CoinPresolveMatrix * prob)857 static bool break2(CoinPresolveMatrix *prob)
858 {
859 int droppedRows = prob->countEmptyRows() - startEmptyRows;
860 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
861 startEmptyRows = prob->countEmptyRows();
862 startEmptyColumns = prob->countEmptyCols();
863 printf("Dropped %d rows and %d columns - current empty %d, %d\n", droppedRows,
864 droppedColumns, startEmptyRows, startEmptyColumns);
865 counter--;
866 if (!counter) {
867 printf("skipping next and all\n");
868 }
869 return (counter <= 0);
870 }
871 #define possibleBreak \
872 if (break2(prob)) \
873 break
874 #define possibleSkip if (!break2(prob))
875 #else
876 #define possibleBreak
877 #define possibleSkip
878 #endif
879 #define SOME_PRESOLVE_DETAIL
880 #ifndef SOME_PRESOLVE_DETAIL
881 #define printProgress(x, y) \
882 { \
883 }
884 #else
885 #define printProgress(x, y) \
886 { \
887 if ((presolveActions_ & 0x80000000) != 0) \
888 printf("%c loop %d %d empty rows, %d empty columns\n", x, y, prob->countEmptyRows(), \
889 prob->countEmptyCols()); \
890 }
891 #endif
892 // This is the presolve loop.
893 // It is a separate virtual function so that it can be easily
894 // customized by subclassing CoinPresolve.
presolve(CoinPresolveMatrix * prob)895 const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
896 {
897 // Messages
898 CoinMessages messages = CoinMessage(prob->messages().language());
899 paction_ = 0;
900 prob->maxSubstLevel_ = CoinMax(3, prob->maxSubstLevel_);
901 #ifndef PRESOLVE_DETAIL
902 if (prob->tuning_) {
903 #endif
904 int numberEmptyRows = 0;
905 for (int i = 0; i < prob->nrows_; i++) {
906 if (!prob->hinrow_[i]) {
907 PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n", i));
908 //printf("pre_empty row %d\n",i);
909 numberEmptyRows++;
910 }
911 }
912 int numberEmptyCols = 0;
913 for (int i = 0; i < prob->ncols_; i++) {
914 if (!prob->hincol_[i]) {
915 PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n", i));
916 //printf("pre_empty col %d\n",i);
917 numberEmptyCols++;
918 }
919 }
920 printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
921 numberEmptyRows, numberEmptyCols);
922 #ifndef PRESOLVE_DETAIL
923 }
924 #endif
925 prob->status_ = 0; // say feasible
926 printProgress('A', 0);
927 paction_ = make_fixed(prob, paction_);
928 paction_ = testRedundant(prob, paction_);
929 printProgress('B', 0);
930 // if integers then switch off dual stuff
931 // later just do individually
932 bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
933 // but allow in some cases
934 if ((presolveActions_ & 512) != 0)
935 doDualStuff = true;
936 if (prob->anyProhibited())
937 doDualStuff = false;
938 if (!doDual())
939 doDualStuff = false;
940 #if PRESOLVE_CONSISTENCY
941 // presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
942 presolve_links_ok(prob, false, true);
943 #endif
944
945 if (!prob->status_) {
946 bool slackSingleton = doSingletonColumn();
947 slackSingleton = true;
948 const bool slackd = doSingleton();
949 const bool doubleton = doDoubleton();
950 const bool tripleton = doTripleton();
951 //#define NO_FORCING
952 #ifndef NO_FORCING
953 const bool forcing = doForcing();
954 #endif
955 const bool ifree = doImpliedFree();
956 const bool zerocost = doTighten();
957 const bool dupcol = doDupcol();
958 const bool duprow = doDuprow();
959 const bool dual = doDualStuff;
960 // Whether we want to allow duplicate intersections
961 if (doIntersection())
962 prob->presolveOptions_ |= 0x10;
963 // zero small elements in aggregation
964 prob->presolveOptions_ |= zeroSmall() * 0x20000;
965 // some things are expensive so just do once (normally)
966
967 int i;
968 // say look at all
969 if (!prob->anyProhibited()) {
970 for (i = 0; i < nrows_; i++)
971 prob->rowsToDo_[i] = i;
972 prob->numberRowsToDo_ = nrows_;
973 for (i = 0; i < ncols_; i++)
974 prob->colsToDo_[i] = i;
975 prob->numberColsToDo_ = ncols_;
976 } else {
977 // some stuff must be left alone
978 prob->numberRowsToDo_ = 0;
979 for (i = 0; i < nrows_; i++)
980 if (!prob->rowProhibited(i))
981 prob->rowsToDo_[prob->numberRowsToDo_++] = i;
982 prob->numberColsToDo_ = 0;
983 for (i = 0; i < ncols_; i++)
984 if (!prob->colProhibited(i))
985 prob->colsToDo_[prob->numberColsToDo_++] = i;
986 }
987
988 // transfer costs (may want to do it in OsiPresolve)
989 // need a transfer back at end of postsolve transferCosts(prob);
990
991 int iLoop;
992 #if PRESOLVE_CHECK_SOL
993 check_sol(prob, 1.0e0);
994 #endif
995 if (dupcol) {
996 // maybe allow integer columns to be checked
997 if ((presolveActions_ & 512) != 0)
998 prob->setPresolveOptions(prob->presolveOptions() | 1);
999 possibleSkip;
1000 paction_ = dupcol_action::presolve(prob, paction_);
1001 printProgress('C', 0);
1002 }
1003 if (doTwoxTwo()) {
1004 possibleSkip;
1005 paction_ = twoxtwo_action::presolve(prob, paction_);
1006 }
1007 if (duprow) {
1008 possibleSkip;
1009 if (doTwoxTwo()) {
1010 int nTightened = tightenDoubletons2(prob);
1011 if (nTightened)
1012 PRESOLVE_DETAIL_PRINT(printf("%d doubletons tightened\n",
1013 nTightened));
1014 }
1015 paction_ = duprow_action::presolve(prob, paction_);
1016 printProgress('D', 0);
1017 //paction_ = doubleton_action::presolve(prob, paction_);
1018 //printProgress('d',0);
1019 //if (doDependency()) {
1020 //paction_ = duprow3_action::presolve(prob, paction_);
1021 //printProgress('Z',0);
1022 //}
1023 }
1024 if (doGubrow()) {
1025 possibleSkip;
1026 paction_ = gubrow_action::presolve(prob, paction_);
1027 printProgress('E', 0);
1028 }
1029 if (ifree) {
1030 int fill_level = CoinMax(10, prob->maxSubstLevel_);
1031 const CoinPresolveAction *lastAction = NULL;
1032 int iPass = 4;
1033 while (lastAction != paction_ && iPass) {
1034 lastAction = paction_;
1035 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1036 printProgress('l', 0);
1037 iPass--;
1038 }
1039 }
1040
1041 if ((presolveActions_ & 16384) != 0)
1042 prob->setPresolveOptions(prob->presolveOptions() | 16384);
1043 // For inaccurate data in implied free
1044 if ((presolveActions_ & 1024) != 0)
1045 prob->setPresolveOptions(prob->presolveOptions() | 0x20000);
1046 // Check number rows dropped
1047 int lastDropped = 0;
1048 prob->pass_ = 0;
1049 #if CLP_INHERIT_MODE > 1
1050 int numberRowsStart = nrows_ - prob->countEmptyRows();
1051 int numberColumnsStart = ncols_ - prob->countEmptyCols();
1052 int numberRowsLeft = numberRowsStart;
1053 int numberColumnsLeft = numberColumnsStart;
1054 bool lastPassWasGood = true;
1055 #if ABC_NORMAL_DEBUG
1056 printf("Original rows,columns %d,%d starting first pass with %d,%d\n",
1057 nrows_, ncols_, numberRowsLeft, numberColumnsLeft);
1058 #endif
1059 #endif
1060 if (numberPasses_ <= 5)
1061 prob->presolveOptions_ |= 0x10000; // say more lightweight
1062 for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
1063 // See if we want statistics
1064 if ((presolveActions_ & 0x80000000) != 0)
1065 printf("Starting major pass %d after %g seconds with %d rows, %d columns\n", iLoop + 1, CoinCpuTime() - prob->startTime_,
1066 nrows_ - prob->countEmptyRows(),
1067 ncols_ - prob->countEmptyCols());
1068 #ifdef PRESOLVE_SUMMARY
1069 printf("Starting major pass %d\n", iLoop + 1);
1070 #endif
1071 const CoinPresolveAction *const paction0 = paction_;
1072 // look for substitutions with no fill
1073 //#define IMPLIED 3
1074 #ifdef IMPLIED
1075 int fill_level = 3;
1076 #define IMPLIED2 99
1077 #if IMPLIED != 3
1078 #if IMPLIED > 2 && IMPLIED < 11
1079 fill_level = IMPLIED;
1080 COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1081 #endif
1082 #if IMPLIED > 11 && IMPLIED < 21
1083 fill_level = -(IMPLIED - 10);
1084 COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
1085 #endif
1086 #endif
1087 #else
1088 int fill_level = prob->maxSubstLevel_;
1089 #endif
1090 int whichPass = 0;
1091 while (1) {
1092 whichPass++;
1093 prob->pass_++;
1094 const CoinPresolveAction *const paction1 = paction_;
1095
1096 if (slackd) {
1097 bool notFinished = true;
1098 while (notFinished) {
1099 possibleBreak;
1100 paction_ = slack_doubleton_action::presolve(prob, paction_,
1101 notFinished);
1102 }
1103 printProgress('F', iLoop + 1);
1104 if (prob->status_)
1105 break;
1106 }
1107 if (zerocost) {
1108 possibleBreak;
1109 paction_ = do_tighten_action::presolve(prob, paction_);
1110 if (prob->status_)
1111 break;
1112 printProgress('J', iLoop + 1);
1113 }
1114 if (dual && whichPass == 1) {
1115 // this can also make E rows so do one bit here
1116 possibleBreak;
1117 paction_ = remove_dual_action::presolve(prob, paction_);
1118 if (prob->status_)
1119 break;
1120 printProgress('G', iLoop + 1);
1121 }
1122
1123 if (doubleton) {
1124 possibleBreak;
1125 paction_ = doubleton_action::presolve(prob, paction_);
1126 if (prob->status_)
1127 break;
1128 printProgress('H', iLoop + 1);
1129 }
1130 if (tripleton) {
1131 possibleBreak;
1132 paction_ = tripleton_action::presolve(prob, paction_);
1133 if (prob->status_)
1134 break;
1135 printProgress('I', iLoop + 1);
1136 }
1137
1138 #ifndef NO_FORCING
1139 if (forcing) {
1140 possibleBreak;
1141 paction_ = forcing_constraint_action::presolve(prob, paction_);
1142 if (prob->status_)
1143 break;
1144 printProgress('K', iLoop + 1);
1145 }
1146 #endif
1147
1148 if (ifree && (whichPass % 5) == 1) {
1149 possibleBreak;
1150 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1151 if (prob->status_)
1152 break;
1153 printProgress('L', iLoop + 1);
1154 }
1155
1156 #if PRESOLVE_CHECK_SOL
1157 check_sol(prob, 1.0e0);
1158 #endif
1159
1160 #if PRESOLVE_CONSISTENCY
1161 // presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
1162 // prob->nrows_);
1163 presolve_links_ok(prob, false, true);
1164 #endif
1165
1166 //#if PRESOLVE_DEBUG
1167 // presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
1168 // prob->ncols_);
1169 //#endif
1170 //#if PRESOLVE_CONSISTENCY
1171 // prob->consistent();
1172 //#endif
1173 #if PRESOLVE_CONSISTENCY
1174 presolve_no_zeros(prob, true, false);
1175 presolve_consistent(prob, true);
1176 #endif
1177
1178 {
1179 // set up for next pass
1180 // later do faster if many changes i.e. memset and memcpy
1181 const int *count = prob->hinrow_;
1182 const int *nextToDo = prob->nextRowsToDo_;
1183 int *toDo = prob->rowsToDo_;
1184 int nNext = prob->numberNextRowsToDo_;
1185 int n = 0;
1186 for (int i = 0; i < nNext; i++) {
1187 int index = nextToDo[i];
1188 prob->unsetRowChanged(index);
1189 if (count[index])
1190 toDo[n++] = index;
1191 }
1192 prob->numberRowsToDo_ = n;
1193 prob->numberNextRowsToDo_ = 0;
1194 count = prob->hincol_;
1195 nextToDo = prob->nextColsToDo_;
1196 toDo = prob->colsToDo_;
1197 nNext = prob->numberNextColsToDo_;
1198 n = 0;
1199 for (int i = 0; i < nNext; i++) {
1200 int index = nextToDo[i];
1201 prob->unsetColChanged(index);
1202 if (count[index])
1203 toDo[n++] = index;
1204 }
1205 prob->numberColsToDo_ = n;
1206 prob->numberNextColsToDo_ = 0;
1207 }
1208 if (paction_ == paction1 && fill_level > 0)
1209 break;
1210 }
1211 // say look at all
1212 int i;
1213 if (!prob->anyProhibited()) {
1214 const int *count = prob->hinrow_;
1215 int *toDo = prob->rowsToDo_;
1216 int n = 0;
1217 for (int i = 0; i < nrows_; i++) {
1218 prob->unsetRowChanged(i);
1219 if (count[i])
1220 toDo[n++] = i;
1221 }
1222 prob->numberRowsToDo_ = n;
1223 prob->numberNextRowsToDo_ = 0;
1224 count = prob->hincol_;
1225 toDo = prob->colsToDo_;
1226 n = 0;
1227 for (int i = 0; i < ncols_; i++) {
1228 prob->unsetColChanged(i);
1229 if (count[i])
1230 toDo[n++] = i;
1231 }
1232 prob->numberColsToDo_ = n;
1233 prob->numberNextColsToDo_ = 0;
1234 } else {
1235 // some stuff must be left alone
1236 prob->numberRowsToDo_ = 0;
1237 for (i = 0; i < nrows_; i++)
1238 if (!prob->rowProhibited(i))
1239 prob->rowsToDo_[prob->numberRowsToDo_++] = i;
1240 prob->numberColsToDo_ = 0;
1241 for (i = 0; i < ncols_; i++)
1242 if (!prob->colProhibited(i))
1243 prob->colsToDo_[prob->numberColsToDo_++] = i;
1244 }
1245 // now expensive things
1246 // this caused world.mps to run into numerical difficulties
1247 #ifdef PRESOLVE_SUMMARY
1248 printf("Starting expensive\n");
1249 #endif
1250
1251 if (dual) {
1252 int itry;
1253 for (itry = 0; itry < 5; itry++) {
1254 possibleBreak;
1255 paction_ = remove_dual_action::presolve(prob, paction_);
1256 if (prob->status_)
1257 break;
1258 printProgress('M', iLoop + 1);
1259 const CoinPresolveAction *const paction2 = paction_;
1260 if (ifree) {
1261 #ifdef IMPLIED
1262 #if IMPLIED2 == 0
1263 int fill_level = 0; // switches off substitution
1264 #elif IMPLIED2 != 99
1265 int fill_level = IMPLIED2;
1266 #endif
1267 #endif
1268 if ((itry & 1) == 0) {
1269 possibleBreak;
1270 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1271 }
1272 if (prob->status_)
1273 break;
1274 printProgress('N', iLoop + 1);
1275 }
1276 if (paction_ == paction2)
1277 break;
1278 }
1279 } else if (ifree) {
1280 // just free
1281 #ifdef IMPLIED
1282 #if IMPLIED2 == 0
1283 int fill_level = 0; // switches off substitution
1284 #elif IMPLIED2 != 99
1285 int fill_level = IMPLIED2;
1286 #endif
1287 #endif
1288 possibleBreak;
1289 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
1290 if (prob->status_)
1291 break;
1292 printProgress('O', iLoop + 1);
1293 }
1294 #if PRESOLVE_CHECK_SOL
1295 check_sol(prob, 1.0e0);
1296 #endif
1297 if (dupcol) {
1298 // maybe allow integer columns to be checked
1299 if ((presolveActions_ & 512) != 0)
1300 prob->setPresolveOptions(prob->presolveOptions() | 1);
1301 possibleBreak;
1302 paction_ = dupcol_action::presolve(prob, paction_);
1303 if (prob->status_)
1304 break;
1305 printProgress('P', iLoop + 1);
1306 }
1307 #if PRESOLVE_CHECK_SOL
1308 check_sol(prob, 1.0e0);
1309 #endif
1310
1311 if (duprow) {
1312 possibleBreak;
1313 paction_ = duprow_action::presolve(prob, paction_);
1314 if (prob->status_)
1315 break;
1316 printProgress('Q', iLoop + 1);
1317 }
1318 // Marginally slower on netlib if this call is enabled.
1319 // paction_ = testRedundant(prob,paction_) ;
1320 #if PRESOLVE_CHECK_SOL
1321 check_sol(prob, 1.0e0);
1322 #endif
1323 bool stopLoop = false;
1324 {
1325 int *hinrow = prob->hinrow_;
1326 int numberDropped = 0;
1327 for (i = 0; i < nrows_; i++)
1328 if (!hinrow[i])
1329 numberDropped++;
1330
1331 prob->messageHandler()->message(COIN_PRESOLVE_PASS,
1332 messages)
1333 << numberDropped << iLoop + 1
1334 << CoinMessageEol;
1335 //printf("%d rows dropped after pass %d\n",numberDropped,
1336 // iLoop+1);
1337 if (numberDropped == lastDropped)
1338 stopLoop = true;
1339 else
1340 lastDropped = numberDropped;
1341 }
1342 // Do this here as not very loopy
1343 if (slackSingleton) {
1344 // On most passes do not touch costed slacks
1345 if (paction_ != paction0 && !stopLoop) {
1346 possibleBreak;
1347 paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
1348 } else {
1349 // do costed if Clp (at end as ruins rest of presolve)
1350 possibleBreak;
1351 #ifndef CLP_MOVE_COSTS
1352 paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
1353 #else
1354 double *fakeRowObjective = new double[prob->nrows_];
1355 memset(fakeRowObjective, 0, prob->nrows_ * sizeof(double));
1356 paction_ = slack_singleton_action::presolve(prob, paction_, fakeRowObjective);
1357 delete[] fakeRowObjective;
1358 #endif
1359 stopLoop = true;
1360 }
1361 printProgress('R', iLoop + 1);
1362 }
1363 #if PRESOLVE_CHECK_SOL
1364 check_sol(prob, 1.0e0);
1365 #endif
1366 if (paction_ == paction0 || stopLoop)
1367 break;
1368 #if CLP_INHERIT_MODE > 1
1369 // see whether to stop anyway
1370 int numberRowsNow = nrows_ - prob->countEmptyRows();
1371 int numberColumnsNow = ncols_ - prob->countEmptyCols();
1372 #if ABC_NORMAL_DEBUG
1373 printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n",
1374 nrows_, ncols_, numberRowsLeft, numberColumnsLeft, iLoop + 1, numberRowsNow,
1375 numberColumnsNow);
1376 #endif
1377 int rowsDeleted = numberRowsLeft - numberRowsNow;
1378 int columnsDeleted = numberColumnsLeft - numberColumnsNow;
1379 if (iLoop > 15) {
1380 if (rowsDeleted * 100 < numberRowsStart && columnsDeleted * 100 < numberColumnsStart)
1381 break;
1382 lastPassWasGood = true;
1383 } else if (rowsDeleted * 100 < numberRowsStart && rowsDeleted < 500 && columnsDeleted * 100 < numberColumnsStart && columnsDeleted < 500) {
1384 if (!lastPassWasGood)
1385 break;
1386 else
1387 lastPassWasGood = false;
1388 } else {
1389 lastPassWasGood = true;
1390 }
1391 numberRowsLeft = numberRowsNow;
1392 numberColumnsLeft = numberColumnsNow;
1393 #endif
1394 }
1395 }
1396 if (!prob->status_ && doDependency()) {
1397 paction_ = duprow3_action::presolve(prob, paction_);
1398 printProgress('Z', 0);
1399 }
1400 prob->presolveOptions_ &= ~0x10000;
1401 if (!prob->status_) {
1402 paction_ = drop_zero_coefficients(prob, paction_);
1403 #if PRESOLVE_CHECK_SOL
1404 check_sol(prob, 1.0e0);
1405 #endif
1406
1407 paction_ = drop_empty_cols_action::presolve(prob, paction_);
1408 paction_ = drop_empty_rows_action::presolve(prob, paction_);
1409 #if PRESOLVE_CHECK_SOL
1410 check_sol(prob, 1.0e0);
1411 #endif
1412 }
1413
1414 if (prob->status_) {
1415 if (prob->status_ == 1)
1416 prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
1417 messages)
1418 << prob->feasibilityTolerance_
1419 << CoinMessageEol;
1420 else if (prob->status_ == 2)
1421 prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
1422 messages)
1423 << CoinMessageEol;
1424 else
1425 prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
1426 messages)
1427 << CoinMessageEol;
1428 // get rid of data
1429 destroyPresolve();
1430 }
1431 return (paction_);
1432 }
1433
1434 void check_djs(CoinPostsolveMatrix *prob);
1435
1436 // We could have implemented this by having each postsolve routine
1437 // directly call the next one, but this may make it easier to add debugging checks.
postsolve(CoinPostsolveMatrix & prob)1438 void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
1439 {
1440 {
1441 // Check activities
1442 double *colels = prob.colels_;
1443 int *hrow = prob.hrow_;
1444 CoinBigIndex *mcstrt = prob.mcstrt_;
1445 int *hincol = prob.hincol_;
1446 CoinBigIndex *link = prob.link_;
1447 int ncols = prob.ncols_;
1448
1449 char *cdone = prob.cdone_;
1450
1451 double *csol = prob.sol_;
1452 int nrows = prob.nrows_;
1453
1454 int colx;
1455
1456 double *rsol = prob.acts_;
1457 memset(rsol, 0, nrows * sizeof(double));
1458
1459 for (colx = 0; colx < ncols; ++colx) {
1460 if (cdone[colx]) {
1461 CoinBigIndex k = mcstrt[colx];
1462 int nx = hincol[colx];
1463 double solutionValue = csol[colx];
1464 for (int i = 0; i < nx; ++i) {
1465 int row = hrow[k];
1466 double coeff = colels[k];
1467 k = link[k];
1468 assert(k != NO_LINK || i == nx - 1);
1469 rsol[row] += solutionValue * coeff;
1470 }
1471 }
1472 }
1473 }
1474 if (prob.maxmin_ < 0) {
1475 //for (int i=0;i<presolvedModel_->numberRows();i++)
1476 //prob.rowduals_[i]=-prob.rowduals_[i];
1477 for (int i = 0; i < ncols_; i++) {
1478 prob.cost_[i] = -prob.cost_[i];
1479 //prob.rcosts_[i]=-prob.rcosts_[i];
1480 }
1481 prob.maxmin_ = 1.0;
1482 }
1483 const CoinPresolveAction *paction = paction_;
1484 //#define PRESOLVE_DEBUG 1
1485 #if PRESOLVE_DEBUG
1486 // Check only works after first one
1487 int checkit = -1;
1488 #endif
1489
1490 while (paction) {
1491 #if PRESOLVE_DEBUG
1492 printf("POSTSOLVING %s\n", paction->name());
1493 #endif
1494
1495 paction->postsolve(&prob);
1496
1497 #if PRESOLVE_DEBUG
1498 #if 0
1499 /*
1500 This check fails (on exmip1 (!) in osiUnitTest) because clp
1501 enters postsolve with a solution that seems to have incorrect
1502 status for a logical. You can see similar behaviour with
1503 column status --- incorrect entering postsolve.
1504 -- lh, 111207 --
1505 */
1506 {
1507 int nr = 0;
1508 int i;
1509 for (i = 0; i < prob.nrows_; i++) {
1510 if ((prob.rowstat_[i] & 7) == 1) {
1511 nr++;
1512 } else if ((prob.rowstat_[i] & 7) == 2) {
1513 // at ub
1514 assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
1515 } else if ((prob.rowstat_[i] & 7) == 3) {
1516 // at lb
1517 assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
1518 }
1519 }
1520 int nc = 0;
1521 for (i = 0; i < prob.ncols_; i++) {
1522 if ((prob.colstat_[i] & 7) == 1)
1523 nc++;
1524 }
1525 printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
1526 }
1527 #endif // if 0
1528 checkit++;
1529 if (prob.colstat_ && checkit > 0) {
1530 presolve_check_nbasic(&prob);
1531 presolve_check_sol(&prob, 2, 2, 1);
1532 }
1533 #endif
1534 paction = paction->next;
1535 #if PRESOLVE_DEBUG
1536 // check_djs(&prob);
1537 if (checkit > 0)
1538 presolve_check_reduced_costs(&prob);
1539 #endif
1540 }
1541 #if PRESOLVE_DEBUG
1542 if (prob.colstat_) {
1543 presolve_check_nbasic(&prob);
1544 presolve_check_sol(&prob, 2, 2, 1);
1545 }
1546 #endif
1547 #undef PRESOLVE_DEBUG
1548
1549 #if 0 && PRESOLVE_DEBUG
1550 for (i = 0; i < ncols0; i++) {
1551 if (!cdone[i]) {
1552 printf("!cdone[%d]\n", i);
1553 abort();
1554 }
1555 }
1556
1557 for (i = 0; i < nrows0; i++) {
1558 if (!rdone[i]) {
1559 printf("!rdone[%d]\n", i);
1560 abort();
1561 }
1562 }
1563
1564
1565 for (i = 0; i < ncols0; i++) {
1566 if (sol[i] < -1e10 || sol[i] > 1e10)
1567 printf("!!!%d %g\n", i, sol[i]);
1568
1569 }
1570
1571 #endif
1572
1573 #if 0 && PRESOLVE_DEBUG
1574 // debug check: make sure we ended up with same original matrix
1575 {
1576 int identical = 1;
1577
1578 for (int i = 0; i < ncols0; i++) {
1579 PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
1580 CoinBigIndex kcs0 = &prob->mcstrt0[i];
1581 CoinBigIndex kcs = mcstrt[i];
1582 int n = hincol[i];
1583 for (int k = 0; k < n; k++) {
1584 CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
1585
1586 if (k1 == kcs + n) {
1587 printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
1588 abort();
1589 }
1590
1591 if (colels[k1] != &prob->dels0[kcs0+k])
1592 printf("BAD COLEL[%d %d %d]: %g\n",
1593 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
1594
1595 if (kcs0 + k != k1)
1596 identical = 0;
1597 }
1598 }
1599 printf("identical? %d\n", identical);
1600 }
1601 #endif
1602 }
1603
getTolerance(const ClpSimplex * si,ClpDblParam key)1604 static inline double getTolerance(const ClpSimplex *si, ClpDblParam key)
1605 {
1606 double tol;
1607 if (!si->getDblParam(key, tol)) {
1608 CoinPresolveAction::throwCoinError("getDblParam failed",
1609 "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
1610 }
1611 return (tol);
1612 }
1613
1614 // Assumptions:
1615 // 1. nrows>=m.getNumRows()
1616 // 2. ncols>=m.getNumCols()
1617 //
1618 // In presolve, these values are equal.
1619 // In postsolve, they may be inequal, since the reduced problem
1620 // may be smaller, but we need room for the large problem.
1621 // ncols may be larger than si.getNumCols() in postsolve,
1622 // this at that point si will be the reduced problem,
1623 // but we need to reserve enough space for the original problem.
CoinPrePostsolveMatrix(const ClpSimplex * si,int ncols_in,int nrows_in,CoinBigIndex nelems_in,double bulkRatio)1624 CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex *si,
1625 int ncols_in,
1626 int nrows_in,
1627 CoinBigIndex nelems_in,
1628 double bulkRatio)
1629 : ncols_(si->getNumCols())
1630 , nrows_(si->getNumRows())
1631 , nelems_(si->getNumElements())
1632 , ncols0_(ncols_in)
1633 , nrows0_(nrows_in)
1634 , bulkRatio_(bulkRatio)
1635 , mcstrt_(new CoinBigIndex[ncols_in + 1])
1636 , hincol_(new int[ncols_in + 1])
1637 , cost_(new double[ncols_in])
1638 , clo_(new double[ncols_in])
1639 , cup_(new double[ncols_in])
1640 , rlo_(new double[nrows_in])
1641 , rup_(new double[nrows_in])
1642 , originalColumn_(new int[ncols_in])
1643 , originalRow_(new int[nrows_in])
1644 , ztolzb_(getTolerance(si, ClpPrimalTolerance))
1645 , ztoldj_(getTolerance(si, ClpDualTolerance))
1646 , maxmin_(si->getObjSense())
1647 , sol_(NULL)
1648 , rowduals_(NULL)
1649 , acts_(NULL)
1650 , rcosts_(NULL)
1651 , colstat_(NULL)
1652 , rowstat_(NULL)
1653 , handler_(NULL)
1654 , defaultHandler_(false)
1655 , messages_()
1656
1657 {
1658 bulk0_ = static_cast< CoinBigIndex >(bulkRatio_ * CoinMax(nelems_in, nelems_)
1659 + ncols_in);
1660 // allow for temporary overflow
1661 hrow_ = new int[bulk0_ + ncols_in];
1662 colels_ = new double[bulk0_ + ncols_in];
1663 si->getDblParam(ClpObjOffset, originalOffset_);
1664 int ncols = si->getNumCols();
1665 int nrows = si->getNumRows();
1666
1667 setMessageHandler(si->messageHandler());
1668
1669 ClpDisjointCopyN(si->getColLower(), ncols, clo_);
1670 ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
1671 //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1672 double offset;
1673 ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
1674 ClpDisjointCopyN(si->getRowLower(), nrows, rlo_);
1675 ClpDisjointCopyN(si->getRowUpper(), nrows, rup_);
1676 int i;
1677 for (i = 0; i < ncols_in; i++)
1678 originalColumn_[i] = i;
1679 for (i = 0; i < nrows_in; i++)
1680 originalRow_[i] = i;
1681 sol_ = NULL;
1682 rowduals_ = NULL;
1683 acts_ = NULL;
1684
1685 rcosts_ = NULL;
1686 colstat_ = NULL;
1687 rowstat_ = NULL;
1688 }
1689
1690 // I am not familiar enough with CoinPackedMatrix to be confident
1691 // that I will implement a row-ordered version of toColumnOrderedGapFree
1692 // properly.
isGapFree(const CoinPackedMatrix & matrix)1693 static bool isGapFree(const CoinPackedMatrix &matrix)
1694 {
1695 const CoinBigIndex *start = matrix.getVectorStarts();
1696 const int *length = matrix.getVectorLengths();
1697 int i = matrix.getSizeVectorLengths() - 1;
1698 // Quick check
1699 if (matrix.getNumElements() == start[i]) {
1700 return true;
1701 } else {
1702 for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1703 if (start[i + 1] - start[i] != length[i])
1704 break;
1705 }
1706 return (!(i >= 0));
1707 }
1708 }
1709 #if PRESOLVE_DEBUG
matrix_bounds_ok(const double * lo,const double * up,int n)1710 static void matrix_bounds_ok(const double *lo, const double *up, int n)
1711 {
1712 int i;
1713 for (i = 0; i < n; i++) {
1714 PRESOLVEASSERT(lo[i] <= up[i]);
1715 PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
1716 PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
1717 }
1718 }
1719 #endif
CoinPresolveMatrix(int ncols0_in,double,ClpSimplex * si,int nrows_in,CoinBigIndex nelems_in,bool doStatus,double nonLinearValue,double bulkRatio)1720 CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1721 double /*maxmin*/,
1722 // end prepost members
1723
1724 ClpSimplex *si,
1725
1726 // rowrep
1727 int nrows_in,
1728 CoinBigIndex nelems_in,
1729 bool doStatus,
1730 double nonLinearValue,
1731 double bulkRatio)
1732 :
1733
1734 CoinPrePostsolveMatrix(si,
1735 ncols0_in, nrows_in, nelems_in, bulkRatio)
1736 , clink_(new presolvehlink[ncols0_in + 1])
1737 , rlink_(new presolvehlink[nrows_in + 1])
1738 ,
1739
1740 dobias_(0.0)
1741 ,
1742
1743 // temporary init
1744 integerType_(new unsigned char[ncols0_in])
1745 , anyInteger_(false)
1746 , tuning_(false)
1747 , startTime_(0.0)
1748 , feasibilityTolerance_(0.0)
1749 , status_(-1)
1750 , pass_(0)
1751 , colsToDo_(new int[ncols0_in])
1752 , numberColsToDo_(0)
1753 , nextColsToDo_(new int[ncols0_in])
1754 , numberNextColsToDo_(0)
1755 , rowsToDo_(new int[nrows_in])
1756 , numberRowsToDo_(0)
1757 , nextRowsToDo_(new int[nrows_in])
1758 , numberNextRowsToDo_(0)
1759 , presolveOptions_(0)
1760 {
1761 const CoinBigIndex bufsize = bulk0_;
1762
1763 nrows_ = si->getNumRows();
1764
1765 // Set up change bits etc
1766 rowChanged_ = new unsigned char[nrows_];
1767 memset(rowChanged_, 0, nrows_);
1768 colChanged_ = new unsigned char[ncols_];
1769 memset(colChanged_, 0, ncols_);
1770 CoinPackedMatrix *m = si->matrix();
1771
1772 // The coefficient matrix is a big hunk of stuff.
1773 // Do the copy here to try to avoid running out of memory.
1774
1775 const CoinBigIndex *start = m->getVectorStarts();
1776 const int *row = m->getIndices();
1777 const double *element = m->getElements();
1778 int icol, nel = 0;
1779 mcstrt_[0] = 0;
1780 ClpDisjointCopyN(m->getVectorLengths(), ncols_, hincol_);
1781 if (si->getObjSense() < 0.0) {
1782 for (int i = 0; i < ncols_; i++)
1783 cost_[i] = -cost_[i];
1784 maxmin_ = 1.0;
1785 }
1786 for (icol = 0; icol < ncols_; icol++) {
1787 CoinBigIndex j;
1788 for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
1789 hrow_[nel] = row[j];
1790 if (fabs(element[j]) > ZTOLDP)
1791 colels_[nel++] = element[j];
1792 }
1793 mcstrt_[icol + 1] = nel;
1794 hincol_[icol] = static_cast< int >(nel - mcstrt_[icol]);
1795 }
1796
1797 // same thing for row rep
1798 CoinPackedMatrix *mRow = new CoinPackedMatrix();
1799 mRow->setExtraGap(0.0);
1800 mRow->setExtraMajor(0.0);
1801 mRow->reverseOrderedCopyOf(*m);
1802 //mRow->removeGaps();
1803 //mRow->setExtraGap(0.0);
1804
1805 // Now get rid of matrix
1806 si->createEmptyMatrix();
1807
1808 double *el = mRow->getMutableElements();
1809 int *ind = mRow->getMutableIndices();
1810 CoinBigIndex *strt = mRow->getMutableVectorStarts();
1811 int *len = mRow->getMutableVectorLengths();
1812 // Do carefully to save memory
1813 rowels_ = new double[bulk0_];
1814 ClpDisjointCopyN(el, nelems_, rowels_);
1815 mRow->nullElementArray();
1816 delete[] el;
1817 hcol_ = new int[bulk0_];
1818 ClpDisjointCopyN(ind, nelems_, hcol_);
1819 mRow->nullIndexArray();
1820 delete[] ind;
1821 mrstrt_ = new CoinBigIndex[nrows_in + 1];
1822 ClpDisjointCopyN(strt, nrows_, mrstrt_);
1823 mRow->nullStartArray();
1824 mrstrt_[nrows_] = nelems_;
1825 delete[] strt;
1826 hinrow_ = new int[nrows_in + 1];
1827 ClpDisjointCopyN(len, nrows_, hinrow_);
1828 if (nelems_ > nel) {
1829 nelems_ = nel;
1830 // Clean any small elements
1831 int irow;
1832 nel = 0;
1833 CoinBigIndex start = 0;
1834 for (irow = 0; irow < nrows_; irow++) {
1835 CoinBigIndex j;
1836 for (j = start; j < start + hinrow_[irow]; j++) {
1837 hcol_[nel] = hcol_[j];
1838 if (fabs(rowels_[j]) > ZTOLDP)
1839 rowels_[nel++] = rowels_[j];
1840 }
1841 start = mrstrt_[irow + 1];
1842 mrstrt_[irow + 1] = nel;
1843 hinrow_[irow] = static_cast< int >(nel - mrstrt_[irow]);
1844 }
1845 }
1846
1847 delete mRow;
1848 if (si->integerInformation()) {
1849 CoinMemcpyN(reinterpret_cast< unsigned char * >(si->integerInformation()), ncols_, integerType_);
1850 } else {
1851 ClpFillN< unsigned char >(integerType_, ncols_, static_cast< unsigned char >(0));
1852 }
1853
1854 #ifndef SLIM_CLP
1855 #ifndef NO_RTTI
1856 ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(si->objectiveAsObject()));
1857 #else
1858 ClpQuadraticObjective *quadraticObj = NULL;
1859 if (si->objectiveAsObject()->type() == 2)
1860 quadraticObj = (static_cast< ClpQuadraticObjective * >(si->objectiveAsObject()));
1861 #endif
1862 #endif
1863 // Set up prohibited bits if needed
1864 if (nonLinearValue) {
1865 anyProhibited_ = true;
1866 for (icol = 0; icol < ncols_; icol++) {
1867 CoinBigIndex j;
1868 bool nonLinearColumn = false;
1869 if (cost_[icol] == nonLinearValue)
1870 nonLinearColumn = true;
1871 for (j = mcstrt_[icol]; j < mcstrt_[icol + 1]; j++) {
1872 if (colels_[j] == nonLinearValue) {
1873 nonLinearColumn = true;
1874 setRowProhibited(hrow_[j]);
1875 }
1876 }
1877 if (nonLinearColumn)
1878 setColProhibited(icol);
1879 }
1880 #ifndef SLIM_CLP
1881 } else if (quadraticObj) {
1882 CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
1883 //const int * columnQuadratic = quadratic->getIndices();
1884 //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1885 const int *columnQuadraticLength = quadratic->getVectorLengths();
1886 //double * quadraticElement = quadratic->getMutableElements();
1887 int numberColumns = quadratic->getNumCols();
1888 anyProhibited_ = true;
1889 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1890 if (columnQuadraticLength[iColumn]) {
1891 setColProhibited(iColumn);
1892 //printf("%d prohib\n",iColumn);
1893 }
1894 }
1895 #endif
1896 } else {
1897 anyProhibited_ = false;
1898 }
1899
1900 if (doStatus) {
1901 // allow for status and solution
1902 sol_ = new double[ncols_];
1903 CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);
1904 ;
1905 acts_ = new double[nrows_];
1906 CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
1907 if (!si->statusArray())
1908 si->createStatus();
1909 colstat_ = new unsigned char[nrows_ + ncols_];
1910 CoinMemcpyN(si->statusArray(), (nrows_ + ncols_), colstat_);
1911 rowstat_ = colstat_ + ncols_;
1912 }
1913
1914 // the original model's fields are now unneeded - free them
1915
1916 si->resize(0, 0);
1917
1918 #if PRESOLVE_DEBUG
1919 matrix_bounds_ok(rlo_, rup_, nrows_);
1920 matrix_bounds_ok(clo_, cup_, ncols_);
1921 #endif
1922
1923 #if 0
1924 for (i = 0; i < nrows; ++i)
1925 printf("NR: %6d\n", hinrow[i]);
1926 for (int i = 0; i < ncols; ++i)
1927 printf("NC: %6d\n", hincol[i]);
1928 #endif
1929
1930 presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
1931 presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
1932
1933 // this allows last col/row to expand up to bufsize-1 (22);
1934 // this must come after the calls to presolve_prefix
1935 mcstrt_[ncols_] = bufsize - 1;
1936 mrstrt_[nrows_] = bufsize - 1;
1937 // Allocate useful arrays
1938 initializeStuff();
1939
1940 #if PRESOLVE_CONSISTENCY
1941 //consistent(false);
1942 presolve_consistent(this, false);
1943 #endif
1944 }
1945
1946 // avoid compiler warnings
1947 #if PRESOLVE_SUMMARY > 0
update_model(ClpSimplex * si,int nrows0,int ncols0,CoinBigIndex nelems0)1948 void CoinPresolveMatrix::update_model(ClpSimplex *si,
1949 int nrows0, int ncols0,
1950 CoinBigIndex nelems0)
1951 #else
1952 void CoinPresolveMatrix::update_model(ClpSimplex *si,
1953 int /*nrows0*/,
1954 int /*ncols0*/,
1955 CoinBigIndex /*nelems0*/)
1956 #endif
1957 {
1958 if (si->getObjSense() < 0.0) {
1959 for (int i = 0; i < ncols_; i++)
1960 cost_[i] = -cost_[i];
1961 dobias_ = -dobias_;
1962 }
1963 si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1964 clo_, cup_, cost_, rlo_, rup_);
1965 //delete [] si->integerInformation();
1966 int numberIntegers = 0;
1967 for (int i = 0; i < ncols_; i++) {
1968 if (integerType_[i])
1969 numberIntegers++;
1970 }
1971 if (numberIntegers)
1972 si->copyInIntegerInformation(reinterpret_cast< const char * >(integerType_));
1973 else
1974 si->copyInIntegerInformation(NULL);
1975
1976 #if PRESOLVE_SUMMARY
1977 printf("NEW NCOL/NROW/NELS: %d(-%d) %d(-%d) %d(-%d)\n",
1978 ncols_, ncols0 - ncols_,
1979 nrows_, nrows0 - nrows_,
1980 si->getNumElements(), nelems0 - si->getNumElements());
1981 #endif
1982 si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
1983 if (si->getObjSense() < 0.0) {
1984 // put back
1985 for (int i = 0; i < ncols_; i++)
1986 cost_[i] = -cost_[i];
1987 dobias_ = -dobias_;
1988 maxmin_ = -1.0;
1989 }
1990 }
1991
1992 //////////////// POSTSOLVE
1993
CoinPostsolveMatrix(ClpSimplex * si,int ncols0_in,int nrows0_in,CoinBigIndex nelems0,double maxmin,double * sol_in,double * acts_in,unsigned char * colstat_in,unsigned char * rowstat_in)1994 CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex *si,
1995 int ncols0_in,
1996 int nrows0_in,
1997 CoinBigIndex nelems0,
1998
1999 double maxmin,
2000 // end prepost members
2001
2002 double *sol_in,
2003 double *acts_in,
2004
2005 unsigned char *colstat_in,
2006 unsigned char *rowstat_in)
2007 : CoinPrePostsolveMatrix(si,
2008 ncols0_in, nrows0_in, nelems0, 2.0)
2009 ,
2010
2011 free_list_(0)
2012 ,
2013 // link, free_list, maxlink
2014 maxlink_(bulk0_)
2015 , link_(new CoinBigIndex[/*maxlink*/ bulk0_])
2016 ,
2017
2018 cdone_(new char[ncols0_])
2019 , rdone_(new char[nrows0_in])
2020
2021 {
2022 bulk0_ = maxlink_;
2023 nrows_ = si->getNumRows();
2024 ncols_ = si->getNumCols();
2025
2026 sol_ = sol_in;
2027 rowduals_ = NULL;
2028 acts_ = acts_in;
2029
2030 rcosts_ = NULL;
2031 colstat_ = colstat_in;
2032 rowstat_ = rowstat_in;
2033
2034 // this is the *reduced* model, which is probably smaller
2035 int ncols1 = ncols_;
2036 int nrows1 = nrows_;
2037
2038 const CoinPackedMatrix *m = si->matrix();
2039
2040 const CoinBigIndex nelemsr = m->getNumElements();
2041 if (m->getNumElements() && !isGapFree(*m)) {
2042 // Odd - gaps
2043 CoinPackedMatrix mm(*m);
2044 mm.removeGaps();
2045 mm.setExtraGap(0.0);
2046
2047 ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
2048 CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
2049 mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
2050 ClpDisjointCopyN(mm.getVectorLengths(), ncols1, hincol_);
2051 ClpDisjointCopyN(mm.getIndices(), nelemsr, hrow_);
2052 ClpDisjointCopyN(mm.getElements(), nelemsr, colels_);
2053 } else {
2054 // No gaps
2055
2056 ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
2057 CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
2058 mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
2059 ClpDisjointCopyN(m->getVectorLengths(), ncols1, hincol_);
2060 ClpDisjointCopyN(m->getIndices(), nelemsr, hrow_);
2061 ClpDisjointCopyN(m->getElements(), nelemsr, colels_);
2062 }
2063
2064 #if 0 && PRESOLVE_DEBUG
2065 presolve_check_costs(model, &colcopy);
2066 #endif
2067
2068 // This determines the size of the data structure that contains
2069 // the matrix being postsolved. Links are taken from the free_list
2070 // to recreate matrix entries that were presolved away,
2071 // and links are added to the free_list when entries created during
2072 // presolve are discarded. There is never a need to gc this list.
2073 // Naturally, it should contain
2074 // exactly nelems0 entries "in use" when postsolving is done,
2075 // but I don't know whether the matrix could temporarily get
2076 // larger during postsolving. Substitution into more than two
2077 // rows could do that, in principle. I am being very conservative
2078 // here by reserving much more than the amount of space I probably need.
2079 // If this guess is wrong, check_free_list may be called.
2080 // int bufsize = 2*nelems0;
2081
2082 memset(cdone_, -1, ncols0_);
2083 memset(rdone_, -1, nrows0_);
2084
2085 rowduals_ = new double[nrows0_];
2086 ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
2087
2088 rcosts_ = new double[ncols0_];
2089 ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
2090 if (maxmin < 0.0) {
2091 // change so will look as if minimize
2092 int i;
2093 for (i = 0; i < nrows1; i++)
2094 rowduals_[i] = -rowduals_[i];
2095 for (i = 0; i < ncols1; i++) {
2096 rcosts_[i] = -rcosts_[i];
2097 }
2098 }
2099
2100 //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
2101 //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
2102
2103 ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
2104 si->setDblParam(ClpObjOffset, originalOffset_);
2105 // Test below only needed for QP ..... but .....
2106 // To switch off define COIN_SLOW_PRESOLVE=0
2107 #ifndef COIN_SLOW_PRESOLVE
2108 #define COIN_SLOW_PRESOLVE 1
2109 #endif
2110 for (int j = 0; j < ncols1; j++) {
2111 #if COIN_SLOW_PRESOLVE
2112 if (hincol_[j]) {
2113 #endif
2114 CoinBigIndex kcs = mcstrt_[j];
2115 CoinBigIndex kce = kcs + hincol_[j];
2116 for (CoinBigIndex k = kcs; k < kce; ++k) {
2117 link_[k] = k + 1;
2118 }
2119 link_[kce - 1] = NO_LINK;
2120 #if COIN_SLOW_PRESOLVE
2121 }
2122 #endif
2123 }
2124 {
2125 CoinBigIndex ml = maxlink_;
2126 for (CoinBigIndex k = nelemsr; k < ml; ++k)
2127 link_[k] = k + 1;
2128 if (ml)
2129 link_[ml - 1] = NO_LINK;
2130 }
2131 free_list_ = nelemsr;
2132 #if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
2133 /*
2134 These are used to track the action of postsolve transforms during debugging.
2135 */
2136 CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED);
2137 CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1);
2138 CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED);
2139 CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1);
2140 #endif
2141 }
2142 /* This is main part of Presolve */
2143 ClpSimplex *
gutsOfPresolvedModel(ClpSimplex * originalModel,double feasibilityTolerance,bool keepIntegers,int numberPasses,bool dropNames,bool doRowObjective,const char * prohibitedRows,const char * prohibitedColumns)2144 ClpPresolve::gutsOfPresolvedModel(ClpSimplex *originalModel,
2145 double feasibilityTolerance,
2146 bool keepIntegers,
2147 int numberPasses,
2148 bool dropNames,
2149 bool doRowObjective,
2150 const char *prohibitedRows,
2151 const char *prohibitedColumns)
2152 {
2153 ncols_ = originalModel->getNumCols();
2154 nrows_ = originalModel->getNumRows();
2155 nelems_ = originalModel->getNumElements();
2156 numberPasses_ = numberPasses;
2157
2158 double maxmin = originalModel->getObjSense();
2159 originalModel_ = originalModel;
2160 delete[] originalColumn_;
2161 originalColumn_ = new int[ncols_];
2162 delete[] originalRow_;
2163 originalRow_ = new int[nrows_];
2164 // and fill in case returns early
2165 int i;
2166 for (i = 0; i < ncols_; i++)
2167 originalColumn_[i] = i;
2168 for (i = 0; i < nrows_; i++)
2169 originalRow_[i] = i;
2170 delete[] rowObjective_;
2171 if (doRowObjective) {
2172 rowObjective_ = new double[nrows_];
2173 memset(rowObjective_, 0, nrows_ * sizeof(double));
2174 } else {
2175 rowObjective_ = NULL;
2176 }
2177
2178 // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
2179 int result = -1;
2180
2181 // User may have deleted - its their responsibility
2182 presolvedModel_ = NULL;
2183 // Messages
2184 CoinMessages messages = originalModel->coinMessages();
2185 // Only go round 100 times even if integer preprocessing
2186 int totalPasses = 100;
2187 while (result == -1) {
2188
2189 #ifndef CLP_NO_STD
2190 // make new copy
2191 if (saveFile_ == "") {
2192 #endif
2193 delete presolvedModel_;
2194 #ifndef CLP_NO_STD
2195 // So won't get names
2196 int lengthNames = originalModel->lengthNames();
2197 originalModel->setLengthNames(0);
2198 #endif
2199 presolvedModel_ = new ClpSimplex(*originalModel);
2200 #ifndef CLP_NO_STD
2201 originalModel->setLengthNames(lengthNames);
2202 presolvedModel_->dropNames();
2203 } else {
2204 presolvedModel_ = originalModel;
2205 if (dropNames)
2206 presolvedModel_->dropNames();
2207 }
2208 #endif
2209
2210 // drop integer information if wanted
2211 if (!keepIntegers)
2212 presolvedModel_->deleteIntegerInformation();
2213 totalPasses--;
2214
2215 double ratio = 2.0;
2216 if (substitution_ > 3)
2217 ratio = sqrt((substitution_ - 3) + 5.0);
2218 else if (substitution_ == 2)
2219 ratio = 1.5;
2220 CoinPresolveMatrix prob(ncols_,
2221 maxmin,
2222 presolvedModel_,
2223 nrows_, nelems_, true, nonLinearValue_, ratio);
2224 if (prohibitedRows) {
2225 prob.setAnyProhibited();
2226 for (int i = 0; i < nrows_; i++) {
2227 if (prohibitedRows[i])
2228 prob.setRowProhibited(i);
2229 }
2230 }
2231 if (prohibitedColumns) {
2232 prob.setAnyProhibited();
2233 for (int i = 0; i < ncols_; i++) {
2234 if (prohibitedColumns[i])
2235 prob.setColProhibited(i);
2236 }
2237 }
2238 prob.setMaximumSubstitutionLevel(substitution_);
2239 if (doRowObjective)
2240 memset(rowObjective_, 0, nrows_ * sizeof(double));
2241 // See if we want statistics
2242 if ((presolveActions_ & 0x80000000) != 0)
2243 prob.statistics();
2244 if (doTransfer())
2245 transferCosts(&prob);
2246 // make sure row solution correct
2247 {
2248 double *colels = prob.colels_;
2249 int *hrow = prob.hrow_;
2250 CoinBigIndex *mcstrt = prob.mcstrt_;
2251 int *hincol = prob.hincol_;
2252 int ncols = prob.ncols_;
2253
2254 double *csol = prob.sol_;
2255 double *acts = prob.acts_;
2256 int nrows = prob.nrows_;
2257
2258 int colx;
2259
2260 memset(acts, 0, nrows * sizeof(double));
2261
2262 for (colx = 0; colx < ncols; ++colx) {
2263 double solutionValue = csol[colx];
2264 for (CoinBigIndex i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
2265 int row = hrow[i];
2266 double coeff = colels[i];
2267 acts[row] += solutionValue * coeff;
2268 }
2269 }
2270 }
2271
2272 // move across feasibility tolerance
2273 prob.feasibilityTolerance_ = feasibilityTolerance;
2274
2275 // Do presolve
2276 paction_ = presolve(&prob);
2277 // Get rid of useful arrays
2278 prob.deleteStuff();
2279
2280 result = 0;
2281
2282 bool fixInfeasibility = (prob.presolveOptions_ & 16384) != 0;
2283 bool hasSolution = (prob.presolveOptions_ & 32768) != 0;
2284 if (prob.status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
2285 // Looks feasible but double check to see if anything slipped through
2286 int n = prob.ncols_;
2287 double *lo = prob.clo_;
2288 double *up = prob.cup_;
2289 int i;
2290
2291 for (i = 0; i < n; i++) {
2292 if (up[i] < lo[i]) {
2293 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2294 // infeasible
2295 prob.status_ = 1;
2296 } else {
2297 up[i] = lo[i];
2298 }
2299 }
2300 }
2301
2302 n = prob.nrows_;
2303 lo = prob.rlo_;
2304 up = prob.rup_;
2305
2306 for (i = 0; i < n; i++) {
2307 if (up[i] < lo[i]) {
2308 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2309 // infeasible
2310 prob.status_ = 1;
2311 } else {
2312 up[i] = lo[i];
2313 }
2314 }
2315 }
2316 }
2317 if (prob.status_ == 0 && paction_) {
2318 // feasible
2319
2320 prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
2321 // copy status and solution
2322 CoinMemcpyN(prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
2323 CoinMemcpyN(prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
2324 CoinMemcpyN(prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
2325 CoinMemcpyN(prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
2326 if (fixInfeasibility && hasSolution) {
2327 // Looks feasible but double check to see if anything slipped through
2328 int n = prob.ncols_;
2329 double *lo = prob.clo_;
2330 double *up = prob.cup_;
2331 double *rsol = prob.acts_;
2332 //memset(prob.acts_,0,prob.nrows_*sizeof(double));
2333 presolvedModel_->matrix()->times(prob.sol_, rsol);
2334 int i;
2335
2336 for (i = 0; i < n; i++) {
2337 double gap = up[i] - lo[i];
2338 if (rsol[i] < lo[i] - feasibilityTolerance && fabs(rsol[i] - lo[i]) < 1.0e-3) {
2339 lo[i] = rsol[i];
2340 if (gap < 1.0e5)
2341 up[i] = lo[i] + gap;
2342 } else if (rsol[i] > up[i] + feasibilityTolerance && fabs(rsol[i] - up[i]) < 1.0e-3) {
2343 up[i] = rsol[i];
2344 if (gap < 1.0e5)
2345 lo[i] = up[i] - gap;
2346 }
2347 if (up[i] < lo[i]) {
2348 up[i] = lo[i];
2349 }
2350 }
2351 }
2352
2353 int n = prob.nrows_;
2354 double *lo = prob.rlo_;
2355 double *up = prob.rup_;
2356
2357 for (i = 0; i < n; i++) {
2358 if (up[i] < lo[i]) {
2359 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
2360 // infeasible
2361 prob.status_ = 1;
2362 } else {
2363 up[i] = lo[i];
2364 }
2365 }
2366 }
2367 delete[] prob.sol_;
2368 delete[] prob.acts_;
2369 delete[] prob.colstat_;
2370 prob.sol_ = NULL;
2371 prob.acts_ = NULL;
2372 prob.colstat_ = NULL;
2373
2374 int ncolsNow = presolvedModel_->getNumCols();
2375 CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
2376 #ifndef SLIM_CLP
2377 #ifndef NO_RTTI
2378 ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(originalModel->objectiveAsObject()));
2379 #else
2380 ClpQuadraticObjective *quadraticObj = NULL;
2381 if (originalModel->objectiveAsObject()->type() == 2)
2382 quadraticObj = (static_cast< ClpQuadraticObjective * >(originalModel->objectiveAsObject()));
2383 #endif
2384 if (quadraticObj) {
2385 // set up for subset
2386 char *mark = new char[ncols_];
2387 memset(mark, 0, ncols_);
2388 CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
2389 //const int * columnQuadratic = quadratic->getIndices();
2390 //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
2391 const int *columnQuadraticLength = quadratic->getVectorLengths();
2392 //double * quadraticElement = quadratic->getMutableElements();
2393 int numberColumns = quadratic->getNumCols();
2394 ClpQuadraticObjective *newObj = new ClpQuadraticObjective(*quadraticObj,
2395 ncolsNow,
2396 originalColumn_);
2397 // and modify linear and check
2398 double *linear = newObj->linearObjective();
2399 CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
2400 int iColumn;
2401 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2402 if (columnQuadraticLength[iColumn])
2403 mark[iColumn] = 1;
2404 }
2405 // and new
2406 quadratic = newObj->quadraticObjective();
2407 columnQuadraticLength = quadratic->getVectorLengths();
2408 int numberColumns2 = quadratic->getNumCols();
2409 for (iColumn = 0; iColumn < numberColumns2; iColumn++) {
2410 if (columnQuadraticLength[iColumn])
2411 mark[originalColumn_[iColumn]] = 0;
2412 }
2413 presolvedModel_->setObjective(newObj);
2414 delete newObj;
2415 // final check
2416 for (iColumn = 0; iColumn < numberColumns; iColumn++)
2417 if (mark[iColumn])
2418 printf("Quadratic column %d modified - may be okay\n", iColumn);
2419 delete[] mark;
2420 }
2421 #endif
2422 delete[] prob.originalColumn_;
2423 prob.originalColumn_ = NULL;
2424 int nrowsNow = presolvedModel_->getNumRows();
2425 CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
2426 delete[] prob.originalRow_;
2427 prob.originalRow_ = NULL;
2428 #ifndef CLP_NO_STD
2429 if (!dropNames && originalModel->lengthNames()) {
2430 // Redo names
2431 int iRow;
2432 std::vector< std::string > rowNames;
2433 rowNames.reserve(nrowsNow);
2434 for (iRow = 0; iRow < nrowsNow; iRow++) {
2435 int kRow = originalRow_[iRow];
2436 rowNames.push_back(originalModel->rowName(kRow));
2437 }
2438
2439 int iColumn;
2440 std::vector< std::string > columnNames;
2441 columnNames.reserve(ncolsNow);
2442 for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
2443 int kColumn = originalColumn_[iColumn];
2444 columnNames.push_back(originalModel->columnName(kColumn));
2445 }
2446 presolvedModel_->copyNames(rowNames, columnNames);
2447 } else {
2448 presolvedModel_->setLengthNames(0);
2449 }
2450 #endif
2451 if (rowObjective_) {
2452 int iRow;
2453 #ifndef NDEBUG
2454 int k = -1;
2455 #endif
2456 int nObj = 0;
2457 for (iRow = 0; iRow < nrowsNow; iRow++) {
2458 int kRow = originalRow_[iRow];
2459 #ifndef NDEBUG
2460 assert(kRow > k);
2461 k = kRow;
2462 #endif
2463 rowObjective_[iRow] = rowObjective_[kRow];
2464 if (rowObjective_[iRow])
2465 nObj++;
2466 }
2467 if (nObj) {
2468 printf("%d costed slacks\n", nObj);
2469 presolvedModel_->setRowObjective(rowObjective_);
2470 }
2471 }
2472 /* now clean up integer variables. This can modify original
2473 Don't do if dupcol added columns together */
2474 int i;
2475 const char *information = presolvedModel_->integerInformation();
2476 if ((prob.presolveOptions_ & 0x80000000) == 0 && information) {
2477 int numberChanges = 0;
2478 double *lower0 = originalModel_->columnLower();
2479 double *upper0 = originalModel_->columnUpper();
2480 double *lower = presolvedModel_->columnLower();
2481 double *upper = presolvedModel_->columnUpper();
2482 for (i = 0; i < ncolsNow; i++) {
2483 if (!information[i])
2484 continue;
2485 int iOriginal = originalColumn_[i];
2486 double lowerValue0 = lower0[iOriginal];
2487 double upperValue0 = upper0[iOriginal];
2488 double lowerValue = ceil(lower[i] - 1.0e-5);
2489 double upperValue = floor(upper[i] + 1.0e-5);
2490 lower[i] = lowerValue;
2491 upper[i] = upperValue;
2492 if (lowerValue > upperValue) {
2493 numberChanges++;
2494 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
2495 messages)
2496 << iOriginal
2497 << lowerValue
2498 << upperValue
2499 << CoinMessageEol;
2500 result = 1;
2501 } else {
2502 if (lowerValue > lowerValue0 + 1.0e-8) {
2503 lower0[iOriginal] = lowerValue;
2504 numberChanges++;
2505 }
2506 if (upperValue < upperValue0 - 1.0e-8) {
2507 upper0[iOriginal] = upperValue;
2508 numberChanges++;
2509 }
2510 }
2511 }
2512 if (numberChanges) {
2513 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
2514 messages)
2515 << numberChanges
2516 << CoinMessageEol;
2517 if (!result && totalPasses > 0) {
2518 result = -1; // round again
2519 const CoinPresolveAction *paction = paction_;
2520 while (paction) {
2521 const CoinPresolveAction *next = paction->next;
2522 delete paction;
2523 paction = next;
2524 }
2525 paction_ = NULL;
2526 }
2527 }
2528 }
2529 } else if (prob.status_) {
2530 // infeasible or unbounded
2531 result = 1;
2532 // Put status in nelems_!
2533 nelems_ = -prob.status_;
2534 originalModel->setProblemStatus(prob.status_);
2535 } else {
2536 // no changes - model needs restoring after Lou's changes
2537 #ifndef CLP_NO_STD
2538 if (saveFile_ == "") {
2539 #endif
2540 delete presolvedModel_;
2541 presolvedModel_ = new ClpSimplex(*originalModel);
2542 // but we need to remove gaps
2543 ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(presolvedModel_->clpMatrix());
2544 if (clpMatrix) {
2545 clpMatrix->getPackedMatrix()->removeGaps();
2546 }
2547 #ifndef CLP_NO_STD
2548 } else {
2549 presolvedModel_ = originalModel;
2550 }
2551 presolvedModel_->dropNames();
2552 #endif
2553
2554 // drop integer information if wanted
2555 if (!keepIntegers)
2556 presolvedModel_->deleteIntegerInformation();
2557 result = 2;
2558 }
2559 }
2560 if (result == 0 || result == 2) {
2561 int nrowsAfter = presolvedModel_->getNumRows();
2562 int ncolsAfter = presolvedModel_->getNumCols();
2563 CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
2564 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
2565 messages)
2566 << nrowsAfter << -(nrows_ - nrowsAfter)
2567 << ncolsAfter << -(ncols_ - ncolsAfter)
2568 << nelsAfter << -(nelems_ - nelsAfter)
2569 << CoinMessageEol;
2570 } else {
2571 destroyPresolve();
2572 if (presolvedModel_ != originalModel_)
2573 delete presolvedModel_;
2574 presolvedModel_ = NULL;
2575 }
2576 return presolvedModel_;
2577 }
2578
2579 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2580 */
2581