Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpPresolve.cpp

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