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

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

bug on tiny problems

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