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

Last change on this file since 1721 was 1721, checked in by forrest, 9 years ago

a bit of debug and try and save memory in OsiClp?

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