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

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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