1/* $Id: ClpPresolve.cpp 1802 2011-10-11 17:08:33Z forrest $ */
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
16#include "CoinPackedMatrix.hpp"
17#include "ClpPackedMatrix.hpp"
18#include "ClpSimplex.hpp"
19#include "ClpSimplexOther.hpp"
20#ifndef SLIM_CLP
21#include "ClpQuadraticObjective.hpp"
22#endif
23
24#include "ClpPresolve.hpp"
25#include "CoinPresolveMatrix.hpp"
26
27#include "CoinPresolveEmpty.hpp"
28#include "CoinPresolveFixed.hpp"
29#include "CoinPresolvePsdebug.hpp"
30#include "CoinPresolveSingleton.hpp"
31#include "CoinPresolveDoubleton.hpp"
32#include "CoinPresolveTripleton.hpp"
33#include "CoinPresolveZeros.hpp"
34#include "CoinPresolveSubst.hpp"
35#include "CoinPresolveForcing.hpp"
36#include "CoinPresolveDual.hpp"
37#include "CoinPresolveTighten.hpp"
38#include "CoinPresolveUseless.hpp"
39#include "CoinPresolveDupcol.hpp"
40#include "CoinPresolveImpliedFree.hpp"
41#include "CoinPresolveIsolated.hpp"
42#include "CoinMessage.hpp"
43
44
45
46ClpPresolve::ClpPresolve() :
47 originalModel_(NULL),
48 presolvedModel_(NULL),
49 nonLinearValue_(0.0),
50 originalColumn_(NULL),
51 originalRow_(NULL),
52 rowObjective_(NULL),
53 paction_(0),
54 ncols_(0),
55 nrows_(0),
56 nelems_(0),
57 numberPasses_(5),
58 substitution_(3),
59#ifndef CLP_NO_STD
60 saveFile_(""),
61#endif
62 presolveActions_(0)
63{
64}
65
66ClpPresolve::~ClpPresolve()
67{
68 destroyPresolve();
69}
70// Gets rid of presolve actions (e.g.when infeasible)
71void
72ClpPresolve::destroyPresolve()
73{
74 const CoinPresolveAction *paction = paction_;
75 while (paction) {
76 const CoinPresolveAction *next = paction->next;
77 delete paction;
78 paction = next;
79 }
80 delete [] originalColumn_;
81 delete [] originalRow_;
82 paction_ = NULL;
83 originalColumn_ = NULL;
84 originalRow_ = NULL;
85 delete [] rowObjective_;
86 rowObjective_ = NULL;
87}
88
89/* This version of presolve returns a pointer to a new presolved
90 model. NULL if infeasible
91*/
92ClpSimplex *
93ClpPresolve::presolvedModel(ClpSimplex & si,
94 double feasibilityTolerance,
95 bool keepIntegers,
96 int numberPasses,
97 bool dropNames,
98 bool doRowObjective)
99{
100 // Check matrix
101 int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
102 if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
103 1.0e20,checkType))
104 return NULL;
105 else
106 return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
107 doRowObjective);
108}
109#ifndef CLP_NO_STD
110/* This version of presolve updates
111 model and saves original data to file. Returns non-zero if infeasible
112*/
113int
114ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
115 double feasibilityTolerance,
116 bool keepIntegers,
117 int numberPasses,
118 bool dropNames,
119 bool doRowObjective)
120{
121 // Check matrix
122 if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
123 1.0e20))
124 return 2;
125 saveFile_ = fileName;
126 si.saveModel(saveFile_.c_str());
127 ClpSimplex * model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
128 doRowObjective);
129 if (model == &si) {
130 return 0;
131 } else {
132 si.restoreModel(saveFile_.c_str());
133 remove(saveFile_.c_str());
134 return 1;
135 }
136}
137#endif
138// Return pointer to presolved model
139ClpSimplex *
140ClpPresolve::model() const
141{
142 return presolvedModel_;
143}
144// Return pointer to original model
145ClpSimplex *
146ClpPresolve::originalModel() const
147{
148 return originalModel_;
149}
150// Return presolve status (0,1,2)
151int
152ClpPresolve::presolveStatus() const
153{
154 if (nelems_>=0) {
155 // feasible (or not done yet)
156 return 0;
157 } else {
158 int presolveStatus = - nelems_;
159 // If both infeasible and unbounded - say infeasible
160 if (presolveStatus>2)
161 presolveStatus = 1;
162 return presolveStatus;
163 }
164}
165void
166ClpPresolve::postsolve(bool updateStatus)
167{
168 // Return at once if no presolved model
169 if (!presolvedModel_)
170 return;
171 // Messages
172 CoinMessages messages = originalModel_->coinMessages();
173 if (!presolvedModel_->isProvenOptimal()) {
174 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
175 messages)
176 << CoinMessageEol;
177 }
178
179 // this is the size of the original problem
180 const int ncols0 = ncols_;
181 const int nrows0 = nrows_;
182 const CoinBigIndex nelems0 = nelems_;
183
184 // this is the reduced problem
185 int ncols = presolvedModel_->getNumCols();
186 int nrows = presolvedModel_->getNumRows();
187
188 double * acts = NULL;
189 double * sol = NULL;
190 unsigned char * rowstat = NULL;
191 unsigned char * colstat = NULL;
192#ifndef CLP_NO_STD
193 if (saveFile_ == "") {
194#endif
195 // reality check
196 assert(ncols0 == originalModel_->getNumCols());
197 assert(nrows0 == originalModel_->getNumRows());
198 acts = originalModel_->primalRowSolution();
199 sol = originalModel_->primalColumnSolution();
200 if (updateStatus) {
201 // postsolve does not know about fixed
202 int i;
203 for (i = 0; i < nrows + ncols; i++) {
204 if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
205 presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
206 }
207 unsigned char *status = originalModel_->statusArray();
208 if (!status) {
209 originalModel_->createStatus();
210 status = originalModel_->statusArray();
211 }
212 rowstat = status + ncols0;
213 colstat = status;
214 CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
215 CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
216 }
217#ifndef CLP_NO_STD
218 } else {
219 // from file
220 acts = new double[nrows0];
221 sol = new double[ncols0];
222 CoinZeroN(acts, nrows0);
223 CoinZeroN(sol, ncols0);
224 if (updateStatus) {
225 unsigned char *status = new unsigned char [nrows0+ncols0];
226 rowstat = status + ncols0;
227 colstat = status;
228 CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
229 CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
230 }
231 }
232#endif
233
234 // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
235 // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
236 // cause duplicate free. In case where saveFile == "", as best I can see
237 // arrays are owned by originalModel_. fix is to
238 // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
239 CoinPostsolveMatrix prob(presolvedModel_,
240 ncols0,
241 nrows0,
242 nelems0,
243 presolvedModel_->getObjSense(),
244 // end prepost
245
246 sol, acts,
247 colstat, rowstat);
248
249 postsolve(prob);
250
251#ifndef CLP_NO_STD
252 if (saveFile_ != "") {
253 // From file
254 assert (originalModel_ == presolvedModel_);
255 originalModel_->restoreModel(saveFile_.c_str());
256 remove(saveFile_.c_str());
257 CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
258 // delete [] acts;
259 CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
260 // delete [] sol;
261 if (updateStatus) {
262 CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
263 // delete [] colstat;
264 }
265 } else {
266#endif
267 prob.sol_ = 0 ;
268 prob.acts_ = 0 ;
269 prob.colstat_ = 0 ;
270#ifndef CLP_NO_STD
271 }
272#endif
273 // put back duals
274 CoinMemcpyN(prob.rowduals_, nrows_, originalModel_->dualRowSolution());
275 double maxmin = originalModel_->getObjSense();
276 if (maxmin < 0.0) {
277 // swap signs
278 int i;
279 double * pi = originalModel_->dualRowSolution();
280 for (i = 0; i < nrows_; i++)
281 pi[i] = -pi[i];
282 }
283 // Now check solution
284 double offset;
285 CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
286 originalModel_->primalColumnSolution(), offset, true),
287 ncols_, originalModel_->dualColumnSolution());
288 originalModel_->clpMatrix()->transposeTimes(-1.0,
289 originalModel_->dualRowSolution(),
290 originalModel_->dualColumnSolution());
291 memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
292 originalModel_->clpMatrix()->times(1.0, originalModel_->primalColumnSolution(),
293 originalModel_->primalRowSolution());
294 originalModel_->checkSolutionInternal();
295 if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
296 // See if we can fix easily
297 static_cast<ClpSimplexOther *> (originalModel_)->cleanupAfterPostsolve();
298 }
299 // Messages
300 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
301 messages)
302 << originalModel_->objectiveValue()
303 << originalModel_->sumDualInfeasibilities()
304 << originalModel_->numberDualInfeasibilities()
305 << originalModel_->sumPrimalInfeasibilities()
306 << originalModel_->numberPrimalInfeasibilities()
307 << CoinMessageEol;
308
309 //originalModel_->objectiveValue_=objectiveValue_;
310 originalModel_->setNumberIterations(presolvedModel_->numberIterations());
311 if (!presolvedModel_->status()) {
312 if (!originalModel_->numberDualInfeasibilities() &&
313 !originalModel_->numberPrimalInfeasibilities()) {
314 originalModel_->setProblemStatus( 0);
315 } else {
316 originalModel_->setProblemStatus( -1);
317 // Say not optimal after presolve
318 originalModel_->setSecondaryStatus(7);
319 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
320 messages)
321 << CoinMessageEol;
322 }
323 } else {
324 originalModel_->setProblemStatus( presolvedModel_->status());
325 // but not if close to feasible
326 if( originalModel_->sumPrimalInfeasibilities()<1.0e-1) {
327 originalModel_->setProblemStatus( -1);
328 // Say not optimal after presolve
329 originalModel_->setSecondaryStatus(7);
330 }
331 }
332#ifndef CLP_NO_STD
333 if (saveFile_ != "")
334 presolvedModel_ = NULL;
335#endif
336}
337
338// return pointer to original columns
339const int *
340ClpPresolve::originalColumns() const
341{
342 return originalColumn_;
343}
344// return pointer to original rows
345const int *
346ClpPresolve::originalRows() const
347{
348 return originalRow_;
349}
350// Set pointer to original model
351void
352ClpPresolve::setOriginalModel(ClpSimplex * model)
353{
354 originalModel_ = model;
355}
356#if 0
357// A lazy way to restrict which transformations are applied
358// during debugging.
359static int ATOI(const char *name)
360{
361 return true;
362#if PRESOLVE_DEBUG || PRESOLVE_SUMMARY
363 if (getenv(name)) {
364 int val = atoi(getenv(name));
365 printf("%s = %d\n", name, val);
366 return (val);
367 } else {
368 if (strcmp(name, "off"))
369 return (true);
370 else
371 return (false);
372 }
373#else
374 return (true);
375#endif
376}
377#endif
378//#define PRESOLVE_DEBUG 1
379#if PRESOLVE_DEBUG
380void check_sol(CoinPresolveMatrix *prob, double tol)
381{
382 double *colels = prob->colels_;
383 int *hrow = prob->hrow_;
384 int *mcstrt = prob->mcstrt_;
385 int *hincol = prob->hincol_;
386 int *hinrow = prob->hinrow_;
387 int ncols = prob->ncols_;
388
389
390 double * csol = prob->sol_;
391 double * acts = prob->acts_;
392 double * clo = prob->clo_;
393 double * cup = prob->cup_;
394 int nrows = prob->nrows_;
395 double * rlo = prob->rlo_;
396 double * rup = prob->rup_;
397
398 int colx;
399
400 double * rsol = new double[nrows];
401 memset(rsol, 0, nrows * sizeof(double));
402
403 for (colx = 0; colx < ncols; ++colx) {
404 if (1) {
405 CoinBigIndex k = mcstrt[colx];
406 int nx = hincol[colx];
407 double solutionValue = csol[colx];
408 for (int i = 0; i < nx; ++i) {
409 int row = hrow[k];
410 double coeff = colels[k];
411 k++;
412 rsol[row] += solutionValue * coeff;
413 }
414 if (csol[colx] < clo[colx] - tol) {
415 printf("low CSOL: %d - %g %g %g\n",
416 colx, clo[colx], csol[colx], cup[colx]);
417 } else if (csol[colx] > cup[colx] + tol) {
418 printf("high CSOL: %d - %g %g %g\n",
419 colx, clo[colx], csol[colx], cup[colx]);
420 }
421 }
422 }
423 int rowx;
424 for (rowx = 0; rowx < nrows; ++rowx) {
425 if (hinrow[rowx]) {
426 if (fabs(rsol[rowx] - acts[rowx]) > tol)
427 printf("inacc RSOL: %d - %g %g (acts_ %g) %g\n",
428 rowx, rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
429 if (rsol[rowx] < rlo[rowx] - tol) {
430 printf("low RSOL: %d - %g %g %g\n",
431 rowx, rlo[rowx], rsol[rowx], rup[rowx]);
432 } else if (rsol[rowx] > rup[rowx] + tol ) {
433 printf("high RSOL: %d - %g %g %g\n",
434 rowx, rlo[rowx], rsol[rowx], rup[rowx]);
435 }
436 }
437 }
438 delete [] rsol;
439}
440#endif
441//#define COIN_PRESOLVE_BUG
442#ifdef COIN_PRESOLVE_BUG
443static int counter=1000000;
444static int startEmptyRows=0;
445static int startEmptyColumns=0;
446static bool break2(CoinPresolveMatrix *prob)
447{
448 int droppedRows = prob->countEmptyRows() - startEmptyRows ;
449 int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
450 startEmptyRows=prob->countEmptyRows();
451 startEmptyColumns=prob->countEmptyCols();
452 printf("Dropped %d rows and %d columns - current empty %d, %d\n",droppedRows,
453 droppedColumns,startEmptyRows,startEmptyColumns);
454 counter--;
455 if (!counter) {
456 printf("skipping next and all\n");
457 }
458 return (counter<=0);
459}
460#define possibleBreak if (break2(prob)) break
461#define possibleSkip if (!break2(prob))
462#else
463#define possibleBreak
464#define possibleSkip
465#endif
466// This is the presolve loop.
467// It is a separate virtual function so that it can be easily
468// customized by subclassing CoinPresolve.
469const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
470{
471 // Messages
472 CoinMessages messages = CoinMessage(prob->messages().language());
473 paction_ = 0;
474#ifndef PRESOLVE_DETAIL
475 if (prob->tuning_) {
476#endif
477 int numberEmptyRows=0;
478 for ( int i=0;i<prob->nrows_;i++) {
479 if (!prob->hinrow_[i]) {
480 PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n",i));
481 //printf("pre_empty row %d\n",i);
482 numberEmptyRows++;
483 }
484 }
485 int numberEmptyCols=0;
486 for ( int i=0;i<prob->ncols_;i++) {
487 if (!prob->hincol_[i]) {
488 PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n",i));
489 //printf("pre_empty col %d\n",i);
490 numberEmptyCols++;
491 }
492 }
493 printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
494 numberEmptyRows,numberEmptyCols);
495#ifndef PRESOLVE_DETAIL
496 }
497#endif
498 prob->status_ = 0; // say feasible
499 paction_ = make_fixed(prob, paction_);
500 // if integers then switch off dual stuff
501 // later just do individually
502 bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
503 // but allow in some cases
504 if ((presolveActions_ & 512) != 0)
505 doDualStuff = true;
506 if (prob->anyProhibited())
507 doDualStuff = false;
508 if (!doDual())
509 doDualStuff = false;
510#if PRESOLVE_CONSISTENCY
511// presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
512 presolve_links_ok(prob, false, true) ;
513#endif
514
515 if (!prob->status_) {
516 bool slackSingleton = doSingletonColumn();
517 slackSingleton = true;
518 const bool slackd = doSingleton();
519 const bool doubleton = doDoubleton();
520 const bool tripleton = doTripleton();
521 //#define NO_FORCING
522#ifndef NO_FORCING
523 const bool forcing = doForcing();
524#endif
525 const bool ifree = doImpliedFree();
526 const bool zerocost = doTighten();
527 const bool dupcol = doDupcol();
528 const bool duprow = doDuprow();
529 const bool dual = doDualStuff;
530
531 // some things are expensive so just do once (normally)
532
533 int i;
534 // say look at all
535 if (!prob->anyProhibited()) {
536 for (i = 0; i < nrows_; i++)
537 prob->rowsToDo_[i] = i;
538 prob->numberRowsToDo_ = nrows_;
539 for (i = 0; i < ncols_; i++)
540 prob->colsToDo_[i] = i;
541 prob->numberColsToDo_ = ncols_;
542 } else {
543 // some stuff must be left alone
544 prob->numberRowsToDo_ = 0;
545 for (i = 0; i < nrows_; i++)
546 if (!prob->rowProhibited(i))
547 prob->rowsToDo_[prob->numberRowsToDo_++] = i;
548 prob->numberColsToDo_ = 0;
549 for (i = 0; i < ncols_; i++)
550 if (!prob->colProhibited(i))
551 prob->colsToDo_[prob->numberColsToDo_++] = i;
552 }
553
554
555 int iLoop;
556#if PRESOLVE_DEBUG
557 check_sol(prob, 1.0e0);
558#endif
559 if (dupcol) {
560 // maybe allow integer columns to be checked
561 if ((presolveActions_ & 512) != 0)
562 prob->setPresolveOptions(prob->presolveOptions() | 1);
563 possibleSkip;
564 paction_ = dupcol_action::presolve(prob, paction_);
565 }
566 if (duprow) {
567 possibleSkip;
568 paction_ = duprow_action::presolve(prob, paction_);
569 }
570 if (doGubrow()) {
571 possibleSkip;
572 paction_ = gubrow_action::presolve(prob, paction_);
573 }
574
575 if ((presolveActions_ & 16384) != 0)
576 prob->setPresolveOptions(prob->presolveOptions() | 16384);
577 // Check number rows dropped
578 int lastDropped = 0;
579 prob->pass_ = 0;
580 if (numberPasses_<=5)
581 prob->presolveOptions_ |= 0x10000; // say more lightweight
582 for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
583 // See if we want statistics
584 if ((presolveActions_ & 0x80000000) != 0)
585 printf("Starting major pass %d after %g seconds\n", iLoop + 1, CoinCpuTime() - prob->startTime_);
586#ifdef PRESOLVE_SUMMARY
587 printf("Starting major pass %d\n", iLoop + 1);
588#endif
589 const CoinPresolveAction * const paction0 = paction_;
590 // look for substitutions with no fill
591 //#define IMPLIED 3
592#ifdef IMPLIED
593 int fill_level = 3;
594#define IMPLIED2 99
595#if IMPLIED!=3
596#if IMPLIED>2&&IMPLIED<11
597 fill_level = IMPLIED;
598 COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
599#endif
600#if IMPLIED>11&&IMPLIED<21
601 fill_level = -(IMPLIED - 10);
602 COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
603#endif
604#endif
605#else
606 int fill_level = 2;
607#endif
608 int whichPass = 0;
609 while (1) {
610 whichPass++;
611 prob->pass_++;
612 const CoinPresolveAction * const paction1 = paction_;
613
614 if (slackd) {
615 bool notFinished = true;
616 while (notFinished) {
617 possibleBreak;
618 paction_ = slack_doubleton_action::presolve(prob, paction_,
619 notFinished);
620 }
621 if (prob->status_)
622 break;
623 }
624 if (dual && whichPass == 1) {
625 // this can also make E rows so do one bit here
626 possibleBreak;
627 paction_ = remove_dual_action::presolve(prob, paction_);
628 if (prob->status_)
629 break;
630 }
631
632 if (doubleton) {
633 possibleBreak;
634 paction_ = doubleton_action::presolve(prob, paction_);
635 if (prob->status_)
636 break;
637 }
638 if (tripleton) {
639 possibleBreak;
640 paction_ = tripleton_action::presolve(prob, paction_);
641 if (prob->status_)
642 break;
643 }
644
645 if (zerocost) {
646 possibleBreak;
647 paction_ = do_tighten_action::presolve(prob, paction_);
648 if (prob->status_)
649 break;
650 }
651#ifndef NO_FORCING
652 if (forcing) {
653 possibleBreak;
654 paction_ = forcing_constraint_action::presolve(prob, paction_);
655 if (prob->status_)
656 break;
657 }
658#endif
659
660 if (ifree && (whichPass % 5) == 1) {
661 possibleBreak;
662 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
663 if (prob->status_)
664 break;
665 }
666
667#if PRESOLVE_DEBUG
668 check_sol(prob, 1.0e0);
669#endif
670
671#if PRESOLVE_CONSISTENCY
672// presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
673// prob->nrows_);
674 presolve_links_ok(prob, false, true) ;
675#endif
676
677//#if PRESOLVE_DEBUG
678// presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
679// prob->ncols_);
680//#endif
681//#if PRESOLVE_CONSISTENCY
682// prob->consistent();
683//#endif
684#if PRESOLVE_CONSISTENCY
685 presolve_no_zeros(prob, true, false) ;
686 presolve_consistent(prob, true) ;
687#endif
688
689 {
690 // set up for next pass
691 // later do faster if many changes i.e. memset and memcpy
692 const int * count = prob->hinrow_;
693 const int * nextToDo = prob->nextRowsToDo_;
694 int * toDo = prob->rowsToDo_;
695 int nNext = prob->numberNextRowsToDo_;
696 int n = 0;
697 for (int i = 0; i < nNext; i++) {
698 int index = nextToDo[i];
699 prob->unsetRowChanged(index);
700 if (count[index])
701 toDo[n++] = index;
702 }
703 prob->numberRowsToDo_ = n;
704 prob->numberNextRowsToDo_ = 0;
705 count = prob->hincol_;
706 nextToDo = prob->nextColsToDo_;
707 toDo = prob->colsToDo_;
708 nNext = prob->numberNextColsToDo_;
709 n = 0;
710 for (int i = 0; i < nNext; i++) {
711 int index = nextToDo[i];
712 prob->unsetColChanged(index);
713 if (count[index])
714 toDo[n++] = index;
715 }
716 prob->numberColsToDo_ = n;
717 prob->numberNextColsToDo_ = 0;
718 }
719 if (paction_ == paction1 && fill_level > 0)
720 break;
721 }
722 // say look at all
723 int i;
724 if (!prob->anyProhibited()) {
725 const int * count = prob->hinrow_;
726 int * toDo = prob->rowsToDo_;
727 int n = 0;
728 for (int i = 0; i < nrows_; i++) {
729 prob->unsetRowChanged(i);
730 if (count[i])
731 toDo[n++] = i;
732 }
733 prob->numberRowsToDo_ = n;
734 prob->numberNextRowsToDo_ = 0;
735 count = prob->hincol_;
736 toDo = prob->colsToDo_;
737 n = 0;
738 for (int i = 0; i < ncols_; i++) {
739 prob->unsetColChanged(i);
740 if (count[i])
741 toDo[n++] = i;
742 }
743 prob->numberColsToDo_ = n;
744 prob->numberNextColsToDo_ = 0;
745 } else {
746 // some stuff must be left alone
747 prob->numberRowsToDo_ = 0;
748 for (i = 0; i < nrows_; i++)
749 if (!prob->rowProhibited(i))
750 prob->rowsToDo_[prob->numberRowsToDo_++] = i;
751 prob->numberColsToDo_ = 0;
752 for (i = 0; i < ncols_; i++)
753 if (!prob->colProhibited(i))
754 prob->colsToDo_[prob->numberColsToDo_++] = i;
755 }
756 // now expensive things
757 // this caused world.mps to run into numerical difficulties
758#ifdef PRESOLVE_SUMMARY
759 printf("Starting expensive\n");
760#endif
761
762 if (dual) {
763 int itry;
764 for (itry = 0; itry < 5; itry++) {
765 possibleBreak;
766 paction_ = remove_dual_action::presolve(prob, paction_);
767 if (prob->status_)
768 break;
769 const CoinPresolveAction * const paction2 = paction_;
770 if (ifree) {
771#ifdef IMPLIED
772#if IMPLIED2 ==0
773 int fill_level = 0; // switches off substitution
774#elif IMPLIED2!=99
775 int fill_level = IMPLIED2;
776#endif
777#endif
778 if ((itry & 1) == 0) {
779 possibleBreak;
780 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
781 }
782 if (prob->status_)
783 break;
784 }
785 if (paction_ == paction2)
786 break;
787 }
788 } else if (ifree) {
789 // just free
790#ifdef IMPLIED
791#if IMPLIED2 ==0
792 int fill_level = 0; // switches off substitution
793#elif IMPLIED2!=99
794 int fill_level = IMPLIED2;
795#endif
796#endif
797 possibleBreak;
798 paction_ = implied_free_action::presolve(prob, paction_, fill_level);
799 if (prob->status_)
800 break;
801 }
802#if PRESOLVE_DEBUG
803 check_sol(prob, 1.0e0);
804#endif
805 if (dupcol) {
806 // maybe allow integer columns to be checked
807 if ((presolveActions_ & 512) != 0)
808 prob->setPresolveOptions(prob->presolveOptions() | 1);
809 possibleBreak;
810 paction_ = dupcol_action::presolve(prob, paction_);
811 if (prob->status_)
812 break;
813 }
814#if PRESOLVE_DEBUG
815 check_sol(prob, 1.0e0);
816#endif
817
818 if (duprow) {
819 possibleBreak;
820 paction_ = duprow_action::presolve(prob, paction_);
821 if (prob->status_)
822 break;
823 }
824#if PRESOLVE_DEBUG
825 check_sol(prob, 1.0e0);
826#endif
827 bool stopLoop = false;
828 {
829 int * hinrow = prob->hinrow_;
830 int numberDropped = 0;
831 for (i = 0; i < nrows_; i++)
832 if (!hinrow[i])
833 numberDropped++;
834
835 prob->messageHandler()->message(COIN_PRESOLVE_PASS,
836 messages)
837 << numberDropped << iLoop + 1
838 << CoinMessageEol;
839 //printf("%d rows dropped after pass %d\n",numberDropped,
840 // iLoop+1);
841 if (numberDropped == lastDropped)
842 stopLoop = true;
843 else
844 lastDropped = numberDropped;
845 }
846 // Do this here as not very loopy
847 if (slackSingleton) {
848 // On most passes do not touch costed slacks
849 if (paction_ != paction0 && !stopLoop) {
850 possibleBreak;
851 paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
852 } else {
853 // do costed if Clp (at end as ruins rest of presolve)
854 possibleBreak;
855 paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
856 stopLoop = true;
857 }
858 }
859#if PRESOLVE_DEBUG
860 check_sol(prob, 1.0e0);
861#endif
862 if (paction_ == paction0 || stopLoop)
863 break;
864
865 }
866 }
867 prob->presolveOptions_ &= ~0x10000;
868 if (!prob->status_) {
869 paction_ = drop_zero_coefficients(prob, paction_);
870#if PRESOLVE_DEBUG
871 check_sol(prob, 1.0e0);
872#endif
873
874 paction_ = drop_empty_cols_action::presolve(prob, paction_);
875 paction_ = drop_empty_rows_action::presolve(prob, paction_);
876#if PRESOLVE_DEBUG
877 check_sol(prob, 1.0e0);
878#endif
879 }
880
881 if (prob->status_) {
882 if (prob->status_ == 1)
883 prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
884 messages)
885 << prob->feasibilityTolerance_
886 << CoinMessageEol;
887 else if (prob->status_ == 2)
888 prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
889 messages)
890 << CoinMessageEol;
891 else
892 prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
893 messages)
894 << CoinMessageEol;
895 // get rid of data
896 destroyPresolve();
897 }
898 return (paction_);
899}
900
901void check_djs(CoinPostsolveMatrix *prob);
902
903
904// We could have implemented this by having each postsolve routine
905// directly call the next one, but this may make it easier to add debugging checks.
906void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
907{
908 {
909 // Check activities
910 double *colels = prob.colels_;
911 int *hrow = prob.hrow_;
912 CoinBigIndex *mcstrt = prob.mcstrt_;
913 int *hincol = prob.hincol_;
914 int *link = prob.link_;
915 int ncols = prob.ncols_;
916
917 char *cdone = prob.cdone_;
918
919 double * csol = prob.sol_;
920 int nrows = prob.nrows_;
921
922 int colx;
923
924 double * rsol = prob.acts_;
925 memset(rsol, 0, nrows * sizeof(double));
926
927 for (colx = 0; colx < ncols; ++colx) {
928 if (cdone[colx]) {
929 CoinBigIndex k = mcstrt[colx];
930 int nx = hincol[colx];
931 double solutionValue = csol[colx];
932 for (int i = 0; i < nx; ++i) {
933 int row = hrow[k];
934 double coeff = colels[k];
935 k = link[k];
936 rsol[row] += solutionValue * coeff;
937 }
938 }
939 }
940 }
941 const CoinPresolveAction *paction = paction_;
942 //#define PRESOLVE_DEBUG 1
943#if PRESOLVE_DEBUG
944 // Check only works after first one
945 int checkit = -1;
946#endif
947
948 while (paction) {
949#if PRESOLVE_DEBUG
950 printf("POSTSOLVING %s\n", paction->name());
951#endif
952
953 paction->postsolve(&prob);
954
955#if PRESOLVE_DEBUG
956 {
957 int nr = 0;
958 int i;
959 for (i = 0; i < prob.nrows_; i++) {
960 if ((prob.rowstat_[i] & 7) == 1) {
961 nr++;
962 } else if ((prob.rowstat_[i] & 7) == 2) {
963 // at ub
964 assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
965 } else if ((prob.rowstat_[i] & 7) == 3) {
966 // at lb
967 assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
968 }
969 }
970 int nc = 0;
971 for (i = 0; i < prob.ncols_; i++) {
972 if ((prob.colstat_[i] & 7) == 1)
973 nc++;
974 }
975 printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
976 }
977 checkit++;
978 if (prob.colstat_ && checkit > 0) {
979 presolve_check_nbasic(&prob) ;
980 presolve_check_sol(&prob, 2, 2, 1) ;
981 }
982#endif
983 paction = paction->next;
984#if PRESOLVE_DEBUG
985// check_djs(&prob);
986 if (checkit > 0)
987 presolve_check_reduced_costs(&prob) ;
988#endif
989 }
990#if PRESOLVE_DEBUG
991 if (prob.colstat_) {
992 presolve_check_nbasic(&prob) ;
993 presolve_check_sol(&prob, 2, 2, 1) ;
994 }
995#endif
996#undef PRESOLVE_DEBUG
997
998#if 0 && PRESOLVE_DEBUG
999 for (i = 0; i < ncols0; i++) {
1000 if (!cdone[i]) {
1001 printf("!cdone[%d]\n", i);
1002 abort();
1003 }
1004 }
1005
1006 for (i = 0; i < nrows0; i++) {
1007 if (!rdone[i]) {
1008 printf("!rdone[%d]\n", i);
1009 abort();
1010 }
1011 }
1012
1013
1014 for (i = 0; i < ncols0; i++) {
1015 if (sol[i] < -1e10 || sol[i] > 1e10)
1016 printf("!!!%d %g\n", i, sol[i]);
1017
1018 }
1019
1020
1021#endif
1022
1023#if 0 && PRESOLVE_DEBUG
1024 // debug check: make sure we ended up with same original matrix
1025 {
1026 int identical = 1;
1027
1028 for (int i = 0; i < ncols0; i++) {
1029 PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
1030 CoinBigIndex kcs0 = &prob->mcstrt0[i];
1031 CoinBigIndex kcs = mcstrt[i];
1032 int n = hincol[i];
1033 for (int k = 0; k < n; k++) {
1034 CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
1035
1036 if (k1 == kcs + n) {
1037 printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
1038 abort();
1039 }
1040
1041 if (colels[k1] != &prob->dels0[kcs0+k])
1042 printf("BAD COLEL[%d %d %d]: %g\n",
1043 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
1044
1045 if (kcs0 + k != k1)
1046 identical = 0;
1047 }
1048 }
1049 printf("identical? %d\n", identical);
1050 }
1051#endif
1052}
1053
1054
1055
1056
1057
1058
1059
1060
1061static inline double getTolerance(const ClpSimplex *si, ClpDblParam key)
1062{
1063 double tol;
1064 if (! si->getDblParam(key, tol)) {
1065 CoinPresolveAction::throwCoinError("getDblParam failed",
1066 "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
1067 }
1068 return (tol);
1069}
1070
1071
1072// Assumptions:
1073// 1. nrows>=m.getNumRows()
1074// 2. ncols>=m.getNumCols()
1075//
1076// In presolve, these values are equal.
1077// In postsolve, they may be inequal, since the reduced problem
1078// may be smaller, but we need room for the large problem.
1079// ncols may be larger than si.getNumCols() in postsolve,
1080// this at that point si will be the reduced problem,
1081// but we need to reserve enough space for the original problem.
1082CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
1083 int ncols_in,
1084 int nrows_in,
1085 CoinBigIndex nelems_in,
1086 double bulkRatio)
1087 : ncols_(si->getNumCols()),
1088 nrows_(si->getNumRows()),
1089 nelems_(si->getNumElements()),
1090 ncols0_(ncols_in),
1091 nrows0_(nrows_in),
1092 bulkRatio_(bulkRatio),
1093 mcstrt_(new CoinBigIndex[ncols_in+1]),
1094 hincol_(new int[ncols_in+1]),
1095 cost_(new double[ncols_in]),
1096 clo_(new double[ncols_in]),
1097 cup_(new double[ncols_in]),
1098 rlo_(new double[nrows_in]),
1099 rup_(new double[nrows_in]),
1100 originalColumn_(new int[ncols_in]),
1101 originalRow_(new int[nrows_in]),
1102 ztolzb_(getTolerance(si, ClpPrimalTolerance)),
1103 ztoldj_(getTolerance(si, ClpDualTolerance)),
1104 maxmin_(si->getObjSense()),
1105 sol_(NULL),
1106 rowduals_(NULL),
1107 acts_(NULL),
1108 rcosts_(NULL),
1109 colstat_(NULL),
1110 rowstat_(NULL),
1111 handler_(NULL),
1112 defaultHandler_(false),
1113 messages_()
1114
1115{
1116 bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
1117 hrow_ = new int [bulk0_];
1118 colels_ = new double[bulk0_];
1119 si->getDblParam(ClpObjOffset, originalOffset_);
1120 int ncols = si->getNumCols();
1121 int nrows = si->getNumRows();
1122
1123 setMessageHandler(si->messageHandler()) ;
1124
1125 ClpDisjointCopyN(si->getColLower(), ncols, clo_);
1126 ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
1127 //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1128 double offset;
1129 ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
1130 ClpDisjointCopyN(si->getRowLower(), nrows, rlo_);
1131 ClpDisjointCopyN(si->getRowUpper(), nrows, rup_);
1132 int i;
1133 for (i = 0; i < ncols_in; i++)
1134 originalColumn_[i] = i;
1135 for (i = 0; i < nrows_in; i++)
1136 originalRow_[i] = i;
1137 sol_ = NULL;
1138 rowduals_ = NULL;
1139 acts_ = NULL;
1140
1141 rcosts_ = NULL;
1142 colstat_ = NULL;
1143 rowstat_ = NULL;
1144}
1145
1146// I am not familiar enough with CoinPackedMatrix to be confident
1147// that I will implement a row-ordered version of toColumnOrderedGapFree
1148// properly.
1149static bool isGapFree(const CoinPackedMatrix& matrix)
1150{
1151 const CoinBigIndex * start = matrix.getVectorStarts();
1152 const int * length = matrix.getVectorLengths();
1153 int i = matrix.getSizeVectorLengths() - 1;
1154 // Quick check
1155 if (matrix.getNumElements() == start[i]) {
1156 return true;
1157 } else {
1158 for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1159 if (start[i+1] - start[i] != length[i])
1160 break;
1161 }
1162 return (! (i >= 0));
1163 }
1164}
1165#if PRESOLVE_DEBUG
1166static void matrix_bounds_ok(const double *lo, const double *up, int n)
1167{
1168 int i;
1169 for (i = 0; i < n; i++) {
1170 PRESOLVEASSERT(lo[i] <= up[i]);
1171 PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
1172 PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
1173 }
1174}
1175#endif
1176CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1177 double /*maxmin*/,
1178 // end prepost members
1179
1180 ClpSimplex * si,
1181
1182 // rowrep
1183 int nrows_in,
1184 CoinBigIndex nelems_in,
1185 bool doStatus,
1186 double nonLinearValue,
1187 double bulkRatio) :
1188
1189 CoinPrePostsolveMatrix(si,
1190 ncols0_in, nrows_in, nelems_in, bulkRatio),
1191 clink_(new presolvehlink[ncols0_in+1]),
1192 rlink_(new presolvehlink[nrows_in+1]),
1193
1194 dobias_(0.0),
1195
1196
1197 // temporary init
1198 integerType_(new unsigned char[ncols0_in]),
1199 tuning_(false),
1200 startTime_(0.0),
1201 feasibilityTolerance_(0.0),
1202 status_(-1),
1203 colsToDo_(new int [ncols0_in]),
1204 numberColsToDo_(0),
1205 nextColsToDo_(new int[ncols0_in]),
1206 numberNextColsToDo_(0),
1207 rowsToDo_(new int [nrows_in]),
1208 numberRowsToDo_(0),
1209 nextRowsToDo_(new int[nrows_in]),
1210 numberNextRowsToDo_(0),
1211 presolveOptions_(0)
1212{
1213 const int bufsize = bulk0_;
1214
1215 nrows_ = si->getNumRows() ;
1216
1217 // Set up change bits etc
1218 rowChanged_ = new unsigned char[nrows_];
1219 memset(rowChanged_, 0, nrows_);
1220 colChanged_ = new unsigned char[ncols_];
1221 memset(colChanged_, 0, ncols_);
1222 CoinPackedMatrix * m = si->matrix();
1223
1224 // The coefficient matrix is a big hunk of stuff.
1225 // Do the copy here to try to avoid running out of memory.
1226
1227 const CoinBigIndex * start = m->getVectorStarts();
1228 const int * row = m->getIndices();
1229 const double * element = m->getElements();
1230 int icol, nel = 0;
1231 mcstrt_[0] = 0;
1232 ClpDisjointCopyN(m->getVectorLengths(), ncols_, hincol_);
1233 for (icol = 0; icol < ncols_; icol++) {
1234 CoinBigIndex j;
1235 for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
1236 hrow_[nel] = row[j];
1237 if (fabs(element[j]) > ZTOLDP)
1238 colels_[nel++] = element[j];
1239 }
1240 mcstrt_[icol+1] = nel;
1241 hincol_[icol] = nel - mcstrt_[icol];
1242 }
1243
1244 // same thing for row rep
1245 CoinPackedMatrix * mRow = new CoinPackedMatrix();
1246 mRow->setExtraGap(0.0);
1247 mRow->setExtraMajor(0.0);
1248 mRow->reverseOrderedCopyOf(*m);
1249 //mRow->removeGaps();
1250 //mRow->setExtraGap(0.0);
1251
1252 // Now get rid of matrix
1253 si->createEmptyMatrix();
1254
1255 double * el = mRow->getMutableElements();
1256 int * ind = mRow->getMutableIndices();
1257 CoinBigIndex * strt = mRow->getMutableVectorStarts();
1258 int * len = mRow->getMutableVectorLengths();
1259 // Do carefully to save memory
1260 rowels_ = new double[bulk0_];
1261 ClpDisjointCopyN(el, nelems_, rowels_);
1262 mRow->nullElementArray();
1263 delete [] el;
1264 hcol_ = new int[bulk0_];
1265 ClpDisjointCopyN(ind, nelems_, hcol_);
1266 mRow->nullIndexArray();
1267 delete [] ind;
1268 mrstrt_ = new CoinBigIndex[nrows_in+1];
1269 ClpDisjointCopyN(strt, nrows_, mrstrt_);
1270 mRow->nullStartArray();
1271 mrstrt_[nrows_] = nelems_;
1272 delete [] strt;
1273 hinrow_ = new int[nrows_in+1];
1274 ClpDisjointCopyN(len, nrows_, hinrow_);
1275 if (nelems_ > nel) {
1276 nelems_ = nel;
1277 // Clean any small elements
1278 int irow;
1279 nel = 0;
1280 CoinBigIndex start = 0;
1281 for (irow = 0; irow < nrows_; irow++) {
1282 CoinBigIndex j;
1283 for (j = start; j < start + hinrow_[irow]; j++) {
1284 hcol_[nel] = hcol_[j];
1285 if (fabs(rowels_[j]) > ZTOLDP)
1286 rowels_[nel++] = rowels_[j];
1287 }
1288 start = mrstrt_[irow+1];
1289 mrstrt_[irow+1] = nel;
1290 hinrow_[irow] = nel - mrstrt_[irow];
1291 }
1292 }
1293
1294 delete mRow;
1295 if (si->integerInformation()) {
1296 CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()), ncols_, integerType_);
1297 } else {
1298 ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
1299 }
1300
1301#ifndef SLIM_CLP
1302#ifndef NO_RTTI
1303 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1304#else
1305 ClpQuadraticObjective * quadraticObj = NULL;
1306 if (si->objectiveAsObject()->type() == 2)
1307 quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
1308#endif
1309#endif
1310 // Set up prohibited bits if needed
1311 if (nonLinearValue) {
1312 anyProhibited_ = true;
1313 for (icol = 0; icol < ncols_; icol++) {
1314 int j;
1315 bool nonLinearColumn = false;
1316 if (cost_[icol] == nonLinearValue)
1317 nonLinearColumn = true;
1318 for (j = mcstrt_[icol]; j < mcstrt_[icol+1]; j++) {
1319 if (colels_[j] == nonLinearValue) {
1320 nonLinearColumn = true;
1321 setRowProhibited(hrow_[j]);
1322 }
1323 }
1324 if (nonLinearColumn)
1325 setColProhibited(icol);
1326 }
1327#ifndef SLIM_CLP
1328 } else if (quadraticObj) {
1329 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1330 //const int * columnQuadratic = quadratic->getIndices();
1331 //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1332 const int * columnQuadraticLength = quadratic->getVectorLengths();
1333 //double * quadraticElement = quadratic->getMutableElements();
1334 int numberColumns = quadratic->getNumCols();
1335 anyProhibited_ = true;
1336 for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1337 if (columnQuadraticLength[iColumn]) {
1338 setColProhibited(iColumn);
1339 //printf("%d prohib\n",iColumn);
1340 }
1341 }
1342#endif
1343 } else {
1344 anyProhibited_ = false;
1345 }
1346
1347 if (doStatus) {
1348 // allow for status and solution
1349 sol_ = new double[ncols_];
1350 CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);
1351 acts_ = new double [nrows_];
1352 CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
1353 if (!si->statusArray())
1354 si->createStatus();
1355 colstat_ = new unsigned char [nrows_+ncols_];
1356 CoinMemcpyN(si->statusArray(), (nrows_ + ncols_), colstat_);
1357 rowstat_ = colstat_ + ncols_;
1358 }
1359
1360 // the original model's fields are now unneeded - free them
1361
1362 si->resize(0, 0);
1363
1364#if PRESOLVE_DEBUG
1365 matrix_bounds_ok(rlo_, rup_, nrows_);
1366 matrix_bounds_ok(clo_, cup_, ncols_);
1367#endif
1368
1369#if 0
1370 for (i = 0; i < nrows; ++i)
1371 printf("NR: %6d\n", hinrow[i]);
1372 for (int i = 0; i < ncols; ++i)
1373 printf("NC: %6d\n", hincol[i]);
1374#endif
1375
1376 presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
1377 presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
1378
1379 // this allows last col/row to expand up to bufsize-1 (22);
1380 // this must come after the calls to presolve_prefix
1381 mcstrt_[ncols_] = bufsize - 1;
1382 mrstrt_[nrows_] = bufsize - 1;
1383 // Allocate useful arrays
1384 initializeStuff();
1385
1386#if PRESOLVE_CONSISTENCY
1387//consistent(false);
1388 presolve_consistent(this, false) ;
1389#endif
1390}
1391
1392void CoinPresolveMatrix::update_model(ClpSimplex * si,
1393 int /*nrows0*/,
1394 int /*ncols0*/,
1395 CoinBigIndex /*nelems0*/)
1396{
1397 si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
1398 clo_, cup_, cost_, rlo_, rup_);
1399 //delete [] si->integerInformation();
1400 int numberIntegers = 0;
1401 for (int i = 0; i < ncols_; i++) {
1402 if (integerType_[i])
1403 numberIntegers++;
1404 }
1405 if (numberIntegers)
1406 si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
1407 else
1408 si->copyInIntegerInformation(NULL);
1409
1410#if PRESOLVE_SUMMARY
1411 printf("NEW NCOL/NROW/NELS: %d(-%d) %d(-%d) %d(-%d)\n",
1412 ncols_, ncols0 - ncols_,
1413 nrows_, nrows0 - nrows_,
1414 si->getNumElements(), nelems0 - si->getNumElements());
1415#endif
1416 si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
1417
1418}
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430//////////////// POSTSOLVE
1431
1432CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex* si,
1433 int ncols0_in,
1434 int nrows0_in,
1435 CoinBigIndex nelems0,
1436
1437 double maxmin,
1438 // end prepost members
1439
1440 double *sol_in,
1441 double *acts_in,
1442
1443 unsigned char *colstat_in,
1444 unsigned char *rowstat_in) :
1445 CoinPrePostsolveMatrix(si,
1446 ncols0_in, nrows0_in, nelems0, 2.0),
1447
1448 free_list_(0),
1449 // link, free_list, maxlink
1450 maxlink_(bulk0_),
1451 link_(new int[/*maxlink*/ bulk0_]),
1452
1453 cdone_(new char[ncols0_]),
1454 rdone_(new char[nrows0_in])
1455
1456{
1457 bulk0_ = maxlink_ ;
1458 nrows_ = si->getNumRows() ;
1459 ncols_ = si->getNumCols() ;
1460
1461 sol_ = sol_in;
1462 rowduals_ = NULL;
1463 acts_ = acts_in;
1464
1465 rcosts_ = NULL;
1466 colstat_ = colstat_in;
1467 rowstat_ = rowstat_in;
1468
1469 // this is the *reduced* model, which is probably smaller
1470 int ncols1 = ncols_ ;
1471 int nrows1 = nrows_ ;
1472
1473 const CoinPackedMatrix * m = si->matrix();
1474
1475 const CoinBigIndex nelemsr = m->getNumElements();
1476 if (m->getNumElements() && !isGapFree(*m)) {
1477 // Odd - gaps
1478 CoinPackedMatrix mm(*m);
1479 mm.removeGaps();
1480 mm.setExtraGap(0.0);
1481
1482 ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
1483 CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1484 mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
1485 ClpDisjointCopyN(mm.getVectorLengths(), ncols1, hincol_);
1486 ClpDisjointCopyN(mm.getIndices(), nelemsr, hrow_);
1487 ClpDisjointCopyN(mm.getElements(), nelemsr, colels_);
1488 } else {
1489 // No gaps
1490
1491 ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1492 CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
1493 mcstrt_[ncols1] = nelems0; // ?? (should point to end of bulk store -- lh --)
1494 ClpDisjointCopyN(m->getVectorLengths(), ncols1, hincol_);
1495 ClpDisjointCopyN(m->getIndices(), nelemsr, hrow_);
1496 ClpDisjointCopyN(m->getElements(), nelemsr, colels_);
1497 }
1498
1499
1500
1501#if 0 && PRESOLVE_DEBUG
1502 presolve_check_costs(model, &colcopy);
1503#endif
1504
1505 // This determines the size of the data structure that contains
1506 // the matrix being postsolved. Links are taken from the free_list
1507 // to recreate matrix entries that were presolved away,
1508 // and links are added to the free_list when entries created during
1509 // presolve are discarded. There is never a need to gc this list.
1510 // Naturally, it should contain
1511 // exactly nelems0 entries "in use" when postsolving is done,
1512 // but I don't know whether the matrix could temporarily get
1513 // larger during postsolving. Substitution into more than two
1514 // rows could do that, in principle. I am being very conservative
1515 // here by reserving much more than the amount of space I probably need.
1516 // If this guess is wrong, check_free_list may be called.
1517 // int bufsize = 2*nelems0;
1518
1519 memset(cdone_, -1, ncols0_);
1520 memset(rdone_, -1, nrows0_);
1521
1522 rowduals_ = new double[nrows0_];
1523 ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1524
1525 rcosts_ = new double[ncols0_];
1526 ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1527 if (maxmin < 0.0) {
1528 // change so will look as if minimize
1529 int i;
1530 for (i = 0; i < nrows1; i++)
1531 rowduals_[i] = - rowduals_[i];
1532 for (i = 0; i < ncols1; i++) {
1533 rcosts_[i] = - rcosts_[i];
1534 }
1535 }
1536
1537 //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
1538 //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
1539
1540 ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
1541 si->setDblParam(ClpObjOffset, originalOffset_);
1542
1543 for (int j = 0; j < ncols1; j++) {
1544 CoinBigIndex kcs = mcstrt_[j];
1545 CoinBigIndex kce = kcs + hincol_[j];
1546 for (CoinBigIndex k = kcs; k < kce; ++k) {
1547 link_[k] = k + 1;
1548 }
1549 link_[kce-1] = NO_LINK ;
1550 }
1551 {
1552 int ml = maxlink_;
1553 for (CoinBigIndex k = nelemsr; k < ml; ++k)
1554 link_[k] = k + 1;
1555 if (ml)
1556 link_[ml-1] = NO_LINK;
1557 }
1558 free_list_ = nelemsr;
1559# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1560 /*
1561 These are used to track the action of postsolve transforms during debugging.
1562 */
1563 CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED) ;
1564 CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1) ;
1565 CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED) ;
1566 CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1) ;
1567# endif
1568}
1569/* This is main part of Presolve */
1570ClpSimplex *
1571ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
1572 double feasibilityTolerance,
1573 bool keepIntegers,
1574 int numberPasses,
1575 bool dropNames,
1576 bool doRowObjective)
1577{
1578 ncols_ = originalModel->getNumCols();
1579 nrows_ = originalModel->getNumRows();
1580 nelems_ = originalModel->getNumElements();
1581 numberPasses_ = numberPasses;
1582
1583 double maxmin = originalModel->getObjSense();
1584 originalModel_ = originalModel;
1585 delete [] originalColumn_;
1586 originalColumn_ = new int[ncols_];
1587 delete [] originalRow_;
1588 originalRow_ = new int[nrows_];
1589 // and fill in case returns early
1590 int i;
1591 for (i = 0; i < ncols_; i++)
1592 originalColumn_[i] = i;
1593 for (i = 0; i < nrows_; i++)
1594 originalRow_[i] = i;
1595 delete [] rowObjective_;
1596 if (doRowObjective) {
1597 rowObjective_ = new double [nrows_];
1598 memset(rowObjective_, 0, nrows_ * sizeof(double));
1599 } else {
1600 rowObjective_ = NULL;
1601 }
1602
1603 // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
1604 int result = -1;
1605
1606 // User may have deleted - its their responsibility
1607 presolvedModel_ = NULL;
1608 // Messages
1609 CoinMessages messages = originalModel->coinMessages();
1610 // Only go round 100 times even if integer preprocessing
1611 int totalPasses = 100;
1612 while (result == -1) {
1613
1614#ifndef CLP_NO_STD
1615 // make new copy
1616 if (saveFile_ == "") {
1617#endif
1618 delete presolvedModel_;
1619#ifndef CLP_NO_STD
1620 // So won't get names
1621 int lengthNames = originalModel->lengthNames();
1622 originalModel->setLengthNames(0);
1623#endif
1624 presolvedModel_ = new ClpSimplex(*originalModel);
1625#ifndef CLP_NO_STD
1626 originalModel->setLengthNames(lengthNames);
1627 presolvedModel_->dropNames();
1628 } else {
1629 presolvedModel_ = originalModel;
1630 if (dropNames)
1631 presolvedModel_->dropNames();
1632 }
1633#endif
1634
1635 // drop integer information if wanted
1636 if (!keepIntegers)
1637 presolvedModel_->deleteIntegerInformation();
1638 totalPasses--;
1639
1640 double ratio = 2.0;
1641 if (substitution_ > 3)
1642 ratio = substitution_;
1643 else if (substitution_ == 2)
1644 ratio = 1.5;
1645 CoinPresolveMatrix prob(ncols_,
1646 maxmin,
1647 presolvedModel_,
1648 nrows_, nelems_, true, nonLinearValue_, ratio);
1649 prob.setMaximumSubstitutionLevel(substitution_);
1650 if (doRowObjective)
1651 memset(rowObjective_, 0, nrows_ * sizeof(double));
1652 // See if we want statistics
1653 if ((presolveActions_ & 0x80000000) != 0)
1654 prob.statistics();
1655 // make sure row solution correct
1656 {
1657 double *colels = prob.colels_;
1658 int *hrow = prob.hrow_;
1659 CoinBigIndex *mcstrt = prob.mcstrt_;
1660 int *hincol = prob.hincol_;
1661 int ncols = prob.ncols_;
1662
1663
1664 double * csol = prob.sol_;
1665 double * acts = prob.acts_;
1666 int nrows = prob.nrows_;
1667
1668 int colx;
1669
1670 memset(acts, 0, nrows * sizeof(double));
1671
1672 for (colx = 0; colx < ncols; ++colx) {
1673 double solutionValue = csol[colx];
1674 for (int i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
1675 int row = hrow[i];
1676 double coeff = colels[i];
1677 acts[row] += solutionValue * coeff;
1678 }
1679 }
1680 }
1681
1682 // move across feasibility tolerance
1683 prob.feasibilityTolerance_ = feasibilityTolerance;
1684
1685 // Do presolve
1686 paction_ = presolve(&prob);
1687 // Get rid of useful arrays
1688 prob.deleteStuff();
1689
1690 result = 0;
1691
1692 bool fixInfeasibility = (prob.presolveOptions_&16384)!=0;
1693 bool hasSolution = (prob.presolveOptions_&32768)!=0;
1694 if (prob.status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
1695 // Looks feasible but double check to see if anything slipped through
1696 int n = prob.ncols_;
1697 double * lo = prob.clo_;
1698 double * up = prob.cup_;
1699 int i;
1700
1701 for (i = 0; i < n; i++) {
1702 if (up[i] < lo[i]) {
1703 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
1704 // infeasible
1705 prob.status_ = 1;
1706 } else {
1707 up[i] = lo[i];
1708 }
1709 }
1710 }
1711
1712 n = prob.nrows_;
1713 lo = prob.rlo_;
1714 up = prob.rup_;
1715
1716 for (i = 0; i < n; i++) {
1717 if (up[i] < lo[i]) {
1718 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
1719 // infeasible
1720 prob.status_ = 1;
1721 } else {
1722 up[i] = lo[i];
1723 }
1724 }
1725 }
1726 }
1727 if (prob.status_ == 0 && paction_) {
1728 // feasible
1729
1730 prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
1731 // copy status and solution
1732 CoinMemcpyN( prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
1733 CoinMemcpyN( prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
1734 CoinMemcpyN( prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
1735 CoinMemcpyN( prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
1736 if (fixInfeasibility && hasSolution) {
1737 // Looks feasible but double check to see if anything slipped through
1738 int n = prob.ncols_;
1739 double * lo = prob.clo_;
1740 double * up = prob.cup_;
1741 double * rsol = prob.acts_;
1742 //memset(prob.acts_,0,prob.nrows_*sizeof(double));
1743 presolvedModel_->matrix()->times(prob.sol_,rsol);
1744 int i;
1745
1746 for (i = 0; i < n; i++) {
1747 double gap=up[i]-lo[i];
1748 if (rsol[i]<lo[i]-feasibilityTolerance&&fabs(rsol[i]-lo[i])<1.0e-3) {
1749 lo[i]=rsol[i];
1750 if (gap<1.0e5)
1751 up[i]=lo[i]+gap;
1752 } else if (rsol[i]>up[i]+feasibilityTolerance&&fabs(rsol[i]-up[i])<1.0e-3) {
1753 up[i]=rsol[i];
1754 if (gap<1.0e5)
1755 lo[i]=up[i]-gap;
1756 }
1757 if (up[i] < lo[i]) {
1758 up[i] = lo[i];
1759 }
1760 }
1761 }
1762
1763 int n = prob.nrows_;
1764 double * lo = prob.rlo_;
1765 double * up = prob.rup_;
1766
1767 for (i = 0; i < n; i++) {
1768 if (up[i] < lo[i]) {
1769 if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
1770 // infeasible
1771 prob.status_ = 1;
1772 } else {
1773 up[i] = lo[i];
1774 }
1775 }
1776 }
1777 delete [] prob.sol_;
1778 delete [] prob.acts_;
1779 delete [] prob.colstat_;
1780 prob.sol_ = NULL;
1781 prob.acts_ = NULL;
1782 prob.colstat_ = NULL;
1783
1784 int ncolsNow = presolvedModel_->getNumCols();
1785 CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
1786#ifndef SLIM_CLP
1787#ifndef NO_RTTI
1788 ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1789#else
1790 ClpQuadraticObjective * quadraticObj = NULL;
1791 if (originalModel->objectiveAsObject()->type() == 2)
1792 quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
1793#endif
1794 if (quadraticObj) {
1795 // set up for subset
1796 char * mark = new char [ncols_];
1797 memset(mark, 0, ncols_);
1798 CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
1799 //const int * columnQuadratic = quadratic->getIndices();
1800 //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
1801 const int * columnQuadraticLength = quadratic->getVectorLengths();
1802 //double * quadraticElement = quadratic->getMutableElements();
1803 int numberColumns = quadratic->getNumCols();
1804 ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
1805 ncolsNow,
1806 originalColumn_);
1807 // and modify linear and check
1808 double * linear = newObj->linearObjective();
1809 CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
1810 int iColumn;
1811 for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
1812 if (columnQuadraticLength[iColumn])
1813 mark[iColumn] = 1;
1814 }
1815 // and new
1816 quadratic = newObj->quadraticObjective();
1817 columnQuadraticLength = quadratic->getVectorLengths();
1818 int numberColumns2 = quadratic->getNumCols();
1819 for ( iColumn = 0; iColumn < numberColumns2; iColumn++) {
1820 if (columnQuadraticLength[iColumn])
1821 mark[originalColumn_[iColumn]] = 0;
1822 }
1823 presolvedModel_->setObjective(newObj);
1824 delete newObj;
1825 // final check
1826 for ( iColumn = 0; iColumn < numberColumns; iColumn++)
1827 if (mark[iColumn])
1828 printf("Quadratic column %d modified - may be okay\n", iColumn);
1829 delete [] mark;
1830 }
1831#endif
1832 delete [] prob.originalColumn_;
1833 prob.originalColumn_ = NULL;
1834 int nrowsNow = presolvedModel_->getNumRows();
1835 CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
1836 delete [] prob.originalRow_;
1837 prob.originalRow_ = NULL;
1838#ifndef CLP_NO_STD
1839 if (!dropNames && originalModel->lengthNames()) {
1840 // Redo names
1841 int iRow;
1842 std::vector<std::string> rowNames;
1843 rowNames.reserve(nrowsNow);
1844 for (iRow = 0; iRow < nrowsNow; iRow++) {
1845 int kRow = originalRow_[iRow];
1846 rowNames.push_back(originalModel->rowName(kRow));
1847 }
1848
1849 int iColumn;
1850 std::vector<std::string> columnNames;
1851 columnNames.reserve(ncolsNow);
1852 for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
1853 int kColumn = originalColumn_[iColumn];
1854 columnNames.push_back(originalModel->columnName(kColumn));
1855 }
1856 presolvedModel_->copyNames(rowNames, columnNames);
1857 } else {
1858 presolvedModel_->setLengthNames(0);
1859 }
1860#endif
1861 if (rowObjective_) {
1862 int iRow;
1863 int k = -1;
1864 int nObj = 0;
1865 for (iRow = 0; iRow < nrowsNow; iRow++) {
1866 int kRow = originalRow_[iRow];
1867 assert (kRow > k);
1868 k = kRow;
1869 rowObjective_[iRow] = rowObjective_[kRow];
1870 if (rowObjective_[iRow])
1871 nObj++;
1872 }
1873 if (nObj) {
1874 printf("%d costed slacks\n", nObj);
1875 presolvedModel_->setRowObjective(rowObjective_);
1876 }
1877 }
1878 /* now clean up integer variables. This can modify original
1879 Don't do if dupcol added columns together */
1880 int i;
1881 const char * information = presolvedModel_->integerInformation();
1882 if ((prob.presolveOptions_ & 0x80000000) == 0 && information) {
1883 int numberChanges = 0;
1884 double * lower0 = originalModel_->columnLower();
1885 double * upper0 = originalModel_->columnUpper();
1886 double * lower = presolvedModel_->columnLower();
1887 double * upper = presolvedModel_->columnUpper();
1888 for (i = 0; i < ncolsNow; i++) {
1889 if (!information[i])
1890 continue;
1891 int iOriginal = originalColumn_[i];
1892 double lowerValue0 = lower0[iOriginal];
1893 double upperValue0 = upper0[iOriginal];
1894 double lowerValue = ceil(lower[i] - 1.0e-5);
1895 double upperValue = floor(upper[i] + 1.0e-5);
1896 lower[i] = lowerValue;
1897 upper[i] = upperValue;
1898 if (lowerValue > upperValue) {
1899 numberChanges++;
1900 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
1901 messages)
1902 << iOriginal
1903 << lowerValue
1904 << upperValue
1905 << CoinMessageEol;
1906 result = 1;
1907 } else {
1908 if (lowerValue > lowerValue0 + 1.0e-8) {
1909 lower0[iOriginal] = lowerValue;
1910 numberChanges++;
1911 }
1912 if (upperValue < upperValue0 - 1.0e-8) {
1913 upper0[iOriginal] = upperValue;
1914 numberChanges++;
1915 }
1916 }
1917 }
1918 if (numberChanges) {
1919 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
1920 messages)
1921 << numberChanges
1922 << CoinMessageEol;
1923 if (!result && totalPasses > 0) {
1924 result = -1; // round again
1925 const CoinPresolveAction *paction = paction_;
1926 while (paction) {
1927 const CoinPresolveAction *next = paction->next;
1928 delete paction;
1929 paction = next;
1930 }
1931 paction_ = NULL;
1932 }
1933 }
1934 }
1935 } else if (prob.status_) {
1936 // infeasible or unbounded
1937 result = 1;
1938 // Put status in nelems_!
1939 nelems_ = - prob.status_;
1940 originalModel->setProblemStatus(prob.status_);
1941 } else {
1942 // no changes - model needs restoring after Lou's changes
1943#ifndef CLP_NO_STD
1944 if (saveFile_ == "") {
1945#endif
1946 delete presolvedModel_;
1947 presolvedModel_ = new ClpSimplex(*originalModel);
1948 // but we need to remove gaps
1949 ClpPackedMatrix* clpMatrix =
1950 dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
1951 if (clpMatrix) {
1952 clpMatrix->getPackedMatrix()->removeGaps();
1953 }
1954#ifndef CLP_NO_STD
1955 } else {
1956 presolvedModel_ = originalModel;
1957 }
1958 presolvedModel_->dropNames();
1959#endif
1960
1961 // drop integer information if wanted
1962 if (!keepIntegers)
1963 presolvedModel_->deleteIntegerInformation();
1964 result = 2;
1965 }
1966 }
1967 if (result == 0 || result == 2) {
1968 int nrowsAfter = presolvedModel_->getNumRows();
1969 int ncolsAfter = presolvedModel_->getNumCols();
1970 CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
1971 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
1972 messages)
1973 << nrowsAfter << -(nrows_ - nrowsAfter)
1974 << ncolsAfter << -(ncols_ - ncolsAfter)
1975 << nelsAfter << -(nelems_ - nelsAfter)
1976 << CoinMessageEol;
1977 } else {
1978 destroyPresolve();
1979 if (presolvedModel_ != originalModel_)
1980 delete presolvedModel_;
1981 presolvedModel_ = NULL;
1982 }
1983 return presolvedModel_;
1984}
1985