source: stable/1.16/Clp/src/ClpPresolve.cpp @ 2357

Last change on this file since 2357 was 2357, checked in by forrest, 6 months ago

bug in tiny problems

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