1// Copyright (C) 2002, International Business Machines
2// Corporation and others. All Rights Reserved.
3// This code is licensed under the terms of the Eclipse Public License (EPL).
4
5/*
6 Debug compile symbols for CoinPresolve.
7
8 DEFINE THE SAME SET OF SYMBOLS when building the various CoinPresolve*.cpp
9 files in CoinUtils. Without consistent symbol definitions, the results will
10 be somewhere between garbage and a core dump.
11*/
12
13// #define PRESOLVE_CONSISTENCY 1
14// #define PRESOLVE_DEBUG 1
15// #define PRESOLVE_SUMMARY 1
16
17#include <stdio.h>
18
19#include <cassert>
20#include <iostream>
21
22#include "CoinHelperFunctions.hpp"
23#include "CoinFinite.hpp"
24
25#include "CoinPackedMatrix.hpp"
26#include "CoinWarmStartBasis.hpp"
27#include "OsiSolverInterface.hpp"
28
29#include "OsiPresolve.hpp"
30#include "CoinPresolveMatrix.hpp"
31
32#include "CoinPresolveEmpty.hpp"
33#include "CoinPresolveFixed.hpp"
34#include "CoinPresolvePsdebug.hpp"
35#include "CoinPresolveSingleton.hpp"
36#include "CoinPresolveDoubleton.hpp"
37#include "CoinPresolveTripleton.hpp"
38#include "CoinPresolveZeros.hpp"
39#include "CoinPresolveSubst.hpp"
40#include "CoinPresolveForcing.hpp"
41#include "CoinPresolveDual.hpp"
42#include "CoinPresolveTighten.hpp"
43#include "CoinPresolveUseless.hpp"
44#include "CoinPresolveDupcol.hpp"
45#include "CoinPresolveImpliedFree.hpp"
46#include "CoinPresolveIsolated.hpp"
47#include "CoinMessage.hpp"
48
49
50OsiPresolve::OsiPresolve() :
51 originalModel_(NULL),
52 presolvedModel_(NULL),
53 nonLinearValue_(0.0),
54 originalColumn_(NULL),
55 originalRow_(NULL),
56 paction_(0),
57 ncols_(0),
58 nrows_(0),
59 nelems_(0),
60 presolveActions_(0),
61 numberPasses_(5)
62{
63}
64
65OsiPresolve::~OsiPresolve()
66{
67 gutsOfDestroy();
68}
69// Gets rid of presolve actions (e.g.when infeasible)
70void
71OsiPresolve::gutsOfDestroy()
72{
73 const CoinPresolveAction *paction = paction_;
74 while (paction) {
75 const CoinPresolveAction *next = paction->next;
76 delete paction;
77 paction = next;
78 }
79 delete [] originalColumn_;
80 delete [] originalRow_;
81 paction_=NULL;
82 originalColumn_=NULL;
83 originalRow_=NULL;
84}
85
86/* This version of presolve returns a pointer to a new presolved
87 model. NULL if infeasible
88
89 doStatus controls activities required to transform an existing
90 solution to match the presolved problem. I'd (lh) argue that this should
91 default to false, but to maintain previous behaviour it defaults to true.
92 Really, this is only useful if you've already optimised before applying
93 presolve and also want to work with the solution after presolve. I think
94 that this is the less common case. The more common situation is to apply
95 presolve before optimising.
96*/
97OsiSolverInterface *
98OsiPresolve::presolvedModel(OsiSolverInterface & si,
99 double feasibilityTolerance,
100 bool keepIntegers,
101 int numberPasses,
102 const char * prohibited,
103 bool doStatus,
104 const char * rowProhibited)
105{
106 ncols_ = si.getNumCols();
107 nrows_ = si.getNumRows();
108 nelems_ = si.getNumElements();
109 numberPasses_ = numberPasses;
110
111 double maxmin = si.getObjSense();
112 originalModel_ = &si;
113 delete [] originalColumn_;
114 originalColumn_ = new int[ncols_];
115 delete [] originalRow_;
116 originalRow_ = new int[nrows_];
117 int i;
118 for (i=0;i<ncols_;i++)
119 originalColumn_[i]=i;
120 for (i=0;i<nrows_;i++)
121 originalRow_[i]=i;
122
123 // result is 0 - okay, 1 infeasible, -1 go round again
124 int result = -1;
125
126 // User may have deleted - its their responsibility
127 presolvedModel_=NULL;
128 // Messages
129 CoinMessages messages = CoinMessage(si.messages().language());
130 // Only go round 100 times even if integer preprocessing
131 int totalPasses=100;
132 while (result==-1) {
133
134 // make new copy
135 delete presolvedModel_;
136 presolvedModel_ = si.clone();
137 totalPasses--;
138
139 // drop integer information if wanted
140 if (!keepIntegers) {
141 int i;
142 for (i=0;i<ncols_;i++)
143 presolvedModel_->setContinuous(i);
144 }
145
146
147 CoinPresolveMatrix prob(ncols_,
148 maxmin,
149 presolvedModel_,
150 nrows_, nelems_,doStatus,nonLinearValue_,prohibited,
151 rowProhibited);
152 // make sure row solution correct
153 if (doStatus) {
154 double *colels = prob.colels_;
155 int *hrow = prob.hrow_;
156 CoinBigIndex *mcstrt = prob.mcstrt_;
157 int *hincol = prob.hincol_;
158 int ncols = prob.ncols_;
159
160
161 double * csol = prob.sol_;
162 double * acts = prob.acts_;
163 int nrows = prob.nrows_;
164
165 int colx;
166
167 memset(acts,0,nrows*sizeof(double));
168
169 for (colx = 0; colx < ncols; ++colx) {
170 double solutionValue = csol[colx];
171 for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
172 int row = hrow[i];
173 double coeff = colels[i];
174 acts[row] += solutionValue*coeff;
175 }
176 }
177 }
178
179 // move across feasibility tolerance
180 prob.feasibilityTolerance_ = feasibilityTolerance;
181
182/*
183 Do presolve. Allow for the possibility that presolve might be ineffective
184 (i.e., we're feasible but no postsolve actions are queued.
185*/
186 paction_ = presolve(&prob) ;
187 result = 0 ;
188 // Get rid of useful arrays
189 prob.deleteStuff();
190/*
191 This we don't need to do unless presolve actually reduced the system.
192*/
193 if (prob.status_==0&&paction_) {
194 // Looks feasible but double check to see if anything slipped through
195 int n = prob.ncols_;
196 double * lo = prob.clo_;
197 double * up = prob.cup_;
198 int i;
199
200 for (i=0;i<n;i++) {
201 if (up[i]<lo[i]) {
202 if (up[i]<lo[i]-1.0e-8) {
203 // infeasible
204 prob.status_=1;
205 } else {
206 up[i]=lo[i];
207 }
208 }
209 }
210
211 n = prob.nrows_;
212 lo = prob.rlo_;
213 up = prob.rup_;
214
215 for (i=0;i<n;i++) {
216 if (up[i]<lo[i]) {
217 if (up[i]<lo[i]-1.0e-8) {
218 // infeasible
219 prob.status_=1;
220 } else {
221 up[i]=lo[i];
222 }
223 }
224 }
225 }
226/*
227 If we're feasible, load the presolved system into the solver. Presumably we
228 could skip model update and copying of status and solution if presolve took
229 no action.
230*/
231 if (prob.status_ == 0) {
232
233 prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
234
235# if PRESOLVE_CONSISTENCY
236 if (doStatus)
237 { int basicCnt = 0 ;
238 int basicColumns = 0;
239 int i ;
240 CoinPresolveMatrix::Status status ;
241 for (i = 0 ; i < prob.ncols_ ; i++)
242 { status = prob.getColumnStatus(i);
243 if (status == CoinPrePostsolveMatrix::basic) basicColumns++ ; }
244 basicCnt = basicColumns;
245 for (i = 0 ; i < prob.nrows_ ; i++)
246 { status = prob.getRowStatus(i);
247 if (status == CoinPrePostsolveMatrix::basic) basicCnt++ ; }
248
249# if PRESOLVE_DEBUG
250 presolve_check_nbasic(&prob) ;
251# endif
252 if (basicCnt>prob.nrows_) {
253 // Take out slacks
254 double * acts = prob.acts_;
255 double * rlo = prob.rlo_;
256 double * rup = prob.rup_;
257 double infinity = si.getInfinity();
258 for (i = 0 ; i < prob.nrows_ ; i++) {
259 status = prob.getRowStatus(i);
260 if (status == CoinPrePostsolveMatrix::basic) {
261 basicCnt-- ;
262 double down = acts[i]-rlo[i];
263 double up = rup[i]-acts[i];
264 if (CoinMin(up,down)<infinity) {
265 if (down<=up)
266 prob.setRowStatus(i,CoinPrePostsolveMatrix::atLowerBound);
267 else
268 prob.setRowStatus(i,CoinPrePostsolveMatrix::atUpperBound);
269 } else {
270 prob.setRowStatus(i,CoinPrePostsolveMatrix::isFree);
271 }
272 }
273 if (basicCnt==prob.nrows_)
274 break;
275 }
276 }
277 }
278#endif
279
280/*
281 Install the status and primal solution, if we've been carrying them along.
282
283 The code that copies status is efficient but brittle. The current definitions
284 for CoinWarmStartBasis::Status and CoinPrePostsolveMatrix::Status are in
285 one-to-one correspondence. This code will fail if that ever changes.
286*/
287 if (doStatus) {
288 presolvedModel_->setColSolution(prob.sol_);
289 CoinWarmStartBasis *basis =
290 dynamic_cast<CoinWarmStartBasis *>(presolvedModel_->getEmptyWarmStart());
291 basis->resize(prob.nrows_,prob.ncols_);
292 int i;
293 for (i=0;i<prob.ncols_;i++) {
294 CoinWarmStartBasis::Status status =
295 static_cast<CoinWarmStartBasis::Status> (prob.getColumnStatus(i));
296 basis->setStructStatus(i,status);
297 }
298 for (i=0;i<prob.nrows_;i++) {
299 CoinWarmStartBasis::Status status =
300 static_cast<CoinWarmStartBasis::Status> (prob.getRowStatus(i));
301 basis->setArtifStatus(i,status);
302 }
303 presolvedModel_->setWarmStart(basis);
304 delete basis ;
305 delete [] prob.sol_;
306 delete [] prob.acts_;
307 delete [] prob.colstat_;
308 prob.sol_=NULL;
309 prob.acts_=NULL;
310 prob.colstat_=NULL;
311 }
312/*
313 Copy original column and row information from the CoinPresolveMatrix object
314 so it'll be available for postsolve.
315*/
316 int ncolsNow = presolvedModel_->getNumCols();
317 memcpy(originalColumn_,prob.originalColumn_,ncolsNow*sizeof(int));
318 delete [] prob.originalColumn_;
319 prob.originalColumn_=NULL;
320 int nrowsNow = presolvedModel_->getNumRows();
321 memcpy(originalRow_,prob.originalRow_,nrowsNow*sizeof(int));
322 delete [] prob.originalRow_;
323 prob.originalRow_=NULL;
324
325 // now clean up integer variables. This can modify original
326 {
327 int numberChanges=0;
328 const double * lower0 = originalModel_->getColLower();
329 const double * upper0 = originalModel_->getColUpper();
330 const double * lower = presolvedModel_->getColLower();
331 const double * upper = presolvedModel_->getColUpper();
332 for (i=0;i<ncolsNow;i++) {
333 if (!presolvedModel_->isInteger(i))
334 continue;
335 int iOriginal = originalColumn_[i];
336 double lowerValue0 = lower0[iOriginal];
337 double upperValue0 = upper0[iOriginal];
338 double lowerValue = ceil(lower[i]-1.0e-5);
339 double upperValue = floor(upper[i]+1.0e-5);
340 presolvedModel_->setColBounds(i,lowerValue,upperValue);
341 if (lowerValue>upperValue) {
342 numberChanges++;
343 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
344 messages)
345 <<iOriginal
346 <<lowerValue
347 <<upperValue
348 <<CoinMessageEol;
349 result=1;
350 } else {
351 if (lowerValue>lowerValue0+1.0e-8) {
352 originalModel_->setColLower(iOriginal,lowerValue);
353 numberChanges++;
354 }
355 if (upperValue<upperValue0-1.0e-8) {
356 originalModel_->setColUpper(iOriginal,upperValue);
357 numberChanges++;
358 }
359 }
360 }
361 if (numberChanges) {
362 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
363 messages)
364 <<numberChanges
365 <<CoinMessageEol;
366 if (!result&&totalPasses>0&&
367 // we can't go round again in integer if dupcols
368 (prob.presolveOptions_ & 0x80000000) == 0) {
369 result = -1; // round again
370 const CoinPresolveAction *paction = paction_;
371 while (paction) {
372 const CoinPresolveAction *next = paction->next;
373 delete paction;
374 paction = next;
375 }
376 paction_=NULL;
377 }
378 }
379 }
380 } else if (prob.status_ != 0) {
381 // infeasible or unbounded
382 result = 1 ;
383 }
384 }
385 if (!result) {
386 int nrowsAfter = presolvedModel_->getNumRows();
387 int ncolsAfter = presolvedModel_->getNumCols();
388 CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
389 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS, messages)
390 <<nrowsAfter<< -(nrows_ - nrowsAfter)
391 << ncolsAfter<< -(ncols_ - ncolsAfter)
392 <<nelsAfter<< -(nelems_ - nelsAfter)
393 <<CoinMessageEol;
394 } else {
395 gutsOfDestroy();
396 delete presolvedModel_;
397 presolvedModel_=NULL;
398 }
399 return presolvedModel_;
400}
401
402// Return pointer to presolved model
403OsiSolverInterface *
404OsiPresolve::model() const
405{
406 return presolvedModel_;
407}
408// Return pointer to original model
409OsiSolverInterface *
410OsiPresolve::originalModel() const
411{
412 return originalModel_;
413}
414void
415OsiPresolve::postsolve(bool updateStatus)
416{
417 // Messages
418 CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
419 if (!presolvedModel_->isProvenOptimal()) {
420 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
421 messages)
422 <<CoinMessageEol;
423 }
424
425 // this is the size of the original problem
426 const int ncols0 = ncols_;
427 const int nrows0 = nrows_;
428 const CoinBigIndex nelems0 = nelems_;
429
430 // reality check
431 assert(ncols0==originalModel_->getNumCols());
432 assert(nrows0==originalModel_->getNumRows());
433
434 // this is the reduced problem
435 int ncols = presolvedModel_->getNumCols();
436 int nrows = presolvedModel_->getNumRows();
437
438 double *acts = new double [nrows0];
439 double *sol = new double [ncols0];
440 CoinZeroN(acts,nrows0);
441 CoinZeroN(sol,ncols0);
442
443 unsigned char * rowstat=NULL;
444 unsigned char * colstat = NULL;
445 CoinWarmStartBasis * presolvedBasis =
446 dynamic_cast<CoinWarmStartBasis*>(presolvedModel_->getWarmStart());
447 if (!presolvedBasis)
448 updateStatus=false;
449 if (updateStatus) {
450 colstat = new unsigned char[ncols0+nrows0];
451# ifdef ZEROFAULT
452 memset(colstat,0,((ncols0+nrows0)*sizeof(char))) ;
453# endif
454 rowstat = colstat + ncols0;
455 int i;
456 for (i=0;i<ncols;i++) {
457 colstat[i] = presolvedBasis->getStructStatus(i);
458 }
459 for (i=0;i<nrows;i++) {
460 rowstat[i] = presolvedBasis->getArtifStatus(i);
461 }
462 }
463 delete presolvedBasis;
464
465# if PRESOLVE_CONSISTENCY > 0
466 if (updateStatus)
467 { int basicCnt = 0 ;
468 int i ;
469 for (i = 0 ; i < ncols ; i++)
470 { if (colstat[i] == CoinWarmStartBasis::basic) basicCnt++ ; }
471 for (i = 0 ; i < nrows ; i++)
472 { if (rowstat[i] == CoinWarmStartBasis::basic) basicCnt++ ; }
473
474 assert (basicCnt == nrows) ;
475 }
476# endif
477
478/*
479 Postsolve back to the original problem. The CoinPostsolveMatrix object
480 assumes ownership of sol, acts, colstat, and rowstat.
481*/
482 CoinPostsolveMatrix prob(presolvedModel_, ncols0, nrows0, nelems0,
483 presolvedModel_->getObjSense(),
484 sol, acts, colstat, rowstat);
485
486 postsolve(prob);
487
488
489# if PRESOLVE_CONSISTENCY > 0
490 if (updateStatus)
491 { int basicCnt = 0 ;
492 int i ;
493 for (i = 0 ; i < ncols0 ; i++)
494 { if (prob.getColumnStatus(i) == CoinWarmStartBasis::basic) basicCnt++ ; }
495 for (i = 0 ; i < nrows0 ; i++)
496 { if (prob.getRowStatus(i) == CoinWarmStartBasis::basic) basicCnt++ ; }
497
498 assert (basicCnt == nrows0) ;
499 }
500# endif
501
502 originalModel_->setColSolution(sol);
503 if (updateStatus) {
504 CoinWarmStartBasis *basis =
505 dynamic_cast<CoinWarmStartBasis *>(presolvedModel_->getEmptyWarmStart());
506 basis->setSize(ncols0,nrows0);
507 int i;
508 for (i=0;i<ncols0;i++) {
509 CoinWarmStartBasis::Status status = static_cast<CoinWarmStartBasis::Status> (prob.getColumnStatus(i));
510 /* FIXME: these asserts seem correct, but seem to reveal some bugs in CoinPresolve */
511 // assert(status != CoinWarmStartBasis::atLowerBound || originalModel_->getColLower()[i] > -originalModel_->getInfinity());
512 // assert(status != CoinWarmStartBasis::atUpperBound || originalModel_->getColUpper()[i] < originalModel_->getInfinity());
513 basis->setStructStatus(i,status);
514 }
515 for (i=0;i<nrows0;i++) {
516 CoinWarmStartBasis::Status status = static_cast<CoinWarmStartBasis::Status> (prob.getRowStatus(i));
517 /* FIXME: these asserts seem correct, but seem to reveal some bugs in CoinPresolve */
518 // assert(status != CoinWarmStartBasis::atUpperBound || originalModel_->getRowLower()[i] > -originalModel_->getInfinity());
519 // assert(status != CoinWarmStartBasis::atLowerBound || originalModel_->getRowUpper()[i] < originalModel_->getInfinity());
520 basis->setArtifStatus(i,status);
521 }
522 originalModel_->setWarmStart(basis);
523 delete basis ;
524 }
525
526}
527
528// return pointer to original columns
529const int *
530OsiPresolve::originalColumns() const
531{
532 return originalColumn_;
533}
534// return pointer to original rows
535const int *
536OsiPresolve::originalRows() const
537{
538 return originalRow_;
539}
540// Set pointer to original model
541void
542OsiPresolve::setOriginalModel(OsiSolverInterface * model)
543{
544 originalModel_=model;
545}
546
547#if 0
548// A lazy way to restrict which transformations are applied
549// during debugging.
550static int ATOI(const char *name)
551{
552 return true;
553#if PRESOLVE_DEBUG || PRESOLVE_SUMMARY
554 if (getenv(name)) {
555 int val = atoi(getenv(name));
556 printf("%s = %d\n", name, val);
557 return (val);
558 } else {
559 if (strcmp(name,"off"))
560 return (true);
561 else
562 return (false);
563 }
564#else
565 return (true);
566#endif
567}
568#endif
569
570#if PRESOLVE_DEBUG
571// Anonymous namespace for debug routines
572namespace {
573
574/*
575 A control routine for debug checks --- keeps down the clutter in doPresolve.
576 Each time it's called, it prints a list of transforms applied since the last
577 call, then does checks.
578*/
579
580void check_and_tell (const CoinPresolveMatrix *const prob,
581 const CoinPresolveAction *first,
582 const CoinPresolveAction *&mark)
583
584{ const CoinPresolveAction *current ;
585
586 if (first != mark)
587 { printf("PRESOLVE: applied") ;
588 for (current = first ;
589 current != mark && current != 0 ;
590 current = current->next)
591 { printf(" %s",current->name()) ; }
592 printf("\n") ;
593
594 presolve_check_sol(prob) ;
595 presolve_check_nbasic(prob) ;
596 mark = first ; }
597
598 return ; }
599
600} // end anonymous namespace for debug routines
601#endif
602//#define COIN_PRESOLVE_BUG
603#ifdef COIN_PRESOLVE_BUG
604double * debugSolution = NULL;
605int debugNumberColumns = -1;
606static int counter=1000000;
607static bool break2(CoinPresolveMatrix *prob)
608{
609 if (counter>0)
610 printf("counter %d\n",counter);
611 counter--;
612 if (debugSolution&&prob->ncols_==debugNumberColumns) {
613 for (int i=0;i<prob->ncols_;i++) {
614 double value = debugSolution[i];
615 if (value<prob->clo_[i]) {
616 printf("%d inf %g %g %g\n",i,prob->clo_[i],value,prob->cup_[i]);
617 } else if (value>prob->cup_[i]) {
618 printf("%d inf %g %g %g\n",i,prob->clo_[i],value,prob->cup_[i]);
619 }
620 }
621 }
622 if (!counter) {
623 printf("skipping next and all\n");
624 }
625 return (counter<=0);
626}
627#define possibleBreak if (break2(prob)) break
628#define possibleSkip if (!break2(prob))
629#else
630#define possibleBreak
631#define possibleSkip
632#endif
633// This is the presolve loop.
634// It is a separate virtual function so that it can be easily
635// customized by subclassing CoinPresolve.
636
637const CoinPresolveAction *OsiPresolve::presolve(CoinPresolveMatrix *prob)
638{
639 paction_ = 0;
640
641 prob->status_=0; // say feasible
642
643# if PRESOLVE_DEBUG
644 const CoinPresolveAction *pactiond = 0 ;
645 presolve_check_sol(prob) ;
646# endif
647 if ((presolveActions_&4)!=0)
648 transferCosts(prob);
649
650/*
651 Fix variables before we get into the main transform loop.
652*/
653 paction_ = make_fixed(prob, paction_);
654
655# if PRESOLVE_DEBUG
656 check_and_tell(prob,paction_,pactiond) ;
657# endif
658
659 // if integers then switch off dual stuff
660 // later just do individually
661 bool doDualStuff = true;
662 if ((presolveActions_&1)==0) {
663 int i;
664 int ncol = presolvedModel_->getNumCols();
665 for (i=0;i<ncol;i++)
666 if (presolvedModel_->isInteger(i))
667 doDualStuff=false;
668 }
669
670
671# if CHECK_CONSISTENCY
672 presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
673# endif
674
675/*
676 If we're feasible, set up for the main presolve transform loop.
677*/
678 if (!prob->status_) {
679# if 0
680/*
681 This block is used during debugging. See ATOI to see how it works. Some
682 editing will be required to turn it all on.
683*/
684 bool slackd = ATOI("SLACKD")!=0;
685 //bool forcing = ATOI("FORCING")!=0;
686 bool doubleton = ATOI("DOUBLETON")!=0;
687 bool forcing = ATOI("off")!=0;
688 bool ifree = ATOI("off")!=0;
689 bool zerocost = ATOI("off")!=0;
690 bool dupcol = ATOI("off")!=0;
691 bool duprow = ATOI("off")!=0;
692 bool dual = ATOI("off")!=0;
693# else
694# if 1
695 // normal operation --- all transforms enabled
696 bool slackSingleton = true;
697 bool slackd = true;
698 bool doubleton = true;
699 bool tripleton = true;
700 //#define NO_FORCING
701#ifndef NO_FORCING
702 bool forcing = true;
703#endif
704 bool ifree = true;
705 bool zerocost = true;
706 bool dupcol = true;
707 bool duprow = true;
708 bool dual = doDualStuff;
709# else
710 // compile time selection of transforms.
711 bool slackd = false;
712 bool doubleton = true;
713 bool tripleton = true;
714 bool forcing = true;
715 bool ifree = false;
716 bool zerocost = false;
717 bool dupcol = false;
718 bool duprow = false;
719 bool dual = false;
720# endif
721# endif
722 // Switch off some stuff if would annoy set partitioning etc
723 if ((presolveActions_&2)!=0) {
724 doubleton = false;
725 tripleton = false;
726 ifree = false;
727 }
728 // stop x+y+z==1
729 if ((presolveActions_&8)!=0)
730 prob->setPresolveOptions(prob->presolveOptions()|4);
731 // switch on stuff which can't be unrolled easily
732 if ((presolveActions_&16)!=0)
733 prob->setPresolveOptions(prob->presolveOptions()|16);
734 // switch on gub stuff
735 if ((presolveActions_&32)!=0)
736 prob->setPresolveOptions(prob->presolveOptions()|32);
737
738 /*
739 The main loop (just below) starts with a minor loop that does
740 inexpensive presolve transforms until convergence. At each iteration
741 of this loop, next[Rows,Cols]ToDo is copied over to [rows,cols]ToDo.
742
743 Then there's a block like the one here, which sets [rows,cols]ToDo for
744 all rows & cols, followed by executions of a set of expensive
745 transforms. Then we come back around for another iteration of the main
746 loop. [rows,cols]ToDo is not reset as we come back around, so we dive
747 into the inexpensive loop set up to process all.
748 */
749
750 int i;
751 // say look at all
752 if (!prob->anyProhibited()) {
753 for (i=0;i<nrows_;i++)
754 prob->rowsToDo_[i]=i;
755 prob->numberRowsToDo_=nrows_;
756 for (i=0;i<ncols_;i++)
757 prob->colsToDo_[i]=i;
758 prob->numberColsToDo_=ncols_;
759 } else {
760 // some stuff must be left alone
761 prob->numberRowsToDo_=0;
762 for (i=0;i<nrows_;i++)
763 if (!prob->rowProhibited(i))
764 prob->rowsToDo_[prob->numberRowsToDo_++]=i;
765 prob->numberColsToDo_=0;
766 for (i=0;i<ncols_;i++)
767 if (!prob->colProhibited(i))
768 prob->colsToDo_[prob->numberColsToDo_++]=i;
769 }
770
771
772 int iLoop;
773 if (dupcol) {
774 // maybe allow integer columns to be checked
775 if ((presolveActions_&1)!=0)
776 prob->setPresolveOptions(prob->presolveOptions()|1);
777 possibleSkip;
778 paction_ = dupcol_action::presolve(prob, paction_);
779 }
780 if (duprow) {
781 possibleSkip;
782 paction_ = duprow_action::presolve(prob, paction_);
783 }
784 // Check number rows dropped
785 int lastDropped=0;
786 /*
787 Note that pass_ is incremented in testRedundant, evoked from
788 implied_free_action. The bulk of testRedundant is executed every other
789 pass.
790 */
791 prob->pass_=0;
792 for (iLoop=0;iLoop<numberPasses_;iLoop++) {
793
794# ifdef PRESOLVE_SUMMARY
795 printf("Starting major pass %d\n",iLoop+1);
796# endif
797
798 const CoinPresolveAction * const paction0 = paction_;
799 // look for substitutions with no fill
800 //#define IMPLIED 3
801#ifdef IMPLIED
802 int fill_level=3;
803#define IMPLIED2 1
804#if IMPLIED!=3
805#if IMPLIED>0&&IMPLIED<11
806 fill_level=IMPLIED;
807 printf("** fill_level == %d !\n",fill_level);
808#endif
809#if IMPLIED>11&&IMPLIED<21
810 fill_level=-(IMPLIED-10);
811 printf("** fill_level == %d !\n",fill_level);
812#endif
813#endif
814#else
815 int fill_level=2;
816#endif
817 int whichPass=0;
818/*
819 Apply inexpensive transforms until convergence.
820*/
821 while (1) {
822 whichPass++;
823 prob->pass_++;
824 const CoinPresolveAction * const paction1 = paction_;
825
826 if (slackd) {
827 bool notFinished = true;
828 while (notFinished) {
829 possibleBreak;
830 paction_ = slack_doubleton_action::presolve(prob, paction_,
831 notFinished);
832 }
833 if (prob->status_)
834 break;
835# if PRESOLVE_DEBUG
836 check_and_tell(prob,paction_,pactiond) ;
837# endif
838 }
839
840 if (dual&&whichPass==1) {
841 possibleBreak;
842 // this can also make E rows so do one bit here
843 paction_ = remove_dual_action::presolve(prob, paction_);
844 if (prob->status_)
845 break;
846 }
847
848 if (doubleton) {
849 possibleBreak;
850 paction_ = doubleton_action::presolve(prob, paction_);
851 if (prob->status_)
852 break;
853# if PRESOLVE_DEBUG
854 check_and_tell(prob,paction_,pactiond) ;
855# endif
856 }
857
858 if (tripleton) {
859 possibleBreak;
860 paction_ = tripleton_action::presolve(prob, paction_);
861 if (prob->status_)
862 break;
863# if PRESOLVE_DEBUG
864 check_and_tell(prob,paction_,pactiond) ;
865# endif
866 }
867
868 if (zerocost) {
869 possibleBreak;
870 paction_ = do_tighten_action::presolve(prob, paction_);
871 if (prob->status_)
872 break;
873# if PRESOLVE_DEBUG
874 check_and_tell(prob,paction_,pactiond) ;
875# endif
876 }
877
878#ifndef NO_FORCING
879 if (forcing) {
880 possibleBreak;
881 paction_ = forcing_constraint_action::presolve(prob, paction_);
882 if (prob->status_)
883 break;
884# if PRESOLVE_DEBUG
885 check_and_tell(prob,paction_,pactiond) ;
886# endif
887 }
888#endif
889
890 if (ifree&&(whichPass%5)==1) {
891 possibleBreak;
892 paction_ = implied_free_action::presolve(prob, paction_,fill_level);
893 if (prob->status_)
894 break;
895# if PRESOLVE_DEBUG
896 check_and_tell(prob,paction_,pactiond) ;
897# endif
898 }
899
900# if CHECK_CONSISTENCY
901 presolve_links_ok(prob->rlink_,prob->mrstrt_,
902 prob->hinrow_,prob->nrows_) ;
903# endif
904
905# if 0 && PRESOLVE_DEBUG
906
907 /*
908 For reasons that escape me just now, the linker is unable to find
909 this function. Copying the code from CoinPresolvePsdebug to the head
910 of this routine works just fine. Library loading order looks ok. Other
911 routines from CoinPresolvePsdebug are found. I'm stumped. -- lh --
912 */
913
914 presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
915 prob->ncols_);
916# endif
917# if CHECK_CONSISTENCY
918 prob->consistent();
919# endif
920
921
922 // set up for next pass
923 // later do faster if many changes i.e. memset and memcpy
924 prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
925 int kcheck; // debug?
926 bool found=false;
927 kcheck=-1;
928 for (i=0;i<prob->numberNextRowsToDo_;i++) {
929 int index = prob->nextRowsToDo_[i];
930 prob->unsetRowChanged(index);
931 prob->rowsToDo_[i] = index;
932 if (index==kcheck) {
933 printf("row %d on list after pass %d\n",kcheck,
934 whichPass);
935 found=true;
936 }
937 }
938 if (!found&&kcheck>=0)
939 prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
940 prob->numberNextRowsToDo_=0;
941 prob->numberColsToDo_ = prob->numberNextColsToDo_;
942 kcheck=-1;
943 found=false;
944 for (i=0;i<prob->numberNextColsToDo_;i++) {
945 int index = prob->nextColsToDo_[i];
946 prob->unsetColChanged(index);
947 prob->colsToDo_[i] = index;
948 if (index==kcheck) {
949 printf("col %d on list after pass %d\n",kcheck,
950 whichPass);
951 found=true;
952 }
953 }
954 if (!found&&kcheck>=0)
955 prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
956 prob->numberNextColsToDo_=0;
957 if (paction_ == paction1&&fill_level>0)
958 break;
959 } // End of inexpensive transform loop
960
961 // say look at all
962 int i;
963 if (!prob->anyProhibited()) {
964 for (i=0;i<nrows_;i++)
965 prob->rowsToDo_[i]=i;
966 prob->numberRowsToDo_=nrows_;
967 for (i=0;i<ncols_;i++)
968 prob->colsToDo_[i]=i;
969 prob->numberColsToDo_=ncols_;
970 } else {
971 // some stuff must be left alone
972 prob->numberRowsToDo_=0;
973 for (i=0;i<nrows_;i++)
974 if (!prob->rowProhibited(i))
975 prob->rowsToDo_[prob->numberRowsToDo_++]=i;
976 prob->numberColsToDo_=0;
977 for (i=0;i<ncols_;i++)
978 if (!prob->colProhibited(i))
979 prob->colsToDo_[prob->numberColsToDo_++]=i;
980 }
981 // now expensive things
982 // this caused world.mps to run into numerical difficulties
983
984# ifdef PRESOLVE_SUMMARY
985 printf("Starting expensive\n");
986# endif
987
988 if (dual) {
989 int itry;
990 for (itry=0;itry<5;itry++) {
991 const CoinPresolveAction * const paction2 = paction_;
992 possibleBreak;
993 paction_ = remove_dual_action::presolve(prob, paction_);
994# if PRESOLVE_DEBUG
995 check_and_tell(prob,paction_,pactiond) ;
996# endif
997 if (prob->status_)
998 break;
999 if (ifree) {
1000#ifdef IMPLIED
1001#if IMPLIED2 ==0
1002 int fill_level=0; // switches off substitution
1003#elif IMPLIED2!=99
1004 int fill_level=IMPLIED2;
1005#endif
1006#endif
1007 if ((itry&1)==0) {
1008 possibleBreak;
1009 paction_ = implied_free_action::presolve(prob, paction_,fill_level);
1010 }
1011# if PRESOLVE_DEBUG
1012 check_and_tell(prob,paction_,pactiond) ;
1013# endif
1014 if (prob->status_)
1015 break;
1016 }
1017 if (paction_ == paction2)
1018 break;
1019 }
1020 } else if (ifree) {
1021 // just free
1022#ifdef IMPLIED
1023#if IMPLIED2 ==0
1024 int fill_level=0; // switches off substitution
1025#elif IMPLIED2!=99
1026 int fill_level=IMPLIED2;
1027#endif
1028#endif
1029 possibleBreak;
1030 paction_ = implied_free_action::presolve(prob, paction_,fill_level);
1031 if (prob->status_)
1032 break;
1033 }
1034
1035 if (dupcol) {
1036 // maybe allow integer columns to be checked
1037 if ((presolveActions_&1)!=0)
1038 prob->setPresolveOptions(prob->presolveOptions()|1);
1039 possibleBreak;
1040 paction_ = dupcol_action::presolve(prob, paction_);
1041# if PRESOLVE_DEBUG
1042 check_and_tell(prob,paction_,pactiond) ;
1043# endif
1044 if (prob->status_)
1045 break;
1046 }
1047
1048 if (duprow) {
1049 possibleBreak;
1050 paction_ = duprow_action::presolve(prob, paction_);
1051# if PRESOLVE_DEBUG
1052 check_and_tell(prob,paction_,pactiond) ;
1053# endif
1054 if (prob->status_)
1055 break;
1056 }
1057 if ((presolveActions_&32)!=0) {
1058 possibleBreak;
1059 paction_ = gubrow_action::presolve(prob, paction_);
1060 }
1061
1062 bool stopLoop=false;
1063 {
1064 int * hinrow = prob->hinrow_;
1065 int numberDropped=0;
1066 for (i=0;i<nrows_;i++)
1067 if (!hinrow[i])
1068 numberDropped++;
1069 //printf("%d rows dropped after pass %d\n",numberDropped,
1070 // iLoop+1);
1071 if (numberDropped==lastDropped)
1072 stopLoop=true;
1073 else
1074 lastDropped = numberDropped;
1075 }
1076 // Do this here as not very loopy
1077 if (slackSingleton) {
1078 // On most passes do not touch costed slacks
1079 if (paction_ != paction0&&!stopLoop) {
1080 possibleBreak;
1081 paction_ = slack_singleton_action::presolve(prob, paction_,NULL);
1082 } else {
1083 // do costed if Clp (at end as ruins rest of presolve)
1084 possibleBreak;
1085 paction_ = slack_singleton_action::presolve(prob, paction_,NULL);
1086 stopLoop=true;
1087 }
1088 }
1089#if PRESOLVE_DEBUG
1090 presolve_check_sol(prob,1);
1091#endif
1092 if (paction_ == paction0||stopLoop)
1093 break;
1094
1095 } // End of major pass loop
1096 }
1097/*
1098 Final cleanup: drop zero coefficients from the matrix, then drop empty rows
1099 and columns.
1100*/
1101 if (!prob->status_) {
1102 paction_ = drop_zero_coefficients(prob, paction_);
1103# if PRESOLVE_DEBUG
1104 check_and_tell(prob,paction_,pactiond) ;
1105# endif
1106
1107 paction_ = drop_empty_cols_action::presolve(prob, paction_);
1108# if PRESOLVE_DEBUG
1109 check_and_tell(prob,paction_,pactiond) ;
1110# endif
1111
1112 paction_ = drop_empty_rows_action::presolve(prob, paction_);
1113# if PRESOLVE_DEBUG
1114 check_and_tell(prob,paction_,pactiond) ;
1115# endif
1116 }
1117 // Messages
1118 CoinMessages messages = CoinMessage(prob->messages().language());
1119 if (prob->status_) {
1120 if (prob->status_==1)
1121 prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
1122 messages)
1123 <<prob->feasibilityTolerance_
1124 <<CoinMessageEol;
1125 else if (prob->status_==2)
1126 prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
1127 messages)
1128 <<CoinMessageEol;
1129 else
1130 prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
1131 messages)
1132 <<CoinMessageEol;
1133 // get rid of data
1134 gutsOfDestroy();
1135 }
1136 return (paction_);
1137}
1138
1139
1140// We could have implemented this by having each postsolve routine
1141// directly call the next one, but this may make it easier to add debugging checks.
1142void OsiPresolve::postsolve(CoinPostsolveMatrix &prob)
1143{
1144 const CoinPresolveAction *paction = paction_;
1145
1146#if PRESOLVE_DEBUG
1147 printf("Begin POSTSOLVING\n") ;
1148 if (prob.colstat_)
1149 { presolve_check_nbasic(&prob);
1150 presolve_check_sol(&prob); }
1151 presolve_check_duals(&prob);
1152#endif
1153
1154
1155 while (paction) {
1156# if PRESOLVE_DEBUG
1157 printf("POSTSOLVING %s\n", paction->name());
1158# endif
1159
1160 paction->postsolve(&prob);
1161
1162# if PRESOLVE_DEBUG
1163 if (prob.colstat_)
1164 { presolve_check_nbasic(&prob);
1165 presolve_check_sol(&prob); }
1166# endif
1167 paction = paction->next;
1168# if PRESOLVE_DEBUG
1169 presolve_check_duals(&prob);
1170# endif
1171 }
1172# if PRESOLVE_DEBUG
1173 printf("End POSTSOLVING\n") ;
1174# endif
1175
1176#if 0 && PRESOLVE_DEBUG
1177
1178 << This block of checks will require some work to get it to compile. >>
1179
1180 for (i=0; i<ncols0; i++) {
1181 if (!cdone_[i]) {
1182 printf("!cdone_[%d]\n", i);
1183 abort();
1184 }
1185 }
1186
1187 for (int i=0; i<nrows0; i++) {
1188 if (!rdone[i]) {
1189 printf("!rdone[%d]\n", i);
1190 abort();
1191 }
1192 }
1193
1194
1195 for (i=0; i<ncols0; i++) {
1196 if (sol[i] < -1e10 || sol[i] > 1e10)
1197 printf("!!!%d %g\n", i, sol[i]);
1198
1199 }
1200
1201
1202#endif
1203
1204#if 0 && PRESOLVE_DEBUG
1205
1206 << This block of checks will require some work to get it to compile. >>
1207
1208 // debug check: make sure we ended up with same original matrix
1209 {
1210 int identical = 1;
1211
1212 for (int i=0; i<ncols0; i++) {
1213 PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
1214 CoinBigIndex kcs0 = &prob->mcstrt0[i];
1215 CoinBigIndex kcs = mcstrt[i];
1216 int n = hincol[i];
1217 for (int k=0; k<n; k++) {
1218 CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
1219
1220 if (k1 == kcs+n) {
1221 printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
1222 abort();
1223 }
1224
1225 if (colels[k1] != &prob->dels0[kcs0+k])
1226 printf("BAD COLEL[%d %d %d]: %g\n",
1227 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
1228
1229 if (kcs0+k != k1)
1230 identical=0;
1231 }
1232 }
1233 printf("identical? %d\n", identical);
1234 }
1235#endif
1236 // put back duals
1237 double maxmin = originalModel_->getObjSense();
1238 if (maxmin<0.0) {
1239 // swap signs
1240 int i;
1241 double * pi = prob.rowduals_;
1242 for (i=0;i<nrows_;i++)
1243 pi[i] = -pi[i];
1244 }
1245 originalModel_->setRowPrice(prob.rowduals_);
1246 // Now check solution
1247 // **** code later - has to be by hand
1248#if 0
1249 memcpy(originalModel_->dualColumnSolution(),
1250 originalModel_->objective(),ncols_*sizeof(double));
1251 originalModel_->transposeTimes(-1.0,
1252 originalModel_->dualRowSolution(),
1253 originalModel_->dualColumnSolution());
1254 memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
1255 originalModel_->times(1.0,originalModel_->primalColumnSolution(),
1256 originalModel_->primalRowSolution());
1257 originalModel_->checkSolution();
1258 // Messages
1259 CoinMessages messages = CoinMessage(presolvedModel_->messages().language());
1260 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
1261 messages)
1262 <<originalModel_->objectiveValue()
1263 <<originalModel_->sumDualInfeasibilities()
1264 <<originalModel_->numberDualInfeasibilities()
1265 <<originalModel_->sumPrimalInfeasibilities()
1266 <<originalModel_->numberPrimalInfeasibilities()
1267 <<CoinMessageEol;
1268 //originalModel_->objectiveValue_=objectiveValue_;
1269 originalModel_->setNumberIterations(presolvedModel_->numberIterations());
1270 if (!presolvedModel_->status()) {
1271 if (!originalModel_->numberDualInfeasibilities()&&
1272 !originalModel_->numberPrimalInfeasibilities()) {
1273 originalModel_->setProblemStatus( 0);
1274 } else {
1275 originalModel_->setProblemStatus( -1);
1276 presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
1277 messages)
1278 <<CoinMessageEol;
1279 }
1280 } else {
1281 originalModel_->setProblemStatus( presolvedModel_->status());
1282 }
1283#endif
1284}
1285
1286
1287static inline double getTolerance(const OsiSolverInterface *si, OsiDblParam key)
1288{
1289 double tol;
1290 if (! si->getDblParam(key, tol)) {
1291 CoinPresolveAction::throwCoinError("getDblParam failed",
1292 "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
1293 }
1294 return (tol);
1295}
1296
1297
1298// Assumptions:
1299// 1. nrows>=m.getNumRows()
1300// 2. ncols>=m.getNumCols()
1301//
1302// In presolve, these values are equal.
1303// In postsolve, they may be inequal, since the reduced problem
1304// may be smaller, but we need room for the large problem.
1305// ncols may be larger than si.getNumCols() in postsolve,
1306// this at that point si will be the reduced problem,
1307// but we need to reserve enough space for the original problem.
1308
1309CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const OsiSolverInterface * si,
1310 int ncols_in,
1311 int nrows_in,
1312 CoinBigIndex nelems_in) :
1313 ncols_(si->getNumCols()),
1314 nelems_(si->getNumElements()),
1315 ncols0_(ncols_in),
1316 nrows0_(nrows_in),
1317 bulkRatio_(2.0),
1318
1319 mcstrt_(new CoinBigIndex[ncols_in+1]),
1320 hincol_(new int[ncols_in+1]),
1321
1322 cost_(new double[ncols_in]),
1323 clo_(new double[ncols_in]),
1324 cup_(new double[ncols_in]),
1325 rlo_(new double[nrows_in]),
1326 rup_(new double[nrows_in]),
1327 originalColumn_(new int[ncols_in]),
1328 originalRow_(new int[nrows_in]),
1329
1330 ztolzb_(getTolerance(si, OsiPrimalTolerance)),
1331 ztoldj_(getTolerance(si, OsiDualTolerance)),
1332
1333 maxmin_(si->getObjSense()),
1334
1335 handler_(0),
1336 defaultHandler_(false),
1337 messages_()
1338
1339{
1340 bulk0_ = static_cast<CoinBigIndex>(bulkRatio_*nelems_in) ;
1341 hrow_ = new int [bulk0_] ;
1342 colels_ = new double[bulk0_] ;
1343
1344 si->getDblParam(OsiObjOffset,originalOffset_);
1345 int ncols = si->getNumCols();
1346 int nrows = si->getNumRows();
1347
1348 setMessageHandler(si->messageHandler()) ;
1349
1350 CoinDisjointCopyN(si->getColLower(), ncols, clo_);
1351 CoinDisjointCopyN(si->getColUpper(), ncols, cup_);
1352 CoinDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
1353 CoinDisjointCopyN(si->getRowLower(), nrows, rlo_);
1354 CoinDisjointCopyN(si->getRowUpper(), nrows, rup_);
1355 int i;
1356 // initialize and clean up bounds
1357 double infinity = si->getInfinity();
1358 if (infinity!=COIN_DBL_MAX) {
1359 for (i=0;i<ncols;i++) {
1360 if (clo_[i]==-infinity)
1361 clo_[i]=-COIN_DBL_MAX;
1362 if (cup_[i]==infinity)
1363 cup_[i]=COIN_DBL_MAX;
1364 }
1365 for (i=0;i<nrows;i++) {
1366 if (rlo_[i]==-infinity)
1367 rlo_[i]=-COIN_DBL_MAX;
1368 if (rup_[i]==infinity)
1369 rup_[i]=COIN_DBL_MAX;
1370 }
1371 }
1372 for (i=0;i<ncols_in;i++)
1373 originalColumn_[i]=i;
1374 for (i=0;i<nrows_in;i++)
1375 originalRow_[i]=i;
1376 sol_=NULL;
1377 rowduals_=NULL;
1378 acts_=NULL;
1379
1380 rcosts_=NULL;
1381 colstat_=NULL;
1382 rowstat_=NULL;
1383}
1384
1385// I am not familiar enough with CoinPackedMatrix to be confident
1386// that I will implement a row-ordered version of toColumnOrderedGapFree
1387// properly.
1388static bool isGapFree(const CoinPackedMatrix& matrix)
1389{
1390 const CoinBigIndex * start = matrix.getVectorStarts();
1391 const int * length = matrix.getVectorLengths();
1392 int i;
1393 for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
1394 if (start[i+1] - start[i] != length[i])
1395 break;
1396 }
1397 return (! (i >= 0));
1398}
1399CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
1400 double /*maxmin_*/,
1401 // end prepost members
1402 OsiSolverInterface * si,
1403 // rowrep
1404 int nrows_in,
1405 CoinBigIndex nelems_in,
1406 bool doStatus,
1407 double nonLinearValue,
1408 const char * prohibited,
1409 const char * rowProhibited) :
1410
1411 CoinPrePostsolveMatrix(si,
1412 ncols0_in, nrows_in, nelems_in),
1413 clink_(new presolvehlink[ncols0_in+1]),
1414 rlink_(new presolvehlink[nrows_in+1]),
1415
1416 dobias_(0.0),
1417
1418 // temporary init
1419 mrstrt_(new CoinBigIndex[nrows_in+1]),
1420 hinrow_(new int[nrows_in+1]),
1421 integerType_(new unsigned char[ncols0_in]),
1422 tuning_(false),
1423 startTime_(0.0),
1424 feasibilityTolerance_(0.0),
1425 status_(-1),
1426 maxSubstLevel_(3),
1427 colsToDo_(new int [ncols0_in]),
1428 numberColsToDo_(0),
1429 nextColsToDo_(new int[ncols0_in]),
1430 numberNextColsToDo_(0),
1431 rowsToDo_(new int [nrows_in]),
1432 numberRowsToDo_(0),
1433 nextRowsToDo_(new int[nrows_in]),
1434 numberNextRowsToDo_(0),
1435 presolveOptions_(0)
1436{
1437
1438 rowels_ = new double [bulk0_] ;
1439 hcol_ = new int [bulk0_] ;
1440
1441 nrows_ = si->getNumRows() ;
1442 const CoinBigIndex bufsize = static_cast<CoinBigIndex>(bulkRatio_*nelems_in) ;
1443
1444 // Set up change bits
1445 rowChanged_ = new unsigned char[nrows_];
1446 memset(rowChanged_,0,nrows_);
1447 colChanged_ = new unsigned char[ncols_];
1448 memset(colChanged_,0,ncols_);
1449 const CoinPackedMatrix * m1 = si->getMatrixByCol();
1450
1451 // The coefficient matrix is a big hunk of stuff.
1452 // Do the copy here to try to avoid running out of memory.
1453
1454 const CoinBigIndex * start = m1->getVectorStarts();
1455 const int * length = m1->getVectorLengths();
1456 const int * row = m1->getIndices();
1457 const double * element = m1->getElements();
1458 int icol,nel=0;
1459 mcstrt_[0]=0;
1460 for (icol=0;icol<ncols_;icol++) {
1461 int j;
1462 for (j=start[icol];j<start[icol]+length[icol];j++) {
1463 if (fabs(element[j])>ZTOLDP) {
1464 hrow_[nel]=row[j];
1465 colels_[nel++]=element[j];
1466 }
1467 }
1468 hincol_[icol]=nel-mcstrt_[icol];
1469 mcstrt_[icol+1]=nel;
1470 }
1471
1472 // same thing for row rep
1473 CoinPackedMatrix * m = new CoinPackedMatrix();
1474 m->reverseOrderedCopyOf(*si->getMatrixByCol());
1475 // do by hand because of zeros m->removeGaps();
1476 CoinDisjointCopyN(m->getVectorStarts(), nrows_, mrstrt_);
1477 mrstrt_[nrows_] = nelems_;
1478 CoinDisjointCopyN(m->getVectorLengths(), nrows_, hinrow_);
1479 CoinDisjointCopyN(m->getIndices(), nelems_, hcol_);
1480 CoinDisjointCopyN(m->getElements(), nelems_, rowels_);
1481 start = m->getVectorStarts();
1482 length = m->getVectorLengths();
1483 const int * column = m->getIndices();
1484 element = m->getElements();
1485 // out zeros
1486 int irow;
1487 nel=0;
1488 mrstrt_[0]=0;
1489 for (irow=0;irow<nrows_;irow++) {
1490 int j;
1491 for (j=start[irow];j<start[irow]+length[irow];j++) {
1492 if (fabs(element[j])>ZTOLDP) {
1493 hcol_[nel]=column[j];
1494 rowels_[nel++]=element[j];
1495 }
1496 }
1497 hinrow_[irow]=nel-mrstrt_[irow];
1498 mrstrt_[irow+1]=nel;
1499 }
1500 nelems_=nel;
1501
1502 delete m;
1503 {
1504 int i;
1505 for (i=0;i<ncols_;i++) {
1506 if (si->isInteger(i))
1507 integerType_[i] = 1;
1508 else
1509 integerType_[i] = 0;
1510 }
1511 }
1512
1513 // Set up prohibited bits if needed
1514 if (nonLinearValue) {
1515 anyProhibited_ = true;
1516 for (icol=0;icol<ncols_;icol++) {
1517 int j;
1518 bool nonLinearColumn = false;
1519 if (cost_[icol]==nonLinearValue)
1520 nonLinearColumn=true;
1521 for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
1522 if (colels_[j]==nonLinearValue) {
1523 nonLinearColumn=true;
1524 setRowProhibited(hrow_[j]);
1525 }
1526 }
1527 if (nonLinearColumn)
1528 setColProhibited(icol);
1529 }
1530 } else if (prohibited) {
1531 anyProhibited_ = true;
1532 for (icol=0;icol<ncols_;icol++) {
1533 if (prohibited[icol])
1534 setColProhibited(icol);
1535 }
1536 } else {
1537 anyProhibited_ = false;
1538 }
1539 // Any rows special?
1540 if (rowProhibited) {
1541 anyProhibited_ = true;
1542 for (int irow=0;irow<nrows_;irow++) {
1543 if (rowProhibited[irow])
1544 setRowProhibited(irow);
1545 }
1546 }
1547 if (doStatus) {
1548 // allow for status and solution
1549 sol_ = new double[ncols_];
1550 const double *presol ;
1551 presol = si->getColSolution() ;
1552 memcpy(sol_,presol,ncols_*sizeof(double));
1553 acts_ = new double [nrows_];
1554 memcpy(acts_,si->getRowActivity(),nrows_*sizeof(double));
1555 CoinWarmStartBasis * basis =
1556 dynamic_cast<CoinWarmStartBasis*>(si->getWarmStart());
1557 colstat_ = new unsigned char [nrows_+ncols_];
1558 rowstat_ = colstat_+ncols_;
1559 // If basis is NULL then put in all slack basis
1560 if (basis&&basis->getNumStructural()==ncols_) {
1561 int i;
1562 for (i=0;i<ncols_;i++) {
1563 colstat_[i] = basis->getStructStatus(i);
1564 }
1565 for (i=0;i<nrows_;i++) {
1566 rowstat_[i] = basis->getArtifStatus(i);
1567 }
1568 } else {
1569 int i;
1570 // no basis
1571 for (i=0;i<ncols_;i++) {
1572 colstat_[i] = 3;
1573 }
1574 for (i=0;i<nrows_;i++) {
1575 rowstat_[i] = 1;
1576 }
1577 }
1578 delete basis;
1579 }
1580
1581#if 0
1582 for (i=0; i<nrows; ++i)
1583 printf("NR: %6d\n", hinrow[i]);
1584 for (int i=0; i<ncols; ++i)
1585 printf("NC: %6d\n", hincol[i]);
1586#endif
1587
1588#if 0 /* for building against CoinUtils 2.6, this #if 1 need to be changed into an #if 0 */
1589 presolve_make_memlists(mcstrt_, hincol_, clink_, ncols_);
1590 presolve_make_memlists(mrstrt_, hinrow_, rlink_, nrows_);
1591#else
1592 presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
1593 presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
1594#endif
1595
1596 // this allows last col/row to expand up to bufsize-1 (22);
1597 // this must come after the calls to presolve_prefix
1598 mcstrt_[ncols_] = bufsize-1;
1599 mrstrt_[nrows_] = bufsize-1;
1600 // Allocate useful arrays
1601 initializeStuff();
1602
1603#if CHECK_CONSISTENCY
1604 consistent(false);
1605#endif
1606}
1607
1608void CoinPresolveMatrix::update_model(OsiSolverInterface * si,
1609 int /*nrows0*/,
1610 int /*ncols0*/,
1611 CoinBigIndex /*nelems0*/)
1612{
1613 int nels=0;
1614 int i;
1615 for ( i=0; i<ncols_; i++)
1616 nels += hincol_[i];
1617 CoinPackedMatrix m(true,nrows_,ncols_,nels, colels_, hrow_,mcstrt_,hincol_);
1618 si->loadProblem(m, clo_, cup_, cost_, rlo_, rup_);
1619
1620 for ( i=0; i<ncols_; i++) {
1621 if (integerType_[i])
1622 si->setInteger(i);
1623 else
1624 si->setContinuous(i);
1625 }
1626
1627#if PRESOLVE_SUMMARY
1628 printf("NEW NCOL/NROW/NELS: %d(-%d) %d(-%d) %d(-%d)\n",
1629 ncols_, ncols0-ncols_,
1630 nrows_, nrows0-nrows_,
1631 si->getNumElements(), nelems0-si->getNumElements());
1632#endif
1633 si->setDblParam(OsiObjOffset,originalOffset_-dobias_);
1634
1635}
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647//////////////// POSTSOLVE
1648
1649CoinPostsolveMatrix::CoinPostsolveMatrix(OsiSolverInterface* si,
1650 int ncols0_in,
1651 int nrows0_in,
1652 CoinBigIndex nelems0,
1653
1654 double maxmin,
1655 // end prepost members
1656
1657 double *sol_in,
1658 double *acts_in,
1659
1660 unsigned char *colstat_in,
1661 unsigned char *rowstat_in) :
1662
1663 CoinPrePostsolveMatrix(si, ncols0_in, nrows0_in, nelems0),
1664/*
1665 Used only to mark processed columns and rows so that debugging routines know
1666 what to check.
1667*/
1668# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1669 cdone_(new char[ncols0_in]),
1670 rdone_(new char[nrows0_in])
1671# else
1672 cdone_(0),
1673 rdone_(0)
1674# endif
1675
1676{
1677/*
1678 The CoinPrePostsolveMatrix constructor will set bulk0_ to bulkRatio_*nelems0.
1679 By default, bulkRatio_ is 2. This is certainly larger than absolutely
1680 necessary, but good for efficiency (minimises the need to compress the bulk
1681 store). The main storage arrays for the threaded column-major representation
1682 (hrow_, colels_, link_) should be allocated to this size.
1683*/
1684 free_list_ = 0 ;
1685 maxlink_ = bulk0_ ;
1686 link_ = new int[maxlink_] ;
1687
1688 nrows_ = si->getNumRows() ;
1689 ncols_ = si->getNumCols() ;
1690
1691 sol_=sol_in;
1692 rowduals_=NULL;
1693 acts_=acts_in;
1694
1695 rcosts_=NULL;
1696 colstat_=colstat_in;
1697 rowstat_=rowstat_in;
1698
1699 // this is the *reduced* model, which is probably smaller
1700 int ncols1 = ncols_ ;
1701 int nrows1 = nrows_ ;
1702
1703 const CoinPackedMatrix * m = si->getMatrixByCol();
1704#if 0
1705 if (! isGapFree(*m)) {
1706 CoinPresolveAction::throwCoinError("Matrix not gap free",
1707 "CoinPostsolveMatrix");
1708 }
1709#endif
1710 const CoinBigIndex nelemsr = m->getNumElements();
1711
1712 if (isGapFree(*m)) {
1713 CoinDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
1714 CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1715 mcstrt_[ncols_] = nelems0; // points to end of bulk store
1716 CoinDisjointCopyN(m->getVectorLengths(),ncols1, hincol_);
1717 CoinDisjointCopyN(m->getIndices(), nelemsr, hrow_);
1718 CoinDisjointCopyN(m->getElements(), nelemsr, colels_);
1719 }
1720 else
1721 {
1722 CoinPackedMatrix* mm = new CoinPackedMatrix(*m);
1723 if( mm->hasGaps())
1724 mm->removeGaps();
1725 assert(nelemsr == mm->getNumElements());
1726 CoinDisjointCopyN(mm->getVectorStarts(), ncols1, mcstrt_);
1727 CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
1728 mcstrt_[ncols_] = nelems0; // points to end of bulk store
1729 CoinDisjointCopyN(mm->getVectorLengths(),ncols1, hincol_);
1730 CoinDisjointCopyN(mm->getIndices(), nelemsr, hrow_);
1731 CoinDisjointCopyN(mm->getElements(), nelemsr, colels_);
1732 }
1733
1734# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1735 memset(cdone_, -1, ncols0_);
1736 memset(rdone_, -1, nrows0_);
1737# endif
1738
1739 rowduals_ = new double[nrows0_];
1740 CoinDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
1741 rcosts_ = new double[ncols0_];
1742 CoinDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
1743
1744#if PRESOLVE_DEBUG
1745 // check accuracy of reduced costs (rcosts_ is recalculated reduced costs)
1746 si->getMatrixByCol()->transposeTimes(rowduals_,rcosts_);
1747 const double * obj =si->getObjCoefficients();
1748 const double * dj =si->getReducedCost();
1749 {
1750 int i;
1751 for (i=0;i<ncols1;i++) {
1752 double newDj = obj[i]-rcosts_[i];
1753 rcosts_[i]=newDj;
1754 assert (fabs(newDj-dj[i])<1.0e-1);
1755 }
1756 }
1757 // check reduced costs are 0 for basic variables
1758 {
1759 int i;
1760 for (i=0;i<ncols1;i++)
1761 if (columnIsBasic(i))
1762 assert (fabs(rcosts_[i])<1.0e-5);
1763 for (i=0;i<nrows1;i++)
1764 if (rowIsBasic(i))
1765 assert (fabs(rowduals_[i])<1.0e-5);
1766 }
1767#endif
1768
1769 if (maxmin<0.0) {
1770 // change so will look as if minimize
1771 int i;
1772 for (i=0;i<nrows1;i++)
1773 rowduals_[i] = - rowduals_[i];
1774 for (i=0;i<ncols1;i++) {
1775 rcosts_[i] = - rcosts_[i];
1776 }
1777 }
1778
1779/*
1780 CoinPresolve requires both column solution and row activity for correct
1781 operation.
1782*/
1783 CoinDisjointCopyN(si->getColSolution(), ncols1, sol_);
1784 CoinDisjointCopyN(si->getRowActivity(), nrows1, acts_) ;
1785 si->setDblParam(OsiObjOffset,originalOffset_);
1786
1787 for (int j=0; j<ncols1; j++) {
1788 CoinBigIndex kcs = mcstrt_[j];
1789 CoinBigIndex kce = kcs + hincol_[j];
1790 for (CoinBigIndex k=kcs; k<kce; ++k) {
1791 link_[k] = k+1;
1792 }
1793 if (kce>0)
1794 link_[kce-1] = NO_LINK ;
1795 }
1796 if (maxlink_>0) {
1797 int ml = maxlink_;
1798 for (CoinBigIndex k=nelemsr; k<ml; ++k)
1799 link_[k] = k+1;
1800 link_[ml-1] = NO_LINK;
1801 }
1802 free_list_ = nelemsr;
1803
1804# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
1805/*
1806 These are used to track the action of postsolve transforms during debugging.
1807*/
1808 CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
1809 CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
1810 CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
1811 CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
1812# endif
1813}
1814