/* $Id$ */
// Copyright (C) 2002, International Business Machines
// Corporation and others.  All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).

//#define	PRESOLVE_CONSISTENCY	1
//#define	PRESOLVE_DEBUG	1

#include <stdio.h>

#include <cassert>
#include <iostream>

#include "CoinHelperFunctions.hpp"
#include "ClpConfig.h"
#ifdef CLP_HAS_ABC
#include "CoinAbcCommon.hpp"
#endif

#include "CoinPackedMatrix.hpp"
#include "ClpPackedMatrix.hpp"
#include "ClpSimplex.hpp"
#include "ClpSimplexOther.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif

#include "ClpPresolve.hpp"
#include "CoinPresolveMatrix.hpp"

#include "CoinPresolveEmpty.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolvePsdebug.hpp"
#include "CoinPresolveSingleton.hpp"
#include "CoinPresolveDoubleton.hpp"
#include "CoinPresolveTripleton.hpp"
#include "CoinPresolveZeros.hpp"
#include "CoinPresolveSubst.hpp"
#include "CoinPresolveForcing.hpp"
#include "CoinPresolveDual.hpp"
#include "CoinPresolveTighten.hpp"
#include "CoinPresolveUseless.hpp"
#include "CoinPresolveDupcol.hpp"
#include "CoinPresolveImpliedFree.hpp"
#include "CoinPresolveIsolated.hpp"
#include "CoinMessage.hpp"

ClpPresolve::ClpPresolve()
  : originalModel_(NULL)
  , presolvedModel_(NULL)
  , nonLinearValue_(0.0)
  , originalColumn_(NULL)
  , originalRow_(NULL)
  , rowObjective_(NULL)
  , paction_(0)
  , ncols_(0)
  , nrows_(0)
  , nelems_(0)
  ,
#ifdef CLP_INHERIT_MODE
  numberPasses_(20)
  ,
#else
  numberPasses_(5)
  ,
#endif
  substitution_(3)
  ,
#ifndef CLP_NO_STD
  saveFile_("")
  ,
#endif
  presolveActions_(0)
{
}

ClpPresolve::~ClpPresolve()
{
  destroyPresolve();
}
// Gets rid of presolve actions (e.g.when infeasible)
void ClpPresolve::destroyPresolve()
{
  const CoinPresolveAction *paction = paction_;
  while (paction) {
    const CoinPresolveAction *next = paction->next;
    delete paction;
    paction = next;
  }
  delete[] originalColumn_;
  delete[] originalRow_;
  paction_ = NULL;
  originalColumn_ = NULL;
  originalRow_ = NULL;
  delete[] rowObjective_;
  rowObjective_ = NULL;
}

/* This version of presolve returns a pointer to a new presolved
   model.  NULL if infeasible
*/
ClpSimplex *
ClpPresolve::presolvedModel(ClpSimplex &si,
  double feasibilityTolerance,
  bool keepIntegers,
  int numberPasses,
  bool dropNames,
  bool doRowObjective,
  const char *prohibitedRows,
  const char *prohibitedColumns)
{
  // Check matrix
  int checkType = ((si.specialOptions() & 128) != 0) ? 14 : 15;
  if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
        1.0e20, checkType))
    return NULL;
  else
    return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
      doRowObjective,
      prohibitedRows,
      prohibitedColumns);
}
#ifndef CLP_NO_STD
/* This version of presolve updates
   model and saves original data to file.  Returns non-zero if infeasible
*/
int ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
  double feasibilityTolerance,
  bool keepIntegers,
  int numberPasses,
  bool dropNames,
  bool doRowObjective)
{
  // Check matrix
  if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
        1.0e20))
    return 2;
  saveFile_ = fileName;
  si.saveModel(saveFile_.c_str());
  ClpSimplex *model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
    doRowObjective);
  if (model == &si) {
    return 0;
  } else {
    si.restoreModel(saveFile_.c_str());
    remove(saveFile_.c_str());
    return 1;
  }
}
#endif
// Return pointer to presolved model
ClpSimplex *
ClpPresolve::model() const
{
  return presolvedModel_;
}
// Return pointer to original model
ClpSimplex *
ClpPresolve::originalModel() const
{
  return originalModel_;
}
// Return presolve status (0,1,2)
int ClpPresolve::presolveStatus() const
{
  if (nelems_ >= 0) {
    // feasible (or not done yet)
    return 0;
  } else {
    int presolveStatus = -static_cast< int >(nelems_);
    // If both infeasible and unbounded - say infeasible
    if (presolveStatus > 2)
      presolveStatus = 1;
    return presolveStatus;
  }
}

static
CoinPostsolveMatrix* create_CoinPostsolveMatrix(ClpSimplex *si,
  int ncols0_in,
  int nrows0_in,
  CoinBigIndex nelems0,

  double maxmin,
  // end prepost members

  double *sol_in,
  double *acts_in,

  unsigned char *colstat_in,
  unsigned char *rowstat_in);

void ClpPresolve::postsolve(bool updateStatus)
{
  // Return at once if no presolved model
  if (!presolvedModel_)
    return;
  // Messages
  CoinMessages messages = originalModel_->coinMessages();
  if (!presolvedModel_->isProvenOptimal()) {
    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
      messages)
      << CoinMessageEol;
  }

  // this is the size of the original problem
  const int ncols0 = ncols_;
  const int nrows0 = nrows_;
  const CoinBigIndex nelems0 = nelems_;

  // this is the reduced problem
  int ncols = presolvedModel_->getNumCols();
  int nrows = presolvedModel_->getNumRows();

  double *acts = NULL;
  double *sol = NULL;
  unsigned char *rowstat = NULL;
  unsigned char *colstat = NULL;
#ifndef CLP_NO_STD
  if (saveFile_ == "") {
#endif
    // reality check
    assert(ncols0 == originalModel_->getNumCols());
    assert(nrows0 == originalModel_->getNumRows());
    acts = originalModel_->primalRowSolution();
    sol = originalModel_->primalColumnSolution();
    if (updateStatus) {
      // postsolve does not know about fixed
      int i;
      for (i = 0; i < nrows + ncols; i++) {
        if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
          presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
      }
      unsigned char *status = originalModel_->statusArray();
      if (!status) {
        originalModel_->createStatus();
        status = originalModel_->statusArray();
      }
      rowstat = status + ncols0;
      colstat = status;
      CoinMemcpyN(presolvedModel_->statusArray(), ncols, colstat);
      CoinMemcpyN(presolvedModel_->statusArray() + ncols, nrows, rowstat);
    }
#ifndef CLP_NO_STD
  } else {
    // from file
    acts = new double[nrows0];
    sol = new double[ncols0];
    CoinZeroN(acts, nrows0);
    CoinZeroN(sol, ncols0);
    if (updateStatus) {
      unsigned char *status = new unsigned char[nrows0 + ncols0];
      rowstat = status + ncols0;
      colstat = status;
      CoinMemcpyN(presolvedModel_->statusArray(), ncols, colstat);
      CoinMemcpyN(presolvedModel_->statusArray() + ncols, nrows, rowstat);
    }
  }
#endif

  // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
  // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
  // cause duplicate free. In case where saveFile == "", as best I can see
  // arrays are owned by originalModel_. fix is to
  // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
  CoinPostsolveMatrix* prob = create_CoinPostsolveMatrix(presolvedModel_,
    ncols0,
    nrows0,
    nelems0,
    presolvedModel_->getObjSense(),
    // end prepost

    sol, acts,
    colstat, rowstat);

  postsolve(*prob);

#ifndef CLP_NO_STD
  if (saveFile_ != "") {
    // From file
    assert(originalModel_ == presolvedModel_);
    originalModel_->restoreModel(saveFile_.c_str());
    remove(saveFile_.c_str());
    CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
    // delete [] acts;
    CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
    // delete [] sol;
    if (updateStatus) {
      CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
      // delete [] colstat;
    }
  } else {
#endif
    prob->sol_ = 0;
    prob->acts_ = 0;
    prob->colstat_ = 0;
#ifndef CLP_NO_STD
  }
#endif
  // put back duals
  CoinMemcpyN(prob->rowduals_, nrows_, originalModel_->dualRowSolution());
  double maxmin = originalModel_->getObjSense();
  if (maxmin < 0.0) {
    // swap signs
    int i;
    double *pi = originalModel_->dualRowSolution();
    for (i = 0; i < nrows_; i++)
      pi[i] = -pi[i];
  }
  // Now check solution
  double offset;
  CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
                originalModel_->primalColumnSolution(), offset, true),
    ncols_, originalModel_->dualColumnSolution());
  originalModel_->clpMatrix()->transposeTimes(-1.0,
    originalModel_->dualRowSolution(),
    originalModel_->dualColumnSolution());
  memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
  originalModel_->clpMatrix()->times(1.0,
    originalModel_->primalColumnSolution(),
    originalModel_->primalRowSolution());
  originalModel_->checkSolutionInternal();
  if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
    // See if we can fix easily
    static_cast< ClpSimplexOther * >(originalModel_)->cleanupAfterPostsolve();
  }
  // Messages
  presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
    messages)
    << originalModel_->objectiveValue()
    << originalModel_->sumDualInfeasibilities()
    << originalModel_->numberDualInfeasibilities()
    << originalModel_->sumPrimalInfeasibilities()
    << originalModel_->numberPrimalInfeasibilities()
    << CoinMessageEol;

  //originalModel_->objectiveValue_=objectiveValue_;
  originalModel_->setNumberIterations(presolvedModel_->numberIterations());
  if (!presolvedModel_->status()) {
    if (!originalModel_->numberDualInfeasibilities() && !originalModel_->numberPrimalInfeasibilities()) {
      originalModel_->setProblemStatus(0);
    } else {
      originalModel_->setProblemStatus(-1);
      // Say not optimal after presolve
      originalModel_->setSecondaryStatus(7);
      presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
        messages)
        << CoinMessageEol;
    }
  } else {
    originalModel_->setProblemStatus(presolvedModel_->status());
    // but not if close to feasible
    if (originalModel_->sumPrimalInfeasibilities() < 1.0e-1) {
      originalModel_->setProblemStatus(-1);
      // Say not optimal after presolve
      originalModel_->setSecondaryStatus(7);
    }
  }
#ifndef CLP_NO_STD
  if (saveFile_ != "")
    presolvedModel_ = NULL;
#endif

  delete prob;
}

// return pointer to original columns
const int *
ClpPresolve::originalColumns() const
{
  return originalColumn_;
}
// return pointer to original rows
const int *
ClpPresolve::originalRows() const
{
  return originalRow_;
}
// Set pointer to original model
void ClpPresolve::setOriginalModel(ClpSimplex *model)
{
  originalModel_ = model;
}
#if 0
// A lazy way to restrict which transformations are applied
// during debugging.
static int ATOI(const char *name)
{
     return true;
#if PRESOLVE_DEBUG || PRESOLVE_SUMMARY
     if (getenv(name)) {
          int val = atoi(getenv(name));
          printf("%s = %d\n", name, val);
          return (val);
     } else {
          if (strcmp(name, "off"))
               return (true);
          else
               return (false);
     }
#else
     return (true);
#endif
}
#endif
#ifdef PRESOLVE_DEBUG
#define PRESOLVE_CHECK_SOL 1
#endif
//#define PRESOLVE_CHECK_SOL 1
#if PRESOLVE_CHECK_SOL
void check_sol(CoinPresolveMatrix *prob, double tol)
{
  double *colels = prob->colels_;
  int *hrow = prob->hrow_;
  int *mcstrt = prob->mcstrt_;
  int *hincol = prob->hincol_;
  int *hinrow = prob->hinrow_;
  int ncols = prob->ncols_;

  double *csol = prob->sol_;
  double *acts = prob->acts_;
  double *clo = prob->clo_;
  double *cup = prob->cup_;
  int nrows = prob->nrows_;
  double *rlo = prob->rlo_;
  double *rup = prob->rup_;

  int colx;

  double *rsol = new double[nrows];
  memset(rsol, 0, nrows * sizeof(double));

  for (colx = 0; colx < ncols; ++colx) {
    if (1) {
      CoinBigIndex k = mcstrt[colx];
      int nx = hincol[colx];
      double solutionValue = csol[colx];
      for (int i = 0; i < nx; ++i) {
        int row = hrow[k];
        double coeff = colels[k];
        k++;
        rsol[row] += solutionValue * coeff;
      }
      if (csol[colx] < clo[colx] - tol) {
        printf("low CSOL:  %d  - %g %g %g\n",
          colx, clo[colx], csol[colx], cup[colx]);
      } else if (csol[colx] > cup[colx] + tol) {
        printf("high CSOL:  %d  - %g %g %g\n",
          colx, clo[colx], csol[colx], cup[colx]);
      }
    }
  }
  int rowx;
  for (rowx = 0; rowx < nrows; ++rowx) {
    if (hinrow[rowx]) {
      if (fabs(rsol[rowx] - acts[rowx]) > tol)
        printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
          rowx, rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
      if (rsol[rowx] < rlo[rowx] - tol) {
        printf("low RSOL:  %d - %g %g %g\n",
          rowx, rlo[rowx], rsol[rowx], rup[rowx]);
      } else if (rsol[rowx] > rup[rowx] + tol) {
        printf("high RSOL:  %d - %g %g %g\n",
          rowx, rlo[rowx], rsol[rowx], rup[rowx]);
      }
    }
  }
  delete[] rsol;
}
#endif
static int tightenDoubletons2(CoinPresolveMatrix *prob)
{
  // column-major representation
  const int ncols = prob->ncols_;
  const CoinBigIndex *const mcstrt = prob->mcstrt_;
  const int *const hincol = prob->hincol_;
  const int *const hrow = prob->hrow_;
  double *colels = prob->colels_;
  double *cost = prob->cost_;

  // column type, bounds, solution, and status
  const unsigned char *const integerType = prob->integerType_;
  double *clo = prob->clo_;
  double *cup = prob->cup_;
  // row-major representation
  //const int nrows = prob->nrows_ ;
  const CoinBigIndex *const mrstrt = prob->mrstrt_;
  const int *const hinrow = prob->hinrow_;
  const int *const hcol = prob->hcol_;
  double *rowels = prob->rowels_;

  // row bounds
  double *const rlo = prob->rlo_;
  double *const rup = prob->rup_;

  // tolerances
  //const double ekkinf2 = PRESOLVE_SMALL_INF ;
  //const double ekkinf = ekkinf2*1.0e8 ;
  //const double ztolcbarj = prob->ztoldj_ ;
  //const CoinRelFltEq relEq(prob->ztolzb_) ;
  int numberChanged = 0;
  double bound[2];
  double alpha[2] = { 0.0, 0.0 };
  double offset = 0.0;

  for (int icol = 0; icol < ncols; icol++) {
    if (hincol[icol] == 2) {
      CoinBigIndex start = mcstrt[icol];
      int row0 = hrow[start];
      if (hinrow[row0] != 2)
        continue;
      int row1 = hrow[start + 1];
      if (hinrow[row1] != 2)
        continue;
      double element0 = colels[start];
      double rowUpper0 = rup[row0];
      bool swapSigns0 = false;
      if (rlo[row0] > -1.0e30) {
        if (rup[row0] > 1.0e30) {
          swapSigns0 = true;
          rowUpper0 = -rlo[row0];
          element0 = -element0;
        } else {
          // range or equality
          continue;
        }
      } else if (rup[row0] > 1.0e30) {
        // free
        continue;
      }
#if 0
      // skip here for speed
      // skip if no cost (should be able to get rid of)
      if (!cost[icol]) {
	printf("should be able to get rid of %d with no cost\n",icol);
	continue;
      }
      // skip if negative cost for now
      if (cost[icol]<0.0) {
	printf("code for negative cost\n");
	continue;
      }
#endif
      double element1 = colels[start + 1];
      double rowUpper1 = rup[row1];
      bool swapSigns1 = false;
      if (rlo[row1] > -1.0e30) {
        if (rup[row1] > 1.0e30) {
          swapSigns1 = true;
          rowUpper1 = -rlo[row1];
          element1 = -element1;
        } else {
          // range or equality
          continue;
        }
      } else if (rup[row1] > 1.0e30) {
        // free
        continue;
      }
      double lowerX = clo[icol];
      double upperX = cup[icol];
      int otherCol = -1;
      CoinBigIndex startRow = mrstrt[row0];
      for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
        int jcol = hcol[j];
        if (jcol != icol) {
          alpha[0] = swapSigns0 ? -rowels[j] : rowels[j];
          otherCol = jcol;
        }
      }
      startRow = mrstrt[row1];
      bool possible = true;
      for (CoinBigIndex j = startRow; j < startRow + 2; j++) {
        int jcol = hcol[j];
        if (jcol != icol) {
          if (jcol == otherCol) {
            alpha[1] = swapSigns1 ? -rowels[j] : rowels[j];
          } else {
            possible = false;
          }
        }
      }
      if (possible) {
        // skip if no cost (should be able to get rid of)
        if (!cost[icol]) {
          PRESOLVE_DETAIL_PRINT(printf("should be able to get rid of %d with no cost\n", icol));
          continue;
        }
        // skip if negative cost for now
        if (cost[icol] < 0.0) {
          PRESOLVE_DETAIL_PRINT(printf("code for negative cost\n"));
          continue;
        }
        bound[0] = clo[otherCol];
        bound[1] = cup[otherCol];
        double lowestLowest = COIN_DBL_MAX;
        double highestLowest = -COIN_DBL_MAX;
        double lowestHighest = COIN_DBL_MAX;
        double highestHighest = -COIN_DBL_MAX;
        int binding0 = 0;
        int binding1 = 0;
        for (int k = 0; k < 2; k++) {
          bool infLow0 = false;
          bool infLow1 = false;
          double sum0 = 0.0;
          double sum1 = 0.0;
          double value = bound[k];
          if (fabs(value) < 1.0e30) {
            sum0 += alpha[0] * value;
            sum1 += alpha[1] * value;
          } else {
            if (alpha[0] > 0.0) {
              if (value < 0.0)
                infLow0 = true;
            } else if (alpha[0] < 0.0) {
              if (value > 0.0)
                infLow0 = true;
            }
            if (alpha[1] > 0.0) {
              if (value < 0.0)
                infLow1 = true;
            } else if (alpha[1] < 0.0) {
              if (value > 0.0)
                infLow1 = true;
            }
          }
          /* Got sums
	   */
          double thisLowest0 = -COIN_DBL_MAX;
          double thisHighest0 = COIN_DBL_MAX;
          if (element0 > 0.0) {
            // upper bound unless inf&2 !=0
            if (!infLow0)
              thisHighest0 = (rowUpper0 - sum0) / element0;
          } else {
            // lower bound unless inf&2 !=0
            if (!infLow0)
              thisLowest0 = (rowUpper0 - sum0) / element0;
          }
          double thisLowest1 = -COIN_DBL_MAX;
          double thisHighest1 = COIN_DBL_MAX;
          if (element1 > 0.0) {
            // upper bound unless inf&2 !=0
            if (!infLow1)
              thisHighest1 = (rowUpper1 - sum1) / element1;
          } else {
            // lower bound unless inf&2 !=0
            if (!infLow1)
              thisLowest1 = (rowUpper1 - sum1) / element1;
          }
          if (thisLowest0 > thisLowest1 + 1.0e-12) {
            if (thisLowest0 > lowerX + 1.0e-12)
              binding0 |= 1 << k;
          } else if (thisLowest1 > thisLowest0 + 1.0e-12) {
            if (thisLowest1 > lowerX + 1.0e-12)
              binding1 |= 1 << k;
            thisLowest0 = thisLowest1;
          }
          if (thisHighest0 < thisHighest1 - 1.0e-12) {
            if (thisHighest0 < upperX - 1.0e-12)
              binding0 |= 1 << k;
          } else if (thisHighest1 < thisHighest0 - 1.0e-12) {
            if (thisHighest1 < upperX - 1.0e-12)
              binding1 |= 1 << k;
            thisHighest0 = thisHighest1;
          }
          lowestLowest = CoinMin(lowestLowest, thisLowest0);
          highestHighest = CoinMax(highestHighest, thisHighest0);
          lowestHighest = CoinMin(lowestHighest, thisHighest0);
          highestLowest = CoinMax(highestLowest, thisLowest0);
        }
        // see if any good
        //#define PRINT_VALUES
        if (!binding0 || !binding1) {
          PRESOLVE_DETAIL_PRINT(printf("Row redundant for column %d\n", icol));
        } else {
#ifdef PRINT_VALUES
          printf("Column %d bounds %g,%g lowest %g,%g highest %g,%g\n",
            icol, lowerX, upperX, lowestLowest, lowestHighest,
            highestLowest, highestHighest);
#endif
          // if integer adjust
          if (integerType[icol]) {
            lowestLowest = ceil(lowestLowest - 1.0e-5);
            highestLowest = ceil(highestLowest - 1.0e-5);
            lowestHighest = floor(lowestHighest + 1.0e-5);
            highestHighest = floor(highestHighest + 1.0e-5);
          }
          // if costed may be able to adjust
          if (cost[icol] >= 0.0) {
            if (highestLowest < upperX && highestLowest >= lowerX && highestHighest < 1.0e30) {
              highestHighest = CoinMin(highestHighest, highestLowest);
            }
          }
          if (cost[icol] <= 0.0) {
            if (lowestHighest > lowerX && lowestHighest <= upperX && lowestHighest > -1.0e30) {
              lowestLowest = CoinMax(lowestLowest, lowestHighest);
            }
          }
#if 1
          if (lowestLowest > lowerX + 1.0e-8) {
#ifdef PRINT_VALUES
            printf("Can increase lower bound on %d from %g to %g\n",
              icol, lowerX, lowestLowest);
#endif
            lowerX = lowestLowest;
          }
          if (highestHighest < upperX - 1.0e-8) {
#ifdef PRINT_VALUES
            printf("Can decrease upper bound on %d from %g to %g\n",
              icol, upperX, highestHighest);
#endif
            upperX = highestHighest;
          }
#endif
          // see if we can move costs
          double xValue;
          double yValue0;
          double yValue1;
          double newLower = COIN_DBL_MAX;
          double newUpper = -COIN_DBL_MAX;
#ifdef PRINT_VALUES
          double ranges0[2];
          double ranges1[2];
#endif
          double costEqual;
          double slope[2];
          assert(binding0 + binding1 == 3);
          // get where equal
          xValue = (rowUpper0 * element1 - rowUpper1 * element0) / (alpha[0] * element1 - alpha[1] * element0);
          yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
          yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
          newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
          newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
          double xValueEqual = xValue;
          double yValueEqual = yValue0;
          costEqual = xValue * cost[otherCol] + yValueEqual * cost[icol];
          if (binding0 == 1) {
#ifdef PRINT_VALUES
            ranges0[0] = bound[0];
            ranges0[1] = yValue0;
            ranges1[0] = yValue0;
            ranges1[1] = bound[1];
#endif
            // take x 1.0 down
            double x = xValue - 1.0;
            double y = (rowUpper0 - x * alpha[0]) / element0;
            double costTotal = x * cost[otherCol] + y * cost[icol];
            slope[0] = costEqual - costTotal;
            // take x 1.0 up
            x = xValue + 1.0;
            y = (rowUpper1 - x * alpha[1]) / element0;
            costTotal = x * cost[otherCol] + y * cost[icol];
            slope[1] = costTotal - costEqual;
          } else {
#ifdef PRINT_VALUES
            ranges1[0] = bound[0];
            ranges1[1] = yValue0;
            ranges0[0] = yValue0;
            ranges0[1] = bound[1];
#endif
            // take x 1.0 down
            double x = xValue - 1.0;
            double y = (rowUpper1 - x * alpha[1]) / element0;
            double costTotal = x * cost[otherCol] + y * cost[icol];
            slope[1] = costEqual - costTotal;
            // take x 1.0 up
            x = xValue + 1.0;
            y = (rowUpper0 - x * alpha[0]) / element0;
            costTotal = x * cost[otherCol] + y * cost[icol];
            slope[0] = costTotal - costEqual;
          }
#ifdef PRINT_VALUES
          printf("equal value of %d is %g, value of %d is max(%g,%g) - %g\n",
            otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
          printf("Cost at equality %g for constraint 0 ranges %g -> %g slope %g for constraint 1 ranges %g -> %g slope %g\n",
            costEqual, ranges0[0], ranges0[1], slope[0], ranges1[0], ranges1[1], slope[1]);
#endif
          xValue = bound[0];
          yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
          yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
#ifdef PRINT_VALUES
          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
            otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
#endif
          newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
          // cost>0 so will be at lower
          //double yValueAtBound0=newLower;
          newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
          xValue = bound[1];
          yValue0 = (rowUpper0 - xValue * alpha[0]) / element0;
          yValue1 = (rowUpper1 - xValue * alpha[1]) / element1;
#ifdef PRINT_VALUES
          printf("value of %d is %g, value of %d is max(%g,%g) - %g\n",
            otherCol, xValue, icol, yValue0, yValue1, CoinMax(yValue0, yValue1));
#endif
          newLower = CoinMin(newLower, CoinMax(yValue0, yValue1));
          // cost>0 so will be at lower
          //double yValueAtBound1=newLower;
          newUpper = CoinMax(newUpper, CoinMax(yValue0, yValue1));
          lowerX = CoinMax(lowerX, newLower - 1.0e-12 * fabs(newLower));
          upperX = CoinMin(upperX, newUpper + 1.0e-12 * fabs(newUpper));
          // Now make duplicate row
          // keep row 0 so need to adjust costs so same
#ifdef PRINT_VALUES
          printf("Costs for x %g,%g,%g are %g,%g,%g\n",
            xValueEqual - 1.0, xValueEqual, xValueEqual + 1.0,
            costEqual - slope[0], costEqual, costEqual + slope[1]);
#endif
          double costOther = cost[otherCol] + slope[1];
          double costThis = cost[icol] + slope[1] * (element0 / alpha[0]);
          xValue = xValueEqual;
          yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
          double thisOffset = costEqual - (costOther * xValue + costThis * yValue0);
          offset += thisOffset;
#ifdef PRINT_VALUES
          printf("new cost at equal %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
#endif
          xValue = xValueEqual - 1.0;
          yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
#ifdef PRINT_VALUES
          printf("new cost at -1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
#endif
          assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual - slope[0])) < 1.0e-5);
          xValue = xValueEqual + 1.0;
          yValue0 = CoinMax((rowUpper0 - xValue * alpha[0]) / element0, lowerX);
#ifdef PRINT_VALUES
          printf("new cost at +1 %g\n", costOther * xValue + costThis * yValue0 + thisOffset);
#endif
          assert(fabs((costOther * xValue + costThis * yValue0 + thisOffset) - (costEqual + slope[1])) < 1.0e-5);
          numberChanged++;
          //	  continue;
          cost[otherCol] = costOther;
          cost[icol] = costThis;
          clo[icol] = lowerX;
          cup[icol] = upperX;
          CoinBigIndex startCol[2];
          CoinBigIndex endCol[2];
          startCol[0] = mcstrt[icol];
          endCol[0] = startCol[0] + 2;
          startCol[1] = mcstrt[otherCol];
          endCol[1] = startCol[1] + hincol[otherCol];
          double values[2] = { 0.0, 0.0 };
          for (int k = 0; k < 2; k++) {
            for (CoinBigIndex i = startCol[k]; i < endCol[k]; i++) {
              if (hrow[i] == row0)
                values[k] = colels[i];
            }
            for (CoinBigIndex i = startCol[k]; i < endCol[k]; i++) {
              if (hrow[i] == row1)
                colels[i] = values[k];
            }
          }
          for (CoinBigIndex i = mrstrt[row1]; i < mrstrt[row1] + 2; i++) {
            if (hcol[i] == icol)
              rowels[i] = values[0];
            else
              rowels[i] = values[1];
          }
        }
      }
    }
  }
#if ABC_NORMAL_DEBUG > 0
  if (offset)
    printf("Cost offset %g\n", offset);
#endif
  return numberChanged;
}
//#define COIN_PRESOLVE_BUG
#ifdef COIN_PRESOLVE_BUG
static int counter = 1000000;
static int startEmptyRows = 0;
static int startEmptyColumns = 0;
static bool break2(CoinPresolveMatrix *prob)
{
  int droppedRows = prob->countEmptyRows() - startEmptyRows;
  int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
  startEmptyRows = prob->countEmptyRows();
  startEmptyColumns = prob->countEmptyCols();
  printf("Dropped %d rows and %d columns - current empty %d, %d\n", droppedRows,
    droppedColumns, startEmptyRows, startEmptyColumns);
  counter--;
  if (!counter) {
    printf("skipping next and all\n");
  }
  return (counter <= 0);
}
#define possibleBreak \
  if (break2(prob))   \
  break
#define possibleSkip if (!break2(prob))
#else
#define possibleBreak
#define possibleSkip
#endif
#define SOME_PRESOLVE_DETAIL
#ifndef SOME_PRESOLVE_DETAIL
#define printProgress(x, y) \
  {                         \
  }
#else
#define printProgress(x, y)                                                                \
  {                                                                                        \
    if ((presolveActions_ & 0x80000000) != 0)                                              \
      printf("%c loop %d %d empty rows, %d empty columns\n", x, y, prob->countEmptyRows(), \
        prob->countEmptyCols());                                                           \
  }
#endif
// This is the presolve loop.
// It is a separate virtual function so that it can be easily
// customized by subclassing CoinPresolve.
const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
{
  // Messages
  CoinMessages messages = CoinMessage(prob->messages().language());
  paction_ = 0;
  prob->maxSubstLevel_ = CoinMax(3, prob->maxSubstLevel_);
#ifndef PRESOLVE_DETAIL
  if (prob->tuning_) {
#endif
    int numberEmptyRows = 0;
    for (int i = 0; i < prob->nrows_; i++) {
      if (!prob->hinrow_[i]) {
        PRESOLVE_DETAIL_PRINT(printf("pre_empty row %d\n", i));
        //printf("pre_empty row %d\n",i);
        numberEmptyRows++;
      }
    }
    int numberEmptyCols = 0;
    for (int i = 0; i < prob->ncols_; i++) {
      if (!prob->hincol_[i]) {
        PRESOLVE_DETAIL_PRINT(printf("pre_empty col %d\n", i));
        //printf("pre_empty col %d\n",i);
        numberEmptyCols++;
      }
    }
    printf("CoinPresolve initial state %d empty rows and %d empty columns\n",
      numberEmptyRows, numberEmptyCols);
#ifndef PRESOLVE_DETAIL
  }
#endif
  prob->status_ = 0; // say feasible
  printProgress('A', 0);
  paction_ = make_fixed(prob, paction_);
  paction_ = testRedundant(prob, paction_);
  printProgress('B', 0);
  // if integers then switch off dual stuff
  // later just do individually
  bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
  // but allow in some cases
  if ((presolveActions_ & 512) != 0)
    doDualStuff = true;
  if (prob->anyProhibited())
    doDualStuff = false;
  if (!doDual())
    doDualStuff = false;
#if PRESOLVE_CONSISTENCY
  //  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
  presolve_links_ok(prob, false, true);
#endif

  if (!prob->status_) {
    bool slackSingleton = doSingletonColumn();
    slackSingleton = true;
    const bool slackd = doSingleton();
    const bool doubleton = doDoubleton();
    const bool tripleton = doTripleton();
    //#define NO_FORCING
#ifndef NO_FORCING
    const bool forcing = doForcing();
#endif
    const bool ifree = doImpliedFree();
    const bool zerocost = doTighten();
    const bool dupcol = doDupcol();
    const bool duprow = doDuprow();
    const bool dual = doDualStuff;
    // Whether we want to allow duplicate intersections
    if (doIntersection())
      prob->presolveOptions_ |= 0x10;
    // zero small elements in aggregation
    prob->presolveOptions_ |= zeroSmall() * 0x20000;
    // some things are expensive so just do once (normally)

    int i;
    // say look at all
    if (!prob->anyProhibited()) {
      for (i = 0; i < nrows_; i++)
        prob->rowsToDo_[i] = i;
      prob->numberRowsToDo_ = nrows_;
      for (i = 0; i < ncols_; i++)
        prob->colsToDo_[i] = i;
      prob->numberColsToDo_ = ncols_;
    } else {
      // some stuff must be left alone
      prob->numberRowsToDo_ = 0;
      for (i = 0; i < nrows_; i++)
        if (!prob->rowProhibited(i))
          prob->rowsToDo_[prob->numberRowsToDo_++] = i;
      prob->numberColsToDo_ = 0;
      for (i = 0; i < ncols_; i++)
        if (!prob->colProhibited(i))
          prob->colsToDo_[prob->numberColsToDo_++] = i;
    }

    // transfer costs (may want to do it in OsiPresolve)
    // need a transfer back at end of postsolve transferCosts(prob);

    int iLoop;
#if PRESOLVE_CHECK_SOL
    check_sol(prob, 1.0e0);
#endif
    if (dupcol) {
      // maybe allow integer columns to be checked
      if ((presolveActions_ & 512) != 0)
        prob->setPresolveOptions(prob->presolveOptions() | 1);
      if ((presolveActions_ & 131072) != 0)
        prob->setPresolveOptions(prob->presolveOptions() | 0x80000);
      possibleSkip;
      paction_ = dupcol_action::presolve(prob, paction_);
      printProgress('C', 0);
    }
    if (doTwoxTwo()) {
      possibleSkip;
      paction_ = twoxtwo_action::presolve(prob, paction_);
    }
    if (duprow) {
      possibleSkip;
      if (doTwoxTwo()) {
        int nTightened = tightenDoubletons2(prob);
        if (nTightened)
          PRESOLVE_DETAIL_PRINT(printf("%d doubletons tightened\n",
            nTightened));
      }
      paction_ = duprow_action::presolve(prob, paction_);
      printProgress('D', 0);
      //paction_ = doubleton_action::presolve(prob, paction_);
      //printProgress('d',0);
      //if (doDependency()) {
      //paction_ = duprow3_action::presolve(prob, paction_);
      //printProgress('Z',0);
      //}
    }
    if (doGubrow()) {
      possibleSkip;
      paction_ = gubrow_action::presolve(prob, paction_);
      printProgress('E', 0);
    }
    if (ifree) {
      int fill_level = CoinMax(10, prob->maxSubstLevel_);
      const CoinPresolveAction *lastAction = NULL;
      int iPass = 4;
      while (lastAction != paction_ && iPass) {
        lastAction = paction_;
        paction_ = implied_free_action::presolve(prob, paction_, fill_level);
        printProgress('l', 0);
        iPass--;
      }
    }

    if ((presolveActions_ & 16384) != 0)
      prob->setPresolveOptions(prob->presolveOptions() | 16384);
    // For inaccurate data in implied free
    if ((presolveActions_ & 1024) != 0)
      prob->setPresolveOptions(prob->presolveOptions() | 0x20000);
    // Check number rows dropped
    int lastDropped = 0;
    prob->pass_ = 0;
#if CLP_INHERIT_MODE > 1
    int numberRowsStart = nrows_ - prob->countEmptyRows();
    int numberColumnsStart = ncols_ - prob->countEmptyCols();
    int numberRowsLeft = numberRowsStart;
    int numberColumnsLeft = numberColumnsStart;
    bool lastPassWasGood = true;
#if ABC_NORMAL_DEBUG
    printf("Original rows,columns %d,%d starting first pass with %d,%d\n",
      nrows_, ncols_, numberRowsLeft, numberColumnsLeft);
#endif
#endif
    if (numberPasses_ <= 5)
      prob->presolveOptions_ |= 0x10000; // say more lightweight
    for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
      // See if we want statistics
      if ((presolveActions_ & 0x80000000) != 0)
        printf("Starting major pass %d after %g seconds with %d rows, %d columns\n", iLoop + 1, CoinCpuTime() - prob->startTime_,
          nrows_ - prob->countEmptyRows(),
          ncols_ - prob->countEmptyCols());
#ifdef PRESOLVE_SUMMARY
      printf("Starting major pass %d\n", iLoop + 1);
#endif
      const CoinPresolveAction *const paction0 = paction_;
      // look for substitutions with no fill
      //#define IMPLIED 3
#ifdef IMPLIED
      int fill_level = 3;
#define IMPLIED2 99
#if IMPLIED != 3
#if IMPLIED > 2 && IMPLIED < 11
      fill_level = IMPLIED;
      COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
#endif
#if IMPLIED > 11 && IMPLIED < 21
      fill_level = -(IMPLIED - 10);
      COIN_DETAIL_PRINT(printf("** fill_level == %d !\n", fill_level));
#endif
#endif
#else
      int fill_level = prob->maxSubstLevel_;
#endif
      int whichPass = 0;
      while (1) {
        whichPass++;
        prob->pass_++;
        const CoinPresolveAction *const paction1 = paction_;

        if (slackd) {
          bool notFinished = true;
          while (notFinished) {
            possibleBreak;
            paction_ = slack_doubleton_action::presolve(prob, paction_,
              notFinished);
          }
          printProgress('F', iLoop + 1);
          if (prob->status_)
            break;
        }
        if (zerocost) {
          possibleBreak;
          paction_ = do_tighten_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('J', iLoop + 1);
        }
        if (dual && whichPass == 1) {
          // this can also make E rows so do one bit here
          possibleBreak;
          paction_ = remove_dual_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('G', iLoop + 1);
        }

        if (doubleton) {
          possibleBreak;
          paction_ = doubleton_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('H', iLoop + 1);
        }
        if (tripleton) {
          possibleBreak;
          paction_ = tripleton_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('I', iLoop + 1);
        }

#ifndef NO_FORCING
        if (forcing) {
          possibleBreak;
          paction_ = forcing_constraint_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('K', iLoop + 1);
        }
#endif

        if (ifree && (whichPass % 5) == 1) {
          possibleBreak;
          paction_ = implied_free_action::presolve(prob, paction_, fill_level);
          if (prob->status_)
            break;
          printProgress('L', iLoop + 1);
        }

#if PRESOLVE_CHECK_SOL
        check_sol(prob, 1.0e0);
#endif

#if PRESOLVE_CONSISTENCY
        //	presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
        //			  prob->nrows_);
        presolve_links_ok(prob, false, true);
#endif

//#if	PRESOLVE_DEBUG
//	presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
//			  prob->ncols_);
//#endif
//#if	PRESOLVE_CONSISTENCY
//	prob->consistent();
//#endif
#if PRESOLVE_CONSISTENCY
        presolve_no_zeros(prob, true, false);
        presolve_consistent(prob, true);
#endif

        {
          // set up for next pass
          // later do faster if many changes i.e. memset and memcpy
          const int *count = prob->hinrow_;
          const int *nextToDo = prob->nextRowsToDo_;
          int *toDo = prob->rowsToDo_;
          int nNext = prob->numberNextRowsToDo_;
          int n = 0;
          for (int i = 0; i < nNext; i++) {
            int index = nextToDo[i];
            prob->unsetRowChanged(index);
            if (count[index])
              toDo[n++] = index;
          }
          prob->numberRowsToDo_ = n;
          prob->numberNextRowsToDo_ = 0;
          count = prob->hincol_;
          nextToDo = prob->nextColsToDo_;
          toDo = prob->colsToDo_;
          nNext = prob->numberNextColsToDo_;
          n = 0;
          for (int i = 0; i < nNext; i++) {
            int index = nextToDo[i];
            prob->unsetColChanged(index);
            if (count[index])
              toDo[n++] = index;
          }
          prob->numberColsToDo_ = n;
          prob->numberNextColsToDo_ = 0;
        }
        if (paction_ == paction1 && fill_level > 0)
          break;
      }
      // say look at all
      int i;
      if (!prob->anyProhibited()) {
        const int *count = prob->hinrow_;
        int *toDo = prob->rowsToDo_;
        int n = 0;
        for (int i = 0; i < nrows_; i++) {
          prob->unsetRowChanged(i);
          if (count[i])
            toDo[n++] = i;
        }
        prob->numberRowsToDo_ = n;
        prob->numberNextRowsToDo_ = 0;
        count = prob->hincol_;
        toDo = prob->colsToDo_;
        n = 0;
        for (int i = 0; i < ncols_; i++) {
          prob->unsetColChanged(i);
          if (count[i])
            toDo[n++] = i;
        }
        prob->numberColsToDo_ = n;
        prob->numberNextColsToDo_ = 0;
      } else {
        // some stuff must be left alone
        prob->numberRowsToDo_ = 0;
        for (i = 0; i < nrows_; i++)
          if (!prob->rowProhibited(i))
            prob->rowsToDo_[prob->numberRowsToDo_++] = i;
        prob->numberColsToDo_ = 0;
        for (i = 0; i < ncols_; i++)
          if (!prob->colProhibited(i))
            prob->colsToDo_[prob->numberColsToDo_++] = i;
      }
      // now expensive things
      // this caused world.mps to run into numerical difficulties
#ifdef PRESOLVE_SUMMARY
      printf("Starting expensive\n");
#endif

      if (dual) {
        int itry;
        for (itry = 0; itry < 5; itry++) {
          possibleBreak;
          paction_ = remove_dual_action::presolve(prob, paction_);
          if (prob->status_)
            break;
          printProgress('M', iLoop + 1);
          const CoinPresolveAction *const paction2 = paction_;
          if (ifree) {
#ifdef IMPLIED
#if IMPLIED2 == 0
            int fill_level = 0; // switches off substitution
#elif IMPLIED2 != 99
            int fill_level = IMPLIED2;
#endif
#endif
            if ((itry & 1) == 0) {
              possibleBreak;
              paction_ = implied_free_action::presolve(prob, paction_, fill_level);
            }
            if (prob->status_)
              break;
            printProgress('N', iLoop + 1);
          }
          if (paction_ == paction2)
            break;
        }
      } else if (ifree) {
        // just free
#ifdef IMPLIED
#if IMPLIED2 == 0
        int fill_level = 0; // switches off substitution
#elif IMPLIED2 != 99
        int fill_level = IMPLIED2;
#endif
#endif
        possibleBreak;
        paction_ = implied_free_action::presolve(prob, paction_, fill_level);
        if (prob->status_)
          break;
        printProgress('O', iLoop + 1);
      }
#if PRESOLVE_CHECK_SOL
      check_sol(prob, 1.0e0);
#endif
      if (dupcol) {
        // maybe allow integer columns to be checked
        if ((presolveActions_ & 512) != 0)
          prob->setPresolveOptions(prob->presolveOptions() | 1);
        possibleBreak;
        paction_ = dupcol_action::presolve(prob, paction_);
        if (prob->status_)
          break;
        printProgress('P', iLoop + 1);
      }
#if PRESOLVE_CHECK_SOL
      check_sol(prob, 1.0e0);
#endif

      if (duprow) {
        possibleBreak;
        paction_ = duprow_action::presolve(prob, paction_);
        if (prob->status_)
          break;
        printProgress('Q', iLoop + 1);
      }
      // Marginally slower on netlib if this call is enabled.
      // paction_ = testRedundant(prob,paction_) ;
#if PRESOLVE_CHECK_SOL
      check_sol(prob, 1.0e0);
#endif
      bool stopLoop = false;
      {
        int *hinrow = prob->hinrow_;
        int numberDropped = 0;
        for (i = 0; i < nrows_; i++)
          if (!hinrow[i])
            numberDropped++;

        prob->messageHandler()->message(COIN_PRESOLVE_PASS,
          messages)
          << numberDropped << iLoop + 1
          << CoinMessageEol;
        //printf("%d rows dropped after pass %d\n",numberDropped,
        //     iLoop+1);
        if (numberDropped == lastDropped)
          stopLoop = true;
        else
          lastDropped = numberDropped;
      }
      // Do this here as not very loopy
      if (slackSingleton) {
        // On most passes do not touch costed slacks
        if (paction_ != paction0 && !stopLoop) {
          possibleBreak;
          paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
        } else {
          // do costed if Clp (at end as ruins rest of presolve)
          possibleBreak;
#ifndef CLP_MOVE_COSTS
          paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
#else
          double *fakeRowObjective = new double[prob->nrows_];
          memset(fakeRowObjective, 0, prob->nrows_ * sizeof(double));
          paction_ = slack_singleton_action::presolve(prob, paction_, fakeRowObjective);
          delete[] fakeRowObjective;
#endif
          stopLoop = true;
        }
        printProgress('R', iLoop + 1);
      }
#if PRESOLVE_CHECK_SOL
      check_sol(prob, 1.0e0);
#endif
      if (paction_ == paction0 || stopLoop)
        break;
#if CLP_INHERIT_MODE > 1
      // see whether to stop anyway
      int numberRowsNow = nrows_ - prob->countEmptyRows();
      int numberColumnsNow = ncols_ - prob->countEmptyCols();
#if ABC_NORMAL_DEBUG
      printf("Original rows,columns %d,%d - last %d,%d end of pass %d has %d,%d\n",
        nrows_, ncols_, numberRowsLeft, numberColumnsLeft, iLoop + 1, numberRowsNow,
        numberColumnsNow);
#endif
      int rowsDeleted = numberRowsLeft - numberRowsNow;
      int columnsDeleted = numberColumnsLeft - numberColumnsNow;
      if (iLoop > 15) {
        if (rowsDeleted * 100 < numberRowsStart && columnsDeleted * 100 < numberColumnsStart)
          break;
        lastPassWasGood = true;
      } else if (rowsDeleted * 100 < numberRowsStart && rowsDeleted < 500 && columnsDeleted * 100 < numberColumnsStart && columnsDeleted < 500) {
        if (!lastPassWasGood)
          break;
        else
          lastPassWasGood = false;
      } else {
        lastPassWasGood = true;
      }
      numberRowsLeft = numberRowsNow;
      numberColumnsLeft = numberColumnsNow;
#endif
    }
  }
  if (!prob->status_ && doDependency()) {
    paction_ = duprow3_action::presolve(prob, paction_);
    printProgress('Z', 0);
  }
  prob->presolveOptions_ &= ~0x10000;
  if (!prob->status_) {
    paction_ = drop_zero_coefficients(prob, paction_);
#if PRESOLVE_CHECK_SOL
    check_sol(prob, 1.0e0);
#endif

    paction_ = drop_empty_cols_action::presolve(prob, paction_);
    paction_ = drop_empty_rows_action::presolve(prob, paction_);
#if PRESOLVE_CHECK_SOL
    check_sol(prob, 1.0e0);
#endif
  }

  if (prob->status_) {
    if (prob->status_ == 1)
      prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
        messages)
        << prob->feasibilityTolerance_
        << CoinMessageEol;
    else if (prob->status_ == 2)
      prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
        messages)
        << CoinMessageEol;
    else
      prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
        messages)
        << CoinMessageEol;
    // get rid of data
    destroyPresolve();
  }

  return (paction_);
}

void check_djs(CoinPostsolveMatrix *prob);

// We could have implemented this by having each postsolve routine
// directly call the next one, but this may make it easier to add debugging checks.
void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
{
  {
    // Check activities
    double *colels = prob.colels_;
    int *hrow = prob.hrow_;
    CoinBigIndex *mcstrt = prob.mcstrt_;
    int *hincol = prob.hincol_;
    CoinBigIndex *link = prob.link_;
    int ncols = prob.ncols_;

    char *cdone = prob.cdone_;

    double *csol = prob.sol_;
    int nrows = prob.nrows_;

    int colx;

    double *rsol = prob.acts_;
    memset(rsol, 0, nrows * sizeof(double));

    for (colx = 0; colx < ncols; ++colx) {
      if (cdone[colx]) {
        CoinBigIndex k = mcstrt[colx];
        int nx = hincol[colx];
        double solutionValue = csol[colx];
        for (int i = 0; i < nx; ++i) {
          int row = hrow[k];
          double coeff = colels[k];
          k = link[k];
          assert(k != NO_LINK || i == nx - 1);
          rsol[row] += solutionValue * coeff;
        }
      }
    }
  }
  if (prob.maxmin_ < 0) {
    //for (int i=0;i<presolvedModel_->numberRows();i++)
    //prob.rowduals_[i]=-prob.rowduals_[i];
    for (int i = 0; i < ncols_; i++) {
      prob.cost_[i] = -prob.cost_[i];
      //prob.rcosts_[i]=-prob.rcosts_[i];
    }
    prob.maxmin_ = 1.0;
  }
  const CoinPresolveAction *paction = paction_;
  //#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG
  // Check only works after first one
  int checkit = -1;
#endif

  while (paction) {
#if PRESOLVE_DEBUG
    printf("POSTSOLVING %s\n", paction->name());
#endif

    paction->postsolve(&prob);

#if PRESOLVE_DEBUG
#if 0
          /*
	    This check fails (on exmip1 (!) in osiUnitTest) because clp
	    enters postsolve with a solution that seems to have incorrect
	    status for a logical. You can see similar behaviour with
	    column status --- incorrect entering postsolve.
	    -- lh, 111207 --
	  */
          {
               int nr = 0;
               int i;
               for (i = 0; i < prob.nrows_; i++) {
                    if ((prob.rowstat_[i] & 7) == 1) {
                         nr++;
                    } else if ((prob.rowstat_[i] & 7) == 2) {
                         // at ub
                         assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
                    } else if ((prob.rowstat_[i] & 7) == 3) {
                         // at lb
                         assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
                    }
               }
               int nc = 0;
               for (i = 0; i < prob.ncols_; i++) {
                    if ((prob.colstat_[i] & 7) == 1)
                         nc++;
               }
               printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
          }
#endif // if 0
    checkit++;
    if (prob.colstat_ && checkit > 0) {
      presolve_check_nbasic(&prob);
      presolve_check_sol(&prob, 2, 2, 1);
    }
#endif
    paction = paction->next;
#if PRESOLVE_DEBUG
    //  check_djs(&prob);
    if (checkit > 0)
      presolve_check_reduced_costs(&prob);
#endif
  }
#if PRESOLVE_DEBUG
  if (prob.colstat_) {
    presolve_check_nbasic(&prob);
    presolve_check_sol(&prob, 2, 2, 1);
  }
#endif
#undef PRESOLVE_DEBUG

#if 0 && PRESOLVE_DEBUG
     for (i = 0; i < ncols0; i++) {
          if (!cdone[i]) {
               printf("!cdone[%d]\n", i);
               abort();
          }
     }

     for (i = 0; i < nrows0; i++) {
          if (!rdone[i]) {
               printf("!rdone[%d]\n", i);
               abort();
          }
     }


     for (i = 0; i < ncols0; i++) {
          if (sol[i] < -1e10 || sol[i] > 1e10)
               printf("!!!%d %g\n", i, sol[i]);

     }

#endif

#if 0 && PRESOLVE_DEBUG
     // debug check:  make sure we ended up with same original matrix
     {
          int identical = 1;

          for (int i = 0; i < ncols0; i++) {
               PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
               CoinBigIndex kcs0 = &prob->mcstrt0[i];
               CoinBigIndex kcs = mcstrt[i];
               int n = hincol[i];
               for (int k = 0; k < n; k++) {
                    CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);

                    if (k1 == kcs + n) {
                         printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
                         abort();
                    }

                    if (colels[k1] != &prob->dels0[kcs0+k])
                         printf("BAD COLEL[%d %d %d]:  %g\n",
                                k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);

                    if (kcs0 + k != k1)
                         identical = 0;
               }
          }
          printf("identical? %d\n", identical);
     }
#endif
}

static inline double getTolerance(const ClpSimplex *si, ClpDblParam key)
{
  double tol;
  if (!si->getDblParam(key, tol)) {
    CoinPresolveAction::throwCoinError("getDblParam failed",
      "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
  }
  return (tol);
}

// Assumptions:
// 1. nrows>=m.getNumRows()
// 2. ncols>=m.getNumCols()
//
// In presolve, these values are equal.
// In postsolve, they may be inequal, since the reduced problem
// may be smaller, but we need room for the large problem.
// ncols may be larger than si.getNumCols() in postsolve,
// this at that point si will be the reduced problem,
// but we need to reserve enough space for the original problem.
static
void init_CoinPrePostsolveMatrix(
  CoinPrePostsolveMatrix* cpm,
  const ClpSimplex *si,
  int ncols_in,
  int nrows_in,
  CoinBigIndex nelems_in,
  double bulkRatio)
{
  cpm->ncols_ = si->getNumCols();
  cpm->nrows_ = si->getNumRows();
  cpm->nelems_ = si->getNumElements();
  cpm->ncols0_ = ncols_in;
  cpm->nrows0_ = nrows_in;
  cpm->bulkRatio_ = bulkRatio;
  cpm->mcstrt_ = new CoinBigIndex[ncols_in + 1];
  cpm->hincol_ = new int[ncols_in + 1];
  cpm->cost_ = new double[ncols_in];
  cpm->clo_ = new double[ncols_in];
  cpm->cup_ = new double[ncols_in];
  cpm->rlo_ = new double[nrows_in];
  cpm->rup_ = new double[nrows_in];
  cpm->originalColumn_ = new int[ncols_in];
  cpm->originalRow_ = new int[nrows_in];
  cpm->ztolzb_ = getTolerance(si, ClpPrimalTolerance);
  cpm->ztoldj_ = getTolerance(si, ClpDualTolerance);
  cpm->maxmin_ = si->getObjSense();
  cpm->sol_ = NULL;
  cpm->rowduals_ = NULL;
  cpm->acts_ = NULL;
  cpm->rcosts_ = NULL;
  cpm->colstat_ = NULL;
  cpm->rowstat_ = NULL;
  delete cpm->handler_;
  cpm->handler_ = NULL;
  cpm->defaultHandler_ = false;

  cpm->bulk0_ = static_cast< CoinBigIndex >(cpm->bulkRatio_ * CoinMax(nelems_in, cpm->nelems_)
    + ncols_in);
  // allow for temporary overflow
  cpm->hrow_ = new int[cpm->bulk0_ + ncols_in];
  cpm->colels_ = new double[cpm->bulk0_ + ncols_in];
  si->getDblParam(ClpObjOffset, cpm->originalOffset_);
  int ncols = si->getNumCols();
  int nrows = si->getNumRows();

  cpm->setMessageHandler(si->messageHandler());

  ClpDisjointCopyN(si->getColLower(), ncols, cpm->clo_);
  ClpDisjointCopyN(si->getColUpper(), ncols, cpm->cup_);
  //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
  double offset;
  ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cpm->cost_);
  ClpDisjointCopyN(si->getRowLower(), nrows, cpm->rlo_);
  ClpDisjointCopyN(si->getRowUpper(), nrows, cpm->rup_);
  int i;
  for (i = 0; i < ncols_in; i++)
	  cpm->originalColumn_[i] = i;
  for (i = 0; i < nrows_in; i++)
	  cpm->originalRow_[i] = i;
  cpm->sol_ = NULL;
  cpm->rowduals_ = NULL;
  cpm->acts_ = NULL;

  cpm->rcosts_ = NULL;
  cpm->colstat_ = NULL;
  cpm->rowstat_ = NULL;
}

// I am not familiar enough with CoinPackedMatrix to be confident
// that I will implement a row-ordered version of toColumnOrderedGapFree
// properly.
static bool isGapFree(const CoinPackedMatrix &matrix)
{
  const CoinBigIndex *start = matrix.getVectorStarts();
  const int *length = matrix.getVectorLengths();
  int i = matrix.getSizeVectorLengths() - 1;
  // Quick check
  if (matrix.getNumElements() == start[i]) {
    return true;
  } else {
    for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
      if (start[i + 1] - start[i] != length[i])
        break;
    }
    return (!(i >= 0));
  }
}
#if PRESOLVE_DEBUG
static void matrix_bounds_ok(const double *lo, const double *up, int n)
{
  int i;
  for (i = 0; i < n; i++) {
    PRESOLVEASSERT(lo[i] <= up[i]);
    PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
    PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
  }
}
#endif
static
CoinPresolveMatrix* create_CoinPresolveMatrix(
  int ncols0_in,
  double /*maxmin*/,
  // end prepost members

  ClpSimplex *si,

  // rowrep
  int nrows_in,
  CoinBigIndex nelems_in,
  bool doStatus,
  double nonLinearValue,
  double bulkRatio)
{
  CoinPresolveMatrix* cpm = new CoinPresolveMatrix(ncols0_in, nrows_in, nelems_in);
  init_CoinPrePostsolveMatrix(cpm, si,
    ncols0_in, nrows_in, nelems_in, bulkRatio);

  cpm->clink_ = new presolvehlink[ncols0_in + 1];
  cpm->rlink_ = new presolvehlink[nrows_in + 1];
  cpm->dobias_= 0.0;

  // temporary init
  cpm->integerType_ = new unsigned char[ncols0_in];
  cpm->anyInteger_ = false;
  cpm->tuning_ = false;
  cpm->startTime_ = 0.0;
  cpm->feasibilityTolerance_ = 0.0;
  cpm->status_= -1;
  cpm->pass_= 0;
  cpm->colsToDo_ = new int[ncols0_in];
  cpm->numberColsToDo_ = 0;
  cpm->nextColsToDo_ = new int[ncols0_in];
  cpm->numberNextColsToDo_ = 0;
  cpm->rowsToDo_ = new int[nrows_in];
  cpm->numberRowsToDo_ = 0;
  cpm->nextRowsToDo_ = new int[nrows_in];
  cpm->numberNextRowsToDo_ = 0;
  cpm->presolveOptions_ = 0;

  const CoinBigIndex bufsize = cpm->bulk0_;

  cpm->nrows_ = si->getNumRows();

  // Set up change bits etc
  cpm->rowChanged_ = new unsigned char[cpm->nrows_];
  memset(cpm->rowChanged_, 0, cpm->nrows_);
  cpm->colChanged_ = new unsigned char[cpm->ncols_];
  memset(cpm->colChanged_, 0, cpm->ncols_);
  CoinPackedMatrix *m = si->matrix();

  // The coefficient matrix is a big hunk of stuff.
  // Do the copy here to try to avoid running out of memory.

  const CoinBigIndex *start = m->getVectorStarts();
  const int *row = m->getIndices();
  const double *element = m->getElements();
  int icol, nel = 0;
  cpm->mcstrt_[0] = 0;
  ClpDisjointCopyN(m->getVectorLengths(), cpm->ncols_, cpm->hincol_);
  if (si->getObjSense() < 0.0) {
    for (int i = 0; i < cpm->ncols_; i++)
    	cpm->cost_[i] = -cpm->cost_[i];
    cpm->maxmin_ = 1.0;
  }
  for (icol = 0; icol < cpm->ncols_; icol++) {
    CoinBigIndex j;
    for (j = start[icol]; j < start[icol] + cpm->hincol_[icol]; j++) {
    	cpm->hrow_[nel] = row[j];
      if (fabs(element[j]) > ZTOLDP)
    	  cpm->colels_[nel++] = element[j];
    }
    cpm->mcstrt_[icol + 1] = nel;
    cpm->hincol_[icol] = static_cast< int >(nel - cpm->mcstrt_[icol]);
  }

  // same thing for row rep
  CoinPackedMatrix *mRow = new CoinPackedMatrix();
  mRow->setExtraGap(0.0);
  mRow->setExtraMajor(0.0);
  mRow->reverseOrderedCopyOf(*m);
  //mRow->removeGaps();
  //mRow->setExtraGap(0.0);

  // Now get rid of matrix
  si->createEmptyMatrix();

  double *el = mRow->getMutableElements();
  int *ind = mRow->getMutableIndices();
  CoinBigIndex *strt = mRow->getMutableVectorStarts();
  int *len = mRow->getMutableVectorLengths();
  // Do carefully to save memory
  cpm->rowels_ = new double[cpm->bulk0_];
  ClpDisjointCopyN(el, cpm->nelems_, cpm->rowels_);
  mRow->nullElementArray();
  delete[] el;
  cpm->hcol_ = new int[cpm->bulk0_];
  ClpDisjointCopyN(ind, cpm->nelems_, cpm->hcol_);
  mRow->nullIndexArray();
  delete[] ind;
  cpm->mrstrt_ = new CoinBigIndex[nrows_in + 1];
  ClpDisjointCopyN(strt, cpm->nrows_, cpm->mrstrt_);
  mRow->nullStartArray();
  cpm->mrstrt_[cpm->nrows_] = cpm->nelems_;
  delete[] strt;
  cpm->hinrow_ = new int[nrows_in + 1];
  ClpDisjointCopyN(len, cpm->nrows_, cpm->hinrow_);
  if (cpm->nelems_ > nel) {
	  cpm->nelems_ = nel;
    // Clean any small elements
    int irow;
    nel = 0;
    CoinBigIndex start = 0;
    for (irow = 0; irow < cpm->nrows_; irow++) {
      CoinBigIndex j;
      for (j = start; j < start + cpm->hinrow_[irow]; j++) {
    	  cpm->hcol_[nel] = cpm->hcol_[j];
        if (fabs(cpm->rowels_[j]) > ZTOLDP)
        	cpm->rowels_[nel++] = cpm->rowels_[j];
      }
      start = cpm->mrstrt_[irow + 1];
      cpm->mrstrt_[irow + 1] = nel;
      cpm->hinrow_[irow] = static_cast< int >(nel - cpm->mrstrt_[irow]);
    }
  }

  delete mRow;
  if (si->integerInformation()) {
    CoinMemcpyN(reinterpret_cast< unsigned char * >(si->integerInformation()), cpm->ncols_, cpm->integerType_);
  } else {
    ClpFillN< unsigned char >(cpm->integerType_, cpm->ncols_, static_cast< unsigned char >(0));
  }

#ifndef SLIM_CLP
#ifndef NO_RTTI
  ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(si->objectiveAsObject()));
#else
  ClpQuadraticObjective *quadraticObj = NULL;
  if (si->objectiveAsObject()->type() == 2)
    quadraticObj = (static_cast< ClpQuadraticObjective * >(si->objectiveAsObject()));
#endif
#endif
  // Set up prohibited bits if needed
  if (nonLinearValue) {
	  cpm->anyProhibited_ = true;
    for (icol = 0; icol < cpm->ncols_; icol++) {
      CoinBigIndex j;
      bool nonLinearColumn = false;
      if (cpm->cost_[icol] == nonLinearValue)
        nonLinearColumn = true;
      for (j = cpm->mcstrt_[icol]; j < cpm->mcstrt_[icol + 1]; j++) {
        if (cpm->colels_[j] == nonLinearValue) {
          nonLinearColumn = true;
          cpm->setRowProhibited(cpm->hrow_[j]);
        }
      }
      if (nonLinearColumn)
    	  cpm->setColProhibited(icol);
    }
#ifndef SLIM_CLP
  } else if (quadraticObj) {
    CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
    //const int * columnQuadratic = quadratic->getIndices();
    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    const int *columnQuadraticLength = quadratic->getVectorLengths();
    //double * quadraticElement = quadratic->getMutableElements();
    int numberColumns = quadratic->getNumCols();
    cpm->anyProhibited_ = true;
    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
      if (columnQuadraticLength[iColumn]) {
    	  cpm->setColProhibited(iColumn);
        //printf("%d prohib\n",iColumn);
      }
    }
#endif
  } else {
	  cpm->anyProhibited_ = false;
  }

  if (doStatus) {
    // allow for status and solution
	  cpm->sol_ = new double[cpm->ncols_];
    CoinMemcpyN(si->primalColumnSolution(), cpm->ncols_, cpm->sol_);
    ;
    cpm->acts_ = new double[cpm->nrows_];
    CoinMemcpyN(si->primalRowSolution(), cpm->nrows_, cpm->acts_);
    if (!si->statusArray())
      si->createStatus();
    cpm->colstat_ = new unsigned char[cpm->nrows_ + cpm->ncols_];
    CoinMemcpyN(si->statusArray(), (cpm->nrows_ + cpm->ncols_), cpm->colstat_);
    cpm->rowstat_ = cpm->colstat_ + cpm->ncols_;
  }

  // the original model's fields are now unneeded - free them

  si->resize(0, 0);

#if PRESOLVE_DEBUG
  matrix_bounds_ok(cpm->rlo_, cpm->rup_, cpm->nrows_);
  matrix_bounds_ok(cpm->clo_, cpm->cup_, cpm->ncols_);
#endif

#if 0
     for (i = 0; i < nrows; ++i)
          printf("NR: %6d\n", cpm->hinrow[i]);
     for (int i = 0; i < ncols; ++i)
          printf("NC: %6d\n", cpm->hincol[i]);
#endif

  presolve_make_memlists(/*mcstrt_,*/ cpm->hincol_, cpm->clink_, cpm->ncols_);
  presolve_make_memlists(/*mrstrt_,*/ cpm->hinrow_, cpm->rlink_, cpm->nrows_);

  // this allows last col/row to expand up to bufsize-1 (22);
  // this must come after the calls to presolve_prefix
  cpm->mcstrt_[cpm->ncols_] = bufsize - 1;
  cpm->mrstrt_[cpm->nrows_] = bufsize - 1;
  // Allocate useful arrays
  cpm->initializeStuff();

#if PRESOLVE_CONSISTENCY
  //consistent(false);
  cpm->presolve_consistent(this, false);
#endif

  return cpm;
}

// avoid compiler warnings
#if PRESOLVE_SUMMARY > 0
static
void update_model(CoinPresolveMatrix& cpm, ClpSimplex *si,
  int nrows0, int ncols0,
  CoinBigIndex nelems0)
#else
static
void update_model(CoinPresolveMatrix& cpm, ClpSimplex *si,
  int /*nrows0*/,
  int /*ncols0*/,
  CoinBigIndex /*nelems0*/)
#endif
{
  if (si->getObjSense() < 0.0) {
    for (int i = 0; i < cpm.ncols_; i++)
    	cpm.cost_[i] = -cpm.cost_[i];
    cpm.dobias_ = -cpm.dobias_;
  }
  si->loadProblem(cpm.ncols_, cpm.nrows_, cpm.mcstrt_, cpm.hrow_, cpm.colels_, cpm.hincol_,
		  cpm.clo_, cpm.cup_, cpm.cost_, cpm.rlo_, cpm.rup_);
  //delete [] si->integerInformation();
  int numberIntegers = 0;
  for (int i = 0; i < cpm.ncols_; i++) {
    if (cpm.integerType_[i])
      numberIntegers++;
  }
  if (numberIntegers)
    si->copyInIntegerInformation(reinterpret_cast< const char * >(cpm.integerType_));
  else
    si->copyInIntegerInformation(NULL);

#if PRESOLVE_SUMMARY
  printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
		  cpm.ncols_, cpm.ncols0 - cpm.ncols_,
		  cpm.nrows_, cpm.nrows0 - cpm.nrows_,
    si->getNumElements(), cpm.nelems0 - si->getNumElements());
#endif
  si->setDblParam(ClpObjOffset, cpm.originalOffset_ - cpm.dobias_);
  if (si->getObjSense() < 0.0) {
    // put back
    for (int i = 0; i < cpm.ncols_; i++)
    	cpm.cost_[i] = -cpm.cost_[i];
    cpm.dobias_ = -cpm.dobias_;
    cpm.maxmin_ = -1.0;
  }
}

////////////////  POSTSOLVE

static
CoinPostsolveMatrix* create_CoinPostsolveMatrix(ClpSimplex *si,
  int ncols0_in,
  int nrows0_in,
  CoinBigIndex nelems0,

  double maxmin,
  // end prepost members

  double *sol_in,
  double *acts_in,

  unsigned char *colstat_in,
  unsigned char *rowstat_in)
{
  CoinPostsolveMatrix* cpm = new CoinPostsolveMatrix(ncols0_in, nrows0_in, nelems0);
  init_CoinPrePostsolveMatrix(cpm, si,
      ncols0_in, nrows0_in, nelems0, 2.0);

  cpm->free_list_ = 0;
  // link, free_list, maxlink
  cpm->maxlink_ = cpm->bulk0_;
  cpm->link_ = new CoinBigIndex[/*maxlink*/ cpm->bulk0_];

  cpm->cdone_ = new char[cpm->ncols0_];
  cpm->rdone_ = new char[nrows0_in];

  cpm->bulk0_ = cpm->maxlink_;
  cpm->nrows_ = si->getNumRows();
  cpm->ncols_ = si->getNumCols();

  cpm->sol_ = sol_in;
  cpm->rowduals_ = NULL;
  cpm->acts_ = acts_in;

  cpm->rcosts_ = NULL;
  cpm->colstat_ = colstat_in;
  cpm->rowstat_ = rowstat_in;

  // this is the *reduced* model, which is probably smaller
  int ncols1 = cpm->ncols_;
  int nrows1 = cpm->nrows_;

  const CoinPackedMatrix *m = si->matrix();

  const CoinBigIndex nelemsr = m->getNumElements();
  if (m->getNumElements() && !isGapFree(*m)) {
    // Odd - gaps
    CoinPackedMatrix mm(*m);
    mm.removeGaps();
    mm.setExtraGap(0.0);

    ClpDisjointCopyN(mm.getVectorStarts(), ncols1, cpm->mcstrt_);
    CoinZeroN(cpm->mcstrt_ + ncols1, cpm->ncols0_ - ncols1);
    cpm->mcstrt_[ncols1] = nelems0; // ??    (should point to end of bulk store   -- lh --)
    ClpDisjointCopyN(mm.getVectorLengths(), ncols1, cpm->hincol_);
    ClpDisjointCopyN(mm.getIndices(), nelemsr, cpm->hrow_);
    ClpDisjointCopyN(mm.getElements(), nelemsr, cpm->colels_);
  } else {
    // No gaps

    ClpDisjointCopyN(m->getVectorStarts(), ncols1, cpm->mcstrt_);
    CoinZeroN(cpm->mcstrt_ + ncols1, cpm->ncols0_ - ncols1);
    cpm->mcstrt_[ncols1] = nelems0; // ??    (should point to end of bulk store   -- lh --)
    ClpDisjointCopyN(m->getVectorLengths(), ncols1, cpm->hincol_);
    ClpDisjointCopyN(m->getIndices(), nelemsr, cpm->hrow_);
    ClpDisjointCopyN(m->getElements(), nelemsr, cpm->colels_);
  }

#if 0 && PRESOLVE_DEBUG
     presolve_check_costs(model, &colcopy);
#endif

  // This determines the size of the data structure that contains
  // the matrix being postsolved.  Links are taken from the free_list
  // to recreate matrix entries that were presolved away,
  // and links are added to the free_list when entries created during
  // presolve are discarded.  There is never a need to gc this list.
  // Naturally, it should contain
  // exactly nelems0 entries "in use" when postsolving is done,
  // but I don't know whether the matrix could temporarily get
  // larger during postsolving.  Substitution into more than two
  // rows could do that, in principle.  I am being very conservative
  // here by reserving much more than the amount of space I probably need.
  // If this guess is wrong, check_free_list may be called.
  //  int bufsize = 2*nelems0;

  memset(cpm->cdone_, -1, cpm->ncols0_);
  memset(cpm->rdone_, -1, cpm->nrows0_);

  cpm->rowduals_ = new double[cpm->nrows0_];
  ClpDisjointCopyN(si->getRowPrice(), nrows1, cpm->rowduals_);

  cpm->rcosts_ = new double[cpm->ncols0_];
  ClpDisjointCopyN(si->getReducedCost(), ncols1, cpm->rcosts_);
  if (maxmin < 0.0) {
    // change so will look as if minimize
    int i;
    for (i = 0; i < nrows1; i++)
    	cpm->rowduals_[i] = -cpm->rowduals_[i];
    for (i = 0; i < ncols1; i++) {
    	cpm->rcosts_[i] = -cpm->rcosts_[i];
    }
  }

  //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
  //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);

  ClpDisjointCopyN(si->getColSolution(), ncols1, cpm->sol_);
  si->setDblParam(ClpObjOffset, cpm->originalOffset_);
  // Test below only needed for QP ..... but .....
  // To switch off define COIN_SLOW_PRESOLVE=0
#ifndef COIN_SLOW_PRESOLVE
#define COIN_SLOW_PRESOLVE 1
#endif
  for (int j = 0; j < ncols1; j++) {
#if COIN_SLOW_PRESOLVE
    if (cpm->hincol_[j]) {
#endif
      CoinBigIndex kcs = cpm->mcstrt_[j];
      CoinBigIndex kce = kcs + cpm->hincol_[j];
      for (CoinBigIndex k = kcs; k < kce; ++k) {
    	  cpm->link_[k] = k + 1;
      }
      cpm->link_[kce - 1] = NO_LINK;
#if COIN_SLOW_PRESOLVE
    }
#endif
  }
  {
    CoinBigIndex ml = cpm->maxlink_;
    for (CoinBigIndex k = nelemsr; k < ml; ++k)
    	cpm->link_[k] = k + 1;
    if (ml)
    	cpm->link_[ml - 1] = NO_LINK;
  }
  cpm->free_list_ = nelemsr;
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
  /*
       These are used to track the action of postsolve transforms during debugging.
     */
  CoinFillN(cpm->cdone_, ncols1, PRESENT_IN_REDUCED);
  CoinZeroN(cpm->cdone_ + ncols1, ncols0_in - ncols1);
  CoinFillN(cpm->rdone_, nrows1, PRESENT_IN_REDUCED);
  CoinZeroN(cpm->rdone_ + nrows1, cpm->nrows0_in - nrows1);
#endif

  return cpm;
}
/* This is main part of Presolve */
ClpSimplex *
ClpPresolve::gutsOfPresolvedModel(ClpSimplex *originalModel,
  double feasibilityTolerance,
  bool keepIntegers,
  int numberPasses,
  bool dropNames,
  bool doRowObjective,
  const char *prohibitedRows,
  const char *prohibitedColumns)
{
  ncols_ = originalModel->getNumCols();
  nrows_ = originalModel->getNumRows();
  nelems_ = originalModel->getNumElements();
  numberPasses_ = numberPasses;

  double maxmin = originalModel->getObjSense();
  originalModel_ = originalModel;
  delete[] originalColumn_;
  originalColumn_ = new int[ncols_];
  delete[] originalRow_;
  originalRow_ = new int[nrows_];
  // and fill in case returns early
  int i;
  for (i = 0; i < ncols_; i++)
    originalColumn_[i] = i;
  for (i = 0; i < nrows_; i++)
    originalRow_[i] = i;
  delete[] rowObjective_;
  if (doRowObjective) {
    rowObjective_ = new double[nrows_];
    memset(rowObjective_, 0, nrows_ * sizeof(double));
  } else {
    rowObjective_ = NULL;
  }

  // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
  int result = -1;

  // User may have deleted - its their responsibility
  presolvedModel_ = NULL;
  // Messages
  CoinMessages messages = originalModel->coinMessages();
  // Only go round 100 times even if integer preprocessing
  int totalPasses = 100;
  while (result == -1) {

#ifndef CLP_NO_STD
    // make new copy
    if (saveFile_ == "") {
#endif
      delete presolvedModel_;
#ifndef CLP_NO_STD
      // So won't get names
      int lengthNames = originalModel->lengthNames();
      originalModel->setLengthNames(0);
#endif
      presolvedModel_ = new ClpSimplex(*originalModel);
#ifndef CLP_NO_STD
      originalModel->setLengthNames(lengthNames);
      presolvedModel_->dropNames();
    } else {
      presolvedModel_ = originalModel;
      if (dropNames)
        presolvedModel_->dropNames();
    }
#endif

    // drop integer information if wanted
    if (!keepIntegers)
      presolvedModel_->deleteIntegerInformation();
    totalPasses--;

    double ratio = 2.0;
    if (substitution_ > 3)
      ratio = sqrt((substitution_ - 3) + 5.0);
    else if (substitution_ == 2)
      ratio = 1.5;
    CoinPresolveMatrix* prob = create_CoinPresolveMatrix(ncols_,
      maxmin,
      presolvedModel_,
      nrows_, nelems_, true, nonLinearValue_, ratio);
    if (prohibitedRows) {
      prob->setAnyProhibited();
      for (int i = 0; i < nrows_; i++) {
        if (prohibitedRows[i])
          prob->setRowProhibited(i);
      }
    }
    if (prohibitedColumns) {
      prob->setAnyProhibited();
      for (int i = 0; i < ncols_; i++) {
        if (prohibitedColumns[i])
          prob->setColProhibited(i);
      }
    }
    prob->setMaximumSubstitutionLevel(substitution_);
    if (doRowObjective)
      memset(rowObjective_, 0, nrows_ * sizeof(double));
    // See if we want statistics
    if ((presolveActions_ & 0x80000000) != 0)
      prob->statistics();
    if (doTransfer())
      transferCosts(prob);
    // make sure row solution correct
    {
      double *colels = prob->colels_;
      int *hrow = prob->hrow_;
      CoinBigIndex *mcstrt = prob->mcstrt_;
      int *hincol = prob->hincol_;
      int ncols = prob->ncols_;

      double *csol = prob->sol_;
      double *acts = prob->acts_;
      int nrows = prob->nrows_;

      int colx;

      memset(acts, 0, nrows * sizeof(double));

      for (colx = 0; colx < ncols; ++colx) {
        double solutionValue = csol[colx];
        for (CoinBigIndex i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
          int row = hrow[i];
          double coeff = colels[i];
          acts[row] += solutionValue * coeff;
        }
      }
    }

    // move across feasibility tolerance
    prob->feasibilityTolerance_ = feasibilityTolerance;

    // Do presolve
    paction_ = presolve(prob);
    // Get rid of useful arrays
    prob->deleteStuff();

    result = 0;

    bool fixInfeasibility = (prob->presolveOptions_ & 16384) != 0;
    bool hasSolution = (prob->presolveOptions_ & 32768) != 0;
    if (prob->status_ == 0 && paction_ && (!hasSolution || !fixInfeasibility)) {
      // Looks feasible but double check to see if anything slipped through
      int n = prob->ncols_;
      double *lo = prob->clo_;
      double *up = prob->cup_;
      int i;

      for (i = 0; i < n; i++) {
        if (up[i] < lo[i]) {
          if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
            // infeasible
            prob->status_ = 1;
          } else {
            up[i] = lo[i];
          }
        }
      }

      n = prob->nrows_;
      lo = prob->rlo_;
      up = prob->rup_;

      for (i = 0; i < n; i++) {
        if (up[i] < lo[i]) {
          if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
            // infeasible
            prob->status_ = 1;
          } else {
            up[i] = lo[i];
          }
        }
      }
    }
    if (prob->status_ == 0 && paction_) {
      // feasible

      update_model(*prob, presolvedModel_, nrows_, ncols_, nelems_);
      // copy status and solution
      CoinMemcpyN(prob->sol_, prob->ncols_, presolvedModel_->primalColumnSolution());
      CoinMemcpyN(prob->acts_, prob->nrows_, presolvedModel_->primalRowSolution());
      CoinMemcpyN(prob->colstat_, prob->ncols_, presolvedModel_->statusArray());
      CoinMemcpyN(prob->rowstat_, prob->nrows_, presolvedModel_->statusArray() + prob->ncols_);
      if (fixInfeasibility && hasSolution) {
        // Looks feasible but double check to see if anything slipped through
        int n = prob->ncols_;
        double *lo = prob->clo_;
        double *up = prob->cup_;
        double *rsol = prob->acts_;
        //memset(prob->acts_,0,prob->nrows_*sizeof(double));
        presolvedModel_->matrix()->times(prob->sol_, rsol);
        int i;

        for (i = 0; i < n; i++) {
          double gap = up[i] - lo[i];
          if (rsol[i] < lo[i] - feasibilityTolerance && fabs(rsol[i] - lo[i]) < 1.0e-3) {
            lo[i] = rsol[i];
            if (gap < 1.0e5)
              up[i] = lo[i] + gap;
          } else if (rsol[i] > up[i] + feasibilityTolerance && fabs(rsol[i] - up[i]) < 1.0e-3) {
            up[i] = rsol[i];
            if (gap < 1.0e5)
              lo[i] = up[i] - gap;
          }
          if (up[i] < lo[i]) {
            up[i] = lo[i];
          }
        }
      }

      int n = prob->nrows_;
      double *lo = prob->rlo_;
      double *up = prob->rup_;

      for (i = 0; i < n; i++) {
        if (up[i] < lo[i]) {
          if (up[i] < lo[i] - feasibilityTolerance && !fixInfeasibility) {
            // infeasible
            prob->status_ = 1;
          } else {
            up[i] = lo[i];
          }
        }
      }
      delete[] prob->sol_;
      delete[] prob->acts_;
      delete[] prob->colstat_;
      prob->sol_ = NULL;
      prob->acts_ = NULL;
      prob->colstat_ = NULL;

      int ncolsNow = presolvedModel_->getNumCols();
      CoinMemcpyN(prob->originalColumn_, ncolsNow, originalColumn_);
#ifndef SLIM_CLP
#ifndef NO_RTTI
      ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(originalModel->objectiveAsObject()));
#else
      ClpQuadraticObjective *quadraticObj = NULL;
      if (originalModel->objectiveAsObject()->type() == 2)
        quadraticObj = (static_cast< ClpQuadraticObjective * >(originalModel->objectiveAsObject()));
#endif
      if (quadraticObj) {
        // set up for subset
        char *mark = new char[ncols_];
        memset(mark, 0, ncols_);
        CoinPackedMatrix *quadratic = quadraticObj->quadraticObjective();
        //const int * columnQuadratic = quadratic->getIndices();
        //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
        const int *columnQuadraticLength = quadratic->getVectorLengths();
        //double * quadraticElement = quadratic->getMutableElements();
        int numberColumns = quadratic->getNumCols();
        ClpQuadraticObjective *newObj = new ClpQuadraticObjective(*quadraticObj,
          ncolsNow,
          originalColumn_);
        // and modify linear and check
        double *linear = newObj->linearObjective();
        CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
        int iColumn;
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
          if (columnQuadraticLength[iColumn])
            mark[iColumn] = 1;
        }
        // and new
        quadratic = newObj->quadraticObjective();
        columnQuadraticLength = quadratic->getVectorLengths();
        int numberColumns2 = quadratic->getNumCols();
        for (iColumn = 0; iColumn < numberColumns2; iColumn++) {
          if (columnQuadraticLength[iColumn])
            mark[originalColumn_[iColumn]] = 0;
        }
        presolvedModel_->setObjective(newObj);
        delete newObj;
        // final check
        for (iColumn = 0; iColumn < numberColumns; iColumn++)
          if (mark[iColumn])
            printf("Quadratic column %d modified - may be okay\n", iColumn);
        delete[] mark;
      }
#endif
      delete[] prob->originalColumn_;
      prob->originalColumn_ = NULL;
      int nrowsNow = presolvedModel_->getNumRows();
      CoinMemcpyN(prob->originalRow_, nrowsNow, originalRow_);
      delete[] prob->originalRow_;
      prob->originalRow_ = NULL;
#ifndef CLP_NO_STD
      if (!dropNames && originalModel->lengthNames()) {
        // Redo names
        int iRow;
        std::vector< std::string > rowNames;
        rowNames.reserve(nrowsNow);
        for (iRow = 0; iRow < nrowsNow; iRow++) {
          int kRow = originalRow_[iRow];
          rowNames.push_back(originalModel->rowName(kRow));
        }

        int iColumn;
        std::vector< std::string > columnNames;
        columnNames.reserve(ncolsNow);
        for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
          int kColumn = originalColumn_[iColumn];
          columnNames.push_back(originalModel->columnName(kColumn));
        }
        presolvedModel_->copyNames(rowNames, columnNames);
      } else {
        presolvedModel_->setLengthNames(0);
      }
#endif
      if (rowObjective_) {
        int iRow;
#ifndef NDEBUG
        int k = -1;
#endif
        int nObj = 0;
        for (iRow = 0; iRow < nrowsNow; iRow++) {
          int kRow = originalRow_[iRow];
#ifndef NDEBUG
          assert(kRow > k);
          k = kRow;
#endif
          rowObjective_[iRow] = rowObjective_[kRow];
          if (rowObjective_[iRow])
            nObj++;
        }
        if (nObj) {
          printf("%d costed slacks\n", nObj);
          presolvedModel_->setRowObjective(rowObjective_);
        }
      }
      /* now clean up integer variables.  This can modify original
               		  Don't do if dupcol added columns together */
      int i;
      const char *information = presolvedModel_->integerInformation();
      if ((prob->presolveOptions_ & 0x80000000) == 0 && information) {
        int numberChanges = 0;
        double *lower0 = originalModel_->columnLower();
        double *upper0 = originalModel_->columnUpper();
        double *lower = presolvedModel_->columnLower();
        double *upper = presolvedModel_->columnUpper();
        for (i = 0; i < ncolsNow; i++) {
          if (!information[i])
            continue;
          int iOriginal = originalColumn_[i];
          double lowerValue0 = lower0[iOriginal];
          double upperValue0 = upper0[iOriginal];
          double lowerValue = ceil(lower[i] - 1.0e-5);
          double upperValue = floor(upper[i] + 1.0e-5);
          lower[i] = lowerValue;
          upper[i] = upperValue;
          if (lowerValue > upperValue) {
            numberChanges++;
            presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
              messages)
              << iOriginal
              << lowerValue
              << upperValue
              << CoinMessageEol;
            result = 1;
          } else {
            if (lowerValue > lowerValue0 + 1.0e-8) {
              lower0[iOriginal] = lowerValue;
              numberChanges++;
            }
            if (upperValue < upperValue0 - 1.0e-8) {
              upper0[iOriginal] = upperValue;
              numberChanges++;
            }
          }
        }
        if (numberChanges) {
          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
            messages)
            << numberChanges
            << CoinMessageEol;
          if (!result && totalPasses > 0) {
            result = -1; // round again
            const CoinPresolveAction *paction = paction_;
            while (paction) {
              const CoinPresolveAction *next = paction->next;
              delete paction;
              paction = next;
            }
            paction_ = NULL;
          }
        }
      }
    } else if (prob->status_) {
      // infeasible or unbounded
      result = 1;
      // Put status in nelems_!
      nelems_ = -prob->status_;
      originalModel->setProblemStatus(prob->status_);
    } else {
      // no changes - model needs restoring after Lou's changes
#ifndef CLP_NO_STD
      if (saveFile_ == "") {
#endif
        delete presolvedModel_;
        presolvedModel_ = new ClpSimplex(*originalModel);
        // but we need to remove gaps
        ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(presolvedModel_->clpMatrix());
        if (clpMatrix) {
          clpMatrix->getPackedMatrix()->removeGaps();
        }
#ifndef CLP_NO_STD
      } else {
        presolvedModel_ = originalModel;
      }
      presolvedModel_->dropNames();
#endif

      // drop integer information if wanted
      if (!keepIntegers)
        presolvedModel_->deleteIntegerInformation();
      result = 2;
    }

    delete prob;
  }
  if (result == 0 || result == 2) {
    int nrowsAfter = presolvedModel_->getNumRows();
    int ncolsAfter = presolvedModel_->getNumCols();
    CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
    presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
      messages)
      << nrowsAfter << -(nrows_ - nrowsAfter)
      << ncolsAfter << -(ncols_ - ncolsAfter)
      << nelsAfter << -(nelems_ - nelsAfter)
      << CoinMessageEol;
  } else {
    destroyPresolve();
    if (presolvedModel_ != originalModel_)
      delete presolvedModel_;
    presolvedModel_ = NULL;
  }

  return presolvedModel_;
}

/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
*/