source: trunk/Clp/src/ClpPresolve.cpp @ 1780

Last change on this file since 1780 was 1780, checked in by forrest, 8 years ago

stuff for parametrics

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