Ignore:
Timestamp:
Feb 26, 2010 12:27:59 PM (10 years ago)
Author:
mjs
Message:

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

File:
1 edited

Legend:

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

    r1464 r1525  
    4444
    4545ClpPresolve::ClpPresolve() :
    46   originalModel_(NULL),
    47   presolvedModel_(NULL),
    48   nonLinearValue_(0.0),
    49   originalColumn_(NULL),
    50   originalRow_(NULL),
    51   rowObjective_(NULL),
    52   paction_(0),   
    53   ncols_(0),
    54   nrows_(0),
    55   nelems_(0),
    56   numberPasses_(5),
    57   substitution_(3),
     46     originalModel_(NULL),
     47     presolvedModel_(NULL),
     48     nonLinearValue_(0.0),
     49     originalColumn_(NULL),
     50     originalRow_(NULL),
     51     rowObjective_(NULL),
     52     paction_(0),
     53     ncols_(0),
     54     nrows_(0),
     55     nelems_(0),
     56     numberPasses_(5),
     57     substitution_(3),
    5858#ifndef CLP_NO_STD
    59   saveFile_(""),
    60 #endif
    61   presolveActions_(0)
     59     saveFile_(""),
     60#endif
     61     presolveActions_(0)
    6262{
    6363}
     
    6565ClpPresolve::~ClpPresolve()
    6666{
    67   destroyPresolve();
     67     destroyPresolve();
    6868}
    6969// Gets rid of presolve actions (e.g.when infeasible)
    70 void 
     70void
    7171ClpPresolve::destroyPresolve()
    7272{
    73   const CoinPresolveAction *paction = paction_;
    74   while (paction) {
    75     const CoinPresolveAction *next = paction->next;
    76     delete paction;
    77     paction = next;
    78   }
    79   delete [] originalColumn_;
    80   delete [] originalRow_;
    81   paction_=NULL;
    82   originalColumn_=NULL;
    83   originalRow_=NULL;
    84   delete [] rowObjective_;
    85   rowObjective_=NULL;
    86 }
    87 
    88 /* This version of presolve returns a pointer to a new presolved 
     73     const CoinPresolveAction *paction = paction_;
     74     while (paction) {
     75          const CoinPresolveAction *next = paction->next;
     76          delete paction;
     77          paction = next;
     78     }
     79     delete [] originalColumn_;
     80     delete [] originalRow_;
     81     paction_ = NULL;
     82     originalColumn_ = NULL;
     83     originalRow_ = NULL;
     84     delete [] rowObjective_;
     85     rowObjective_ = NULL;
     86}
     87
     88/* This version of presolve returns a pointer to a new presolved
    8989   model.  NULL if infeasible
    9090*/
    91 ClpSimplex * 
     91ClpSimplex *
    9292ClpPresolve::presolvedModel(ClpSimplex & si,
    9393                            double feasibilityTolerance,
     
    9797                            bool doRowObjective)
    9898{
    99   // Check matrix
    100   if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
    101                                           1.0e20))
    102     return NULL;
    103   else
    104     return gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,dropNames,
    105                                 doRowObjective);
     99     // Check matrix
     100     if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
     101                                             1.0e20))
     102          return NULL;
     103     else
     104          return gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, dropNames,
     105                                      doRowObjective);
    106106}
    107107#ifndef CLP_NO_STD
     
    110110*/
    111111int
    112 ClpPresolve::presolvedModelToFile(ClpSimplex &si,std::string fileName,
     112ClpPresolve::presolvedModelToFile(ClpSimplex &si, std::string fileName,
    113113                                  double feasibilityTolerance,
    114114                                  bool keepIntegers,
     
    116116                                  bool doRowObjective)
    117117{
    118   // Check matrix
    119   if (!si.clpMatrix()->allElementsInRange(&si,si.getSmallElementValue(),
    120                                           1.0e20))
    121     return 2;
    122   saveFile_=fileName;
    123   si.saveModel(saveFile_.c_str());
    124   ClpSimplex * model = gutsOfPresolvedModel(&si,feasibilityTolerance,keepIntegers,numberPasses,true,
    125                                             doRowObjective);
    126   if (model==&si) {
    127     return 0;
    128   } else {
    129     si.restoreModel(saveFile_.c_str());
    130     remove(saveFile_.c_str());
    131     return 1;
    132   }
     118     // Check matrix
     119     if (!si.clpMatrix()->allElementsInRange(&si, si.getSmallElementValue(),
     120                                             1.0e20))
     121          return 2;
     122     saveFile_ = fileName;
     123     si.saveModel(saveFile_.c_str());
     124     ClpSimplex * model = gutsOfPresolvedModel(&si, feasibilityTolerance, keepIntegers, numberPasses, true,
     125                          doRowObjective);
     126     if (model == &si) {
     127          return 0;
     128     } else {
     129          si.restoreModel(saveFile_.c_str());
     130          remove(saveFile_.c_str());
     131          return 1;
     132     }
    133133}
    134134#endif
    135135// Return pointer to presolved model
    136 ClpSimplex * 
     136ClpSimplex *
    137137ClpPresolve::model() const
    138138{
    139   return presolvedModel_;
     139     return presolvedModel_;
    140140}
    141141// Return pointer to original model
    142 ClpSimplex * 
     142ClpSimplex *
    143143ClpPresolve::originalModel() const
    144144{
    145   return originalModel_;
    146 }
    147 void 
     145     return originalModel_;
     146}
     147void
    148148ClpPresolve::postsolve(bool updateStatus)
    149149{
    150   // Return at once if no presolved model
    151   if (!presolvedModel_)
    152     return;
    153   // Messages
    154   CoinMessages messages = originalModel_->coinMessages();
    155   if (!presolvedModel_->isProvenOptimal()) {
    156     presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
    157                                              messages)
    158                                                <<CoinMessageEol;
    159   }
    160 
    161   // this is the size of the original problem
    162   const int ncols0  = ncols_;
    163   const int nrows0  = nrows_;
    164   const CoinBigIndex nelems0 = nelems_;
    165  
    166   // this is the reduced problem
    167   int ncols = presolvedModel_->getNumCols();
    168   int nrows = presolvedModel_->getNumRows();
    169 
    170   double * acts=NULL;
    171   double * sol =NULL;
    172   unsigned char * rowstat=NULL;
    173   unsigned char * colstat = NULL;
     150     // Return at once if no presolved model
     151     if (!presolvedModel_)
     152          return;
     153     // Messages
     154     CoinMessages messages = originalModel_->coinMessages();
     155     if (!presolvedModel_->isProvenOptimal()) {
     156          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NONOPTIMAL,
     157                    messages)
     158                    << CoinMessageEol;
     159     }
     160
     161     // this is the size of the original problem
     162     const int ncols0  = ncols_;
     163     const int nrows0  = nrows_;
     164     const CoinBigIndex nelems0 = nelems_;
     165
     166     // this is the reduced problem
     167     int ncols = presolvedModel_->getNumCols();
     168     int nrows = presolvedModel_->getNumRows();
     169
     170     double * acts = NULL;
     171     double * sol = NULL;
     172     unsigned char * rowstat = NULL;
     173     unsigned char * colstat = NULL;
    174174#ifndef CLP_NO_STD
    175   if (saveFile_=="") {
    176 #endif
    177     // reality check
    178     assert(ncols0==originalModel_->getNumCols());
    179     assert(nrows0==originalModel_->getNumRows());
    180     acts = originalModel_->primalRowSolution();
    181     sol  = originalModel_->primalColumnSolution();
    182     if (updateStatus) {
    183       // postsolve does not know about fixed
    184       int i;
    185       for (i=0;i<nrows+ncols;i++) {
    186         if (presolvedModel_->getColumnStatus(i)==ClpSimplex::isFixed)
    187           presolvedModel_->setColumnStatus(i,ClpSimplex::atLowerBound);
    188       }
    189       unsigned char *status = originalModel_->statusArray();
    190       if (!status) {
    191         originalModel_->createStatus();
    192         status = originalModel_->statusArray();
    193       }
    194       rowstat = status + ncols0;
    195       colstat = status;
    196       CoinMemcpyN( presolvedModel_->statusArray(), ncols,colstat);
    197       CoinMemcpyN( presolvedModel_->statusArray()+ncols, nrows,rowstat);
    198     }
     175     if (saveFile_ == "") {
     176#endif
     177          // reality check
     178          assert(ncols0 == originalModel_->getNumCols());
     179          assert(nrows0 == originalModel_->getNumRows());
     180          acts = originalModel_->primalRowSolution();
     181          sol  = originalModel_->primalColumnSolution();
     182          if (updateStatus) {
     183               // postsolve does not know about fixed
     184               int i;
     185               for (i = 0; i < nrows + ncols; i++) {
     186                    if (presolvedModel_->getColumnStatus(i) == ClpSimplex::isFixed)
     187                         presolvedModel_->setColumnStatus(i, ClpSimplex::atLowerBound);
     188               }
     189               unsigned char *status = originalModel_->statusArray();
     190               if (!status) {
     191                    originalModel_->createStatus();
     192                    status = originalModel_->statusArray();
     193               }
     194               rowstat = status + ncols0;
     195               colstat = status;
     196               CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
     197               CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
     198          }
    199199#ifndef CLP_NO_STD
    200   } else {
    201     // from file
    202     acts = new double[nrows0];
    203     sol  = new double[ncols0];
    204     CoinZeroN(acts,nrows0);
    205     CoinZeroN(sol,ncols0);
    206     if (updateStatus) {
    207       unsigned char *status = new unsigned char [nrows0+ncols0];
    208       rowstat = status + ncols0;
    209       colstat = status;
    210       CoinMemcpyN( presolvedModel_->statusArray(), ncols,colstat);
    211       CoinMemcpyN( presolvedModel_->statusArray()+ncols, nrows,rowstat);
    212     }
    213   }
    214 #endif
    215 
    216   // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
    217   // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
    218   // cause duplicate free. In case where saveFile == "", as best I can see
    219   // arrays are owned by originalModel_. fix is to
    220   // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
    221   CoinPostsolveMatrix prob(presolvedModel_,
    222                        ncols0,
    223                        nrows0,
    224                        nelems0,
    225                        presolvedModel_->getObjSense(),
    226                        // end prepost
    227                        
    228                        sol, acts,
    229                        colstat, rowstat);
    230    
    231   postsolve(prob);
     200     } else {
     201          // from file
     202          acts = new double[nrows0];
     203          sol  = new double[ncols0];
     204          CoinZeroN(acts, nrows0);
     205          CoinZeroN(sol, ncols0);
     206          if (updateStatus) {
     207               unsigned char *status = new unsigned char [nrows0+ncols0];
     208               rowstat = status + ncols0;
     209               colstat = status;
     210               CoinMemcpyN( presolvedModel_->statusArray(), ncols, colstat);
     211               CoinMemcpyN( presolvedModel_->statusArray() + ncols, nrows, rowstat);
     212          }
     213     }
     214#endif
     215
     216     // CoinPostsolveMatrix object assumes ownership of sol, acts, colstat;
     217     // will be deleted by ~CoinPostsolveMatrix. delete[] operations below
     218     // cause duplicate free. In case where saveFile == "", as best I can see
     219     // arrays are owned by originalModel_. fix is to
     220     // clear fields in prob to avoid delete[] in ~CoinPostsolveMatrix.
     221     CoinPostsolveMatrix prob(presolvedModel_,
     222                              ncols0,
     223                              nrows0,
     224                              nelems0,
     225                              presolvedModel_->getObjSense(),
     226                              // end prepost
     227
     228                              sol, acts,
     229                              colstat, rowstat);
     230
     231     postsolve(prob);
    232232
    233233#ifndef CLP_NO_STD
    234   if (saveFile_!="") {
    235     // From file
    236     assert (originalModel_==presolvedModel_);
    237     originalModel_->restoreModel(saveFile_.c_str());
    238     remove(saveFile_.c_str());
    239     CoinMemcpyN(acts,nrows0,originalModel_->primalRowSolution());
    240     // delete [] acts;
    241     CoinMemcpyN(sol,ncols0,originalModel_->primalColumnSolution());
    242     // delete [] sol;
    243     if (updateStatus) {
    244       CoinMemcpyN(colstat,nrows0+ncols0,originalModel_->statusArray());
    245       // delete [] colstat;
    246     }
    247   } else {
    248 #endif
    249     prob.sol_ = 0 ;
    250     prob.acts_ = 0 ;
    251     prob.colstat_ = 0 ;
     234     if (saveFile_ != "") {
     235          // From file
     236          assert (originalModel_ == presolvedModel_);
     237          originalModel_->restoreModel(saveFile_.c_str());
     238          remove(saveFile_.c_str());
     239          CoinMemcpyN(acts, nrows0, originalModel_->primalRowSolution());
     240          // delete [] acts;
     241          CoinMemcpyN(sol, ncols0, originalModel_->primalColumnSolution());
     242          // delete [] sol;
     243          if (updateStatus) {
     244               CoinMemcpyN(colstat, nrows0 + ncols0, originalModel_->statusArray());
     245               // delete [] colstat;
     246          }
     247     } else {
     248#endif
     249          prob.sol_ = 0 ;
     250          prob.acts_ = 0 ;
     251          prob.colstat_ = 0 ;
    252252#ifndef CLP_NO_STD
    253   }
    254 #endif
    255   // put back duals
    256   CoinMemcpyN(prob.rowduals_,   nrows_,originalModel_->dualRowSolution());
    257   double maxmin = originalModel_->getObjSense();
    258   if (maxmin<0.0) {
    259     // swap signs
    260     int i;
    261     double * pi = originalModel_->dualRowSolution();
    262     for (i=0;i<nrows_;i++)
    263       pi[i] = -pi[i];
    264   }
    265   // Now check solution
    266   double offset;
    267   CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
    268                                                             originalModel_->primalColumnSolution(),offset,true),
    269               ncols_,originalModel_->dualColumnSolution());
    270   originalModel_->transposeTimes(-1.0,
    271                                 originalModel_->dualRowSolution(),
    272                                 originalModel_->dualColumnSolution());
    273   memset(originalModel_->primalRowSolution(),0,nrows_*sizeof(double));
    274   originalModel_->times(1.0,originalModel_->primalColumnSolution(),
    275                         originalModel_->primalRowSolution());
    276   originalModel_->checkSolutionInternal(); 
    277   if (originalModel_->sumDualInfeasibilities()>1.0e-1) {
    278     // See if we can fix easily
    279     static_cast<ClpSimplexOther *> (originalModel_)->cleanupAfterPostsolve();
    280   }
    281   // Messages
    282   presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
    283                                             messages)
    284                                               <<originalModel_->objectiveValue()
    285                                               <<originalModel_->sumDualInfeasibilities()
    286                                               <<originalModel_->numberDualInfeasibilities()
    287                                               <<originalModel_->sumPrimalInfeasibilities()
    288                                               <<originalModel_->numberPrimalInfeasibilities()
    289                                                <<CoinMessageEol;
    290  
    291   //originalModel_->objectiveValue_=objectiveValue_;
    292   originalModel_->setNumberIterations(presolvedModel_->numberIterations());
    293   if (!presolvedModel_->status()) {
    294     if (!originalModel_->numberDualInfeasibilities()&&
    295         !originalModel_->numberPrimalInfeasibilities()) {
    296       originalModel_->setProblemStatus( 0);
    297     } else {
    298       originalModel_->setProblemStatus( -1);
    299       // Say not optimal after presolve
    300       originalModel_->setSecondaryStatus(7);
    301       presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
    302                                             messages)
    303                                               <<CoinMessageEol;
    304     }
    305   } else {
    306     originalModel_->setProblemStatus( presolvedModel_->status());
    307   }
     253     }
     254#endif
     255     // put back duals
     256     CoinMemcpyN(prob.rowduals_,        nrows_, originalModel_->dualRowSolution());
     257     double maxmin = originalModel_->getObjSense();
     258     if (maxmin < 0.0) {
     259          // swap signs
     260          int i;
     261          double * pi = originalModel_->dualRowSolution();
     262          for (i = 0; i < nrows_; i++)
     263               pi[i] = -pi[i];
     264     }
     265     // Now check solution
     266     double offset;
     267     CoinMemcpyN(originalModel_->objectiveAsObject()->gradient(originalModel_,
     268                 originalModel_->primalColumnSolution(), offset, true),
     269                 ncols_, originalModel_->dualColumnSolution());
     270     originalModel_->transposeTimes(-1.0,
     271                                    originalModel_->dualRowSolution(),
     272                                    originalModel_->dualColumnSolution());
     273     memset(originalModel_->primalRowSolution(), 0, nrows_ * sizeof(double));
     274     originalModel_->times(1.0, originalModel_->primalColumnSolution(),
     275                           originalModel_->primalRowSolution());
     276     originalModel_->checkSolutionInternal();
     277     if (originalModel_->sumDualInfeasibilities() > 1.0e-1) {
     278          // See if we can fix easily
     279          static_cast<ClpSimplexOther *> (originalModel_)->cleanupAfterPostsolve();
     280     }
     281     // Messages
     282     presolvedModel_->messageHandler()->message(COIN_PRESOLVE_POSTSOLVE,
     283               messages)
     284               << originalModel_->objectiveValue()
     285               << originalModel_->sumDualInfeasibilities()
     286               << originalModel_->numberDualInfeasibilities()
     287               << originalModel_->sumPrimalInfeasibilities()
     288               << originalModel_->numberPrimalInfeasibilities()
     289               << CoinMessageEol;
     290
     291     //originalModel_->objectiveValue_=objectiveValue_;
     292     originalModel_->setNumberIterations(presolvedModel_->numberIterations());
     293     if (!presolvedModel_->status()) {
     294          if (!originalModel_->numberDualInfeasibilities() &&
     295                    !originalModel_->numberPrimalInfeasibilities()) {
     296               originalModel_->setProblemStatus( 0);
     297          } else {
     298               originalModel_->setProblemStatus( -1);
     299               // Say not optimal after presolve
     300               originalModel_->setSecondaryStatus(7);
     301               presolvedModel_->messageHandler()->message(COIN_PRESOLVE_NEEDS_CLEANING,
     302                         messages)
     303                         << CoinMessageEol;
     304          }
     305     } else {
     306          originalModel_->setProblemStatus( presolvedModel_->status());
     307     }
    308308#ifndef CLP_NO_STD
    309   if (saveFile_!="")
    310     presolvedModel_=NULL;
     309     if (saveFile_ != "")
     310          presolvedModel_ = NULL;
    311311#endif
    312312}
    313313
    314314// return pointer to original columns
    315 const int * 
     315const int *
    316316ClpPresolve::originalColumns() const
    317317{
    318   return originalColumn_;
     318     return originalColumn_;
    319319}
    320320// return pointer to original rows
    321 const int * 
     321const int *
    322322ClpPresolve::originalRows() const
    323323{
    324   return originalRow_;
     324     return originalRow_;
    325325}
    326326// Set pointer to original model
    327 void 
     327void
    328328ClpPresolve::setOriginalModel(ClpSimplex * model)
    329329{
    330   originalModel_=model;
     330     originalModel_ = model;
    331331}
    332332#if 0
     
    335335static int ATOI(const char *name)
    336336{
    337  return true;
     337     return true;
    338338#if     PRESOLVE_DEBUG || PRESOLVE_SUMMARY
    339   if (getenv(name)) {
    340     int val = atoi(getenv(name));
    341     printf("%s = %d\n", name, val);
    342     return (val);
    343   } else {
    344     if (strcmp(name,"off"))
    345       return (true);
    346     else
    347       return (false);
    348   }
     339     if (getenv(name)) {
     340          int val = atoi(getenv(name));
     341          printf("%s = %d\n", name, val);
     342          return (val);
     343     } else {
     344          if (strcmp(name, "off"))
     345               return (true);
     346          else
     347               return (false);
     348     }
    349349#else
    350   return (true);
     350     return (true);
    351351#endif
    352352}
     
    354354//#define PRESOLVE_DEBUG 1
    355355#if PRESOLVE_DEBUG
    356 void check_sol(CoinPresolveMatrix *prob,double tol)
    357 {
    358   double *colels        = prob->colels_;
    359   int *hrow             = prob->hrow_;
    360   int *mcstrt           = prob->mcstrt_;
    361   int *hincol           = prob->hincol_;
    362   int *hinrow           = prob->hinrow_;
    363   int ncols             = prob->ncols_;
    364 
    365 
    366   double * csol = prob->sol_;
    367   double * acts = prob->acts_;
    368   double * clo = prob->clo_;
    369   double * cup = prob->cup_;
    370   int nrows = prob->nrows_;
    371   double * rlo = prob->rlo_;
    372   double * rup = prob->rup_;
    373 
    374   int colx;
    375 
    376   double * rsol = new double[nrows];
    377   memset(rsol,0,nrows*sizeof(double));
    378 
    379   for (colx = 0; colx < ncols; ++colx) {
    380     if (1) {
    381       CoinBigIndex k = mcstrt[colx];
    382       int nx = hincol[colx];
    383       double solutionValue = csol[colx];
    384       for (int i=0; i<nx; ++i) {
    385         int row = hrow[k];
    386         double coeff = colels[k];
    387         k++;
    388         rsol[row] += solutionValue*coeff;
    389       }
    390       if (csol[colx]<clo[colx]-tol) {
    391         printf("low CSOL:  %d  - %g %g %g\n",
    392                    colx, clo[colx], csol[colx], cup[colx]);
    393       } else if (csol[colx]>cup[colx]+tol) {
    394         printf("high CSOL:  %d  - %g %g %g\n",
    395                    colx, clo[colx], csol[colx], cup[colx]);
    396       }
    397     }
    398   }
    399   int rowx;
    400   for (rowx = 0; rowx < nrows; ++rowx) {
    401     if (hinrow[rowx]) {
    402       if (fabs(rsol[rowx]-acts[rowx])>tol)
    403         printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
    404                    rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
    405       if (rsol[rowx]<rlo[rowx]-tol) {
    406         printf("low RSOL:  %d - %g %g %g\n",
    407                    rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
    408       } else if (rsol[rowx]>rup[rowx]+tol ) {
    409         printf("high RSOL:  %d - %g %g %g\n",
    410                    rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
    411       }
    412     }
    413   }
    414   delete [] rsol;
     356void check_sol(CoinPresolveMatrix *prob, double tol)
     357{
     358     double *colels     = prob->colels_;
     359     int *hrow          = prob->hrow_;
     360     int *mcstrt                = prob->mcstrt_;
     361     int *hincol                = prob->hincol_;
     362     int *hinrow                = prob->hinrow_;
     363     int ncols          = prob->ncols_;
     364
     365
     366     double * csol = prob->sol_;
     367     double * acts = prob->acts_;
     368     double * clo = prob->clo_;
     369     double * cup = prob->cup_;
     370     int nrows = prob->nrows_;
     371     double * rlo = prob->rlo_;
     372     double * rup = prob->rup_;
     373
     374     int colx;
     375
     376     double * rsol = new double[nrows];
     377     memset(rsol, 0, nrows * sizeof(double));
     378
     379     for (colx = 0; colx < ncols; ++colx) {
     380          if (1) {
     381               CoinBigIndex k = mcstrt[colx];
     382               int nx = hincol[colx];
     383               double solutionValue = csol[colx];
     384               for (int i = 0; i < nx; ++i) {
     385                    int row = hrow[k];
     386                    double coeff = colels[k];
     387                    k++;
     388                    rsol[row] += solutionValue * coeff;
     389               }
     390               if (csol[colx] < clo[colx] - tol) {
     391                    printf("low CSOL:  %d  - %g %g %g\n",
     392                           colx, clo[colx], csol[colx], cup[colx]);
     393               } else if (csol[colx] > cup[colx] + tol) {
     394                    printf("high CSOL:  %d  - %g %g %g\n",
     395                           colx, clo[colx], csol[colx], cup[colx]);
     396               }
     397          }
     398     }
     399     int rowx;
     400     for (rowx = 0; rowx < nrows; ++rowx) {
     401          if (hinrow[rowx]) {
     402               if (fabs(rsol[rowx] - acts[rowx]) > tol)
     403                    printf("inacc RSOL:  %d - %g %g (acts_ %g) %g\n",
     404                           rowx,  rlo[rowx], rsol[rowx], acts[rowx], rup[rowx]);
     405               if (rsol[rowx] < rlo[rowx] - tol) {
     406                    printf("low RSOL:  %d - %g %g %g\n",
     407                           rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
     408               } else if (rsol[rowx] > rup[rowx] + tol ) {
     409                    printf("high RSOL:  %d - %g %g %g\n",
     410                           rowx,  rlo[rowx], rsol[rowx], rup[rowx]);
     411               }
     412          }
     413     }
     414     delete [] rsol;
    415415}
    416416#endif
     
    420420const CoinPresolveAction *ClpPresolve::presolve(CoinPresolveMatrix *prob)
    421421{
    422   // Messages
    423   CoinMessages messages = CoinMessage(prob->messages().language());
    424   paction_ = 0;
    425 
    426   prob->status_=0; // say feasible
    427   paction_ = make_fixed(prob, paction_);
    428   // if integers then switch off dual stuff
    429   // later just do individually
    430   bool doDualStuff = (presolvedModel_->integerInformation()==NULL);
    431   // but allow in some cases
    432   if ((presolveActions_&512)!=0)
    433     doDualStuff=true;
    434   if (prob->anyProhibited())
    435     doDualStuff=false;
    436   if (!doDual())
    437     doDualStuff=false;
     422     // Messages
     423     CoinMessages messages = CoinMessage(prob->messages().language());
     424     paction_ = 0;
     425
     426     prob->status_ = 0; // say feasible
     427     paction_ = make_fixed(prob, paction_);
     428     // if integers then switch off dual stuff
     429     // later just do individually
     430     bool doDualStuff = (presolvedModel_->integerInformation() == NULL);
     431     // but allow in some cases
     432     if ((presolveActions_ & 512) != 0)
     433          doDualStuff = true;
     434     if (prob->anyProhibited())
     435          doDualStuff = false;
     436     if (!doDual())
     437          doDualStuff = false;
    438438#if     PRESOLVE_CONSISTENCY
    439439//  presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, prob->nrows_);
    440     presolve_links_ok(prob,false,true) ;
    441 #endif
    442 
    443   if (!prob->status_) {
    444     bool slackSingleton = doSingletonColumn();
    445     slackSingleton = true;
    446     const bool slackd = doSingleton();
    447     const bool doubleton = doDoubleton();
    448     const bool tripleton = doTripleton();
    449     //#define NO_FORCING
     440     presolve_links_ok(prob, false, true) ;
     441#endif
     442
     443     if (!prob->status_) {
     444          bool slackSingleton = doSingletonColumn();
     445          slackSingleton = true;
     446          const bool slackd = doSingleton();
     447          const bool doubleton = doDoubleton();
     448          const bool tripleton = doTripleton();
     449          //#define NO_FORCING
    450450#ifndef NO_FORCING
    451     const bool forcing = doForcing();
    452 #endif
    453     const bool ifree = doImpliedFree();
    454     const bool zerocost = doTighten();
    455     const bool dupcol = doDupcol();
    456     const bool duprow = doDuprow();
    457     const bool dual = doDualStuff;
    458    
    459     // some things are expensive so just do once (normally)
    460 
    461     int i;
    462     // say look at all
    463     if (!prob->anyProhibited()) {
    464       for (i=0;i<nrows_;i++)
    465         prob->rowsToDo_[i]=i;
    466       prob->numberRowsToDo_=nrows_;
    467       for (i=0;i<ncols_;i++)
    468         prob->colsToDo_[i]=i;
    469       prob->numberColsToDo_=ncols_;
    470     } else {
    471       // some stuff must be left alone
    472       prob->numberRowsToDo_=0;
    473       for (i=0;i<nrows_;i++)
    474         if (!prob->rowProhibited(i))
    475             prob->rowsToDo_[prob->numberRowsToDo_++]=i;
    476       prob->numberColsToDo_=0;
    477       for (i=0;i<ncols_;i++)
    478         if (!prob->colProhibited(i))
    479             prob->colsToDo_[prob->numberColsToDo_++]=i;
    480     }
    481 
    482 
    483     int iLoop;
     451          const bool forcing = doForcing();
     452#endif
     453          const bool ifree = doImpliedFree();
     454          const bool zerocost = doTighten();
     455          const bool dupcol = doDupcol();
     456          const bool duprow = doDuprow();
     457          const bool dual = doDualStuff;
     458
     459          // some things are expensive so just do once (normally)
     460
     461          int i;
     462          // say look at all
     463          if (!prob->anyProhibited()) {
     464               for (i = 0; i < nrows_; i++)
     465                    prob->rowsToDo_[i] = i;
     466               prob->numberRowsToDo_ = nrows_;
     467               for (i = 0; i < ncols_; i++)
     468                    prob->colsToDo_[i] = i;
     469               prob->numberColsToDo_ = ncols_;
     470          } else {
     471               // some stuff must be left alone
     472               prob->numberRowsToDo_ = 0;
     473               for (i = 0; i < nrows_; i++)
     474                    if (!prob->rowProhibited(i))
     475                         prob->rowsToDo_[prob->numberRowsToDo_++] = i;
     476               prob->numberColsToDo_ = 0;
     477               for (i = 0; i < ncols_; i++)
     478                    if (!prob->colProhibited(i))
     479                         prob->colsToDo_[prob->numberColsToDo_++] = i;
     480          }
     481
     482
     483          int iLoop;
    484484#if     PRESOLVE_DEBUG
    485     check_sol(prob,1.0e0);
    486 #endif
    487     if (dupcol) {
    488       // maybe allow integer columns to be checked
    489       if ((presolveActions_&512)!=0)
    490         prob->setPresolveOptions(prob->presolveOptions()|1);
    491       paction_ = dupcol_action::presolve(prob, paction_);
    492     }
    493     if (duprow) {
    494       paction_ = duprow_action::presolve(prob, paction_);
    495     }
    496     if (doGubrow()) {
    497       paction_ = gubrow_action::presolve(prob, paction_);
    498     }
    499 
    500     if ((presolveActions_&16384)!=0)
    501       prob->setPresolveOptions(prob->presolveOptions()|16384);
    502     // Check number rows dropped
    503     int lastDropped=0;
    504     prob->pass_=0;
    505     for (iLoop=0;iLoop<numberPasses_;iLoop++) {
    506     // See if we want statistics
    507     if ((presolveActions_&0x80000000)!=0)
    508       printf("Starting major pass %d after %g seconds\n",iLoop+1,CoinCpuTime()-prob->startTime_);
     485          check_sol(prob, 1.0e0);
     486#endif
     487          if (dupcol) {
     488               // maybe allow integer columns to be checked
     489               if ((presolveActions_ & 512) != 0)
     490                    prob->setPresolveOptions(prob->presolveOptions() | 1);
     491               paction_ = dupcol_action::presolve(prob, paction_);
     492          }
     493          if (duprow) {
     494               paction_ = duprow_action::presolve(prob, paction_);
     495          }
     496          if (doGubrow()) {
     497               paction_ = gubrow_action::presolve(prob, paction_);
     498          }
     499
     500          if ((presolveActions_ & 16384) != 0)
     501               prob->setPresolveOptions(prob->presolveOptions() | 16384);
     502          // Check number rows dropped
     503          int lastDropped = 0;
     504          prob->pass_ = 0;
     505          for (iLoop = 0; iLoop < numberPasses_; iLoop++) {
     506               // See if we want statistics
     507               if ((presolveActions_ & 0x80000000) != 0)
     508                    printf("Starting major pass %d after %g seconds\n", iLoop + 1, CoinCpuTime() - prob->startTime_);
    509509#ifdef PRESOLVE_SUMMARY
    510       printf("Starting major pass %d\n",iLoop+1);
    511 #endif
    512       const CoinPresolveAction * const paction0 = paction_;
    513       // look for substitutions with no fill
    514       //#define IMPLIED 3
     510               printf("Starting major pass %d\n", iLoop + 1);
     511#endif
     512               const CoinPresolveAction * const paction0 = paction_;
     513               // look for substitutions with no fill
     514               //#define IMPLIED 3
    515515#ifdef IMPLIED
    516       int fill_level=3;
     516               int fill_level = 3;
    517517#define IMPLIED2 99
    518518#if IMPLIED!=3
    519519#if IMPLIED>2&&IMPLIED<11
    520       fill_level=IMPLIED;
    521       printf("** fill_level == %d !\n",fill_level);
     520               fill_level = IMPLIED;
     521               printf("** fill_level == %d !\n", fill_level);
    522522#endif
    523523#if IMPLIED>11&&IMPLIED<21
    524       fill_level=-(IMPLIED-10);
    525       printf("** fill_level == %d !\n",fill_level);
     524               fill_level = -(IMPLIED - 10);
     525               printf("** fill_level == %d !\n", fill_level);
    526526#endif
    527527#endif
    528528#else
    529       int fill_level=2;
    530 #endif
    531       int whichPass=0;
    532       while (1) {
    533         whichPass++;
    534         prob->pass_++;
    535         const CoinPresolveAction * const paction1 = paction_;
    536 
    537         if (slackd) {
    538           bool notFinished = true;
    539           while (notFinished)
    540             paction_ = slack_doubleton_action::presolve(prob, paction_,
    541                                                         notFinished);
    542           if (prob->status_)
    543             break;
    544         }
    545         if (dual&&whichPass==1) {
    546           // this can also make E rows so do one bit here
    547           paction_ = remove_dual_action::presolve(prob, paction_);
    548           if (prob->status_)
    549             break;
    550         }
    551 
    552         if (doubleton) {
    553           paction_ = doubleton_action::presolve(prob, paction_);
    554           if (prob->status_)
    555             break;
    556         }
    557         if (tripleton) {
    558           paction_ = tripleton_action::presolve(prob, paction_);
    559           if (prob->status_)
    560             break;
    561         }
    562 
    563         if (zerocost) {
    564           paction_ = do_tighten_action::presolve(prob, paction_);
    565           if (prob->status_)
    566             break;
    567         }
     529               int fill_level = 2;
     530#endif
     531               int whichPass = 0;
     532               while (1) {
     533                    whichPass++;
     534                    prob->pass_++;
     535                    const CoinPresolveAction * const paction1 = paction_;
     536
     537                    if (slackd) {
     538                         bool notFinished = true;
     539                         while (notFinished)
     540                              paction_ = slack_doubleton_action::presolve(prob, paction_,
     541                                         notFinished);
     542                         if (prob->status_)
     543                              break;
     544                    }
     545                    if (dual && whichPass == 1) {
     546                         // this can also make E rows so do one bit here
     547                         paction_ = remove_dual_action::presolve(prob, paction_);
     548                         if (prob->status_)
     549                              break;
     550                    }
     551
     552                    if (doubleton) {
     553                         paction_ = doubleton_action::presolve(prob, paction_);
     554                         if (prob->status_)
     555                              break;
     556                    }
     557                    if (tripleton) {
     558                         paction_ = tripleton_action::presolve(prob, paction_);
     559                         if (prob->status_)
     560                              break;
     561                    }
     562
     563                    if (zerocost) {
     564                         paction_ = do_tighten_action::presolve(prob, paction_);
     565                         if (prob->status_)
     566                              break;
     567                    }
    568568#ifndef NO_FORCING
    569         if (forcing) {
    570           paction_ = forcing_constraint_action::presolve(prob, paction_);
    571           if (prob->status_)
    572             break;
    573         }
    574 #endif
    575 
    576         if (ifree&&(whichPass%5)==1) {
    577           paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    578         if (prob->status_)
    579           break;
    580         }
     569                    if (forcing) {
     570                         paction_ = forcing_constraint_action::presolve(prob, paction_);
     571                         if (prob->status_)
     572                              break;
     573                    }
     574#endif
     575
     576                    if (ifree && (whichPass % 5) == 1) {
     577                         paction_ = implied_free_action::presolve(prob, paction_, fill_level);
     578                         if (prob->status_)
     579                              break;
     580                    }
    581581
    582582#if     PRESOLVE_DEBUG
    583         check_sol(prob,1.0e0);
     583                    check_sol(prob, 1.0e0);
    584584#endif
    585585
    586586#if     PRESOLVE_CONSISTENCY
    587 //      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_, 
     587//      presolve_links_ok(prob->rlink_, prob->mrstrt_, prob->hinrow_,
    588588//                        prob->nrows_);
    589         presolve_links_ok(prob,false,true) ;
     589                    presolve_links_ok(prob, false, true) ;
    590590#endif
    591591
    592592//#if   PRESOLVE_DEBUG
    593 //      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_, 
     593//      presolve_no_zeros(prob->mcstrt_, prob->colels_, prob->hincol_,
    594594//                        prob->ncols_);
    595595//#endif
     
    598598//#endif
    599599#if     PRESOLVE_CONSISTENCY
    600         presolve_no_zeros(prob,true,false) ;
    601         presolve_consistent(prob,true) ;
    602 #endif
    603 
    604          
    605         // set up for next pass
    606         // later do faster if many changes i.e. memset and memcpy
    607         prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
    608         //int kcheck;
    609         //bool found=false;
    610         //kcheck=-1;
    611         for (i=0;i<prob->numberNextRowsToDo_;i++) {
    612           int index = prob->nextRowsToDo_[i];
    613           prob->unsetRowChanged(index);
    614           prob->rowsToDo_[i] = index;
    615           //if (index==kcheck) {
    616           //printf("row %d on list after pass %d\n",kcheck,
    617           //   whichPass);
    618           //found=true;
    619           //}
    620         }
    621         //if (!found&&kcheck>=0)
    622         //prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
    623         prob->numberNextRowsToDo_=0;
    624         prob->numberColsToDo_ = prob->numberNextColsToDo_;
    625         //kcheck=-1;
    626         //found=false;
    627         for (i=0;i<prob->numberNextColsToDo_;i++) {
    628           int index = prob->nextColsToDo_[i];
    629           prob->unsetColChanged(index);
    630           prob->colsToDo_[i] = index;
    631           //if (index==kcheck) {
    632           //printf("col %d on list after pass %d\n",kcheck,
    633           //   whichPass);
    634           //found=true;
    635           //}
    636         }
    637         //if (!found&&kcheck>=0)
    638         //prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
    639         prob->numberNextColsToDo_=0;
    640         if (paction_ == paction1&&fill_level>0)
    641           break;
    642       }
    643       // say look at all
    644       int i;
    645       if (!prob->anyProhibited()) {
    646         for (i=0;i<nrows_;i++)
    647           prob->rowsToDo_[i]=i;
    648         prob->numberRowsToDo_=nrows_;
    649         for (i=0;i<ncols_;i++)
    650           prob->colsToDo_[i]=i;
    651         prob->numberColsToDo_=ncols_;
    652       } else {
    653         // some stuff must be left alone
    654         prob->numberRowsToDo_=0;
    655         for (i=0;i<nrows_;i++)
    656           if (!prob->rowProhibited(i))
    657             prob->rowsToDo_[prob->numberRowsToDo_++]=i;
    658         prob->numberColsToDo_=0;
    659         for (i=0;i<ncols_;i++)
    660           if (!prob->colProhibited(i))
    661             prob->colsToDo_[prob->numberColsToDo_++]=i;
    662       }
    663       // now expensive things
    664       // this caused world.mps to run into numerical difficulties
     600                    presolve_no_zeros(prob, true, false) ;
     601                    presolve_consistent(prob, true) ;
     602#endif
     603
     604
     605                    // set up for next pass
     606                    // later do faster if many changes i.e. memset and memcpy
     607                    prob->numberRowsToDo_ = prob->numberNextRowsToDo_;
     608                    //int kcheck;
     609                    //bool found=false;
     610                    //kcheck=-1;
     611                    for (i = 0; i < prob->numberNextRowsToDo_; i++) {
     612                         int index = prob->nextRowsToDo_[i];
     613                         prob->unsetRowChanged(index);
     614                         prob->rowsToDo_[i] = index;
     615                         //if (index==kcheck) {
     616                         //printf("row %d on list after pass %d\n",kcheck,
     617                         //   whichPass);
     618                         //found=true;
     619                         //}
     620                    }
     621                    //if (!found&&kcheck>=0)
     622                    //prob->rowsToDo_[prob->numberRowsToDo_++]=kcheck;
     623                    prob->numberNextRowsToDo_ = 0;
     624                    prob->numberColsToDo_ = prob->numberNextColsToDo_;
     625                    //kcheck=-1;
     626                    //found=false;
     627                    for (i = 0; i < prob->numberNextColsToDo_; i++) {
     628                         int index = prob->nextColsToDo_[i];
     629                         prob->unsetColChanged(index);
     630                         prob->colsToDo_[i] = index;
     631                         //if (index==kcheck) {
     632                         //printf("col %d on list after pass %d\n",kcheck,
     633                         //   whichPass);
     634                         //found=true;
     635                         //}
     636                    }
     637                    //if (!found&&kcheck>=0)
     638                    //prob->colsToDo_[prob->numberColsToDo_++]=kcheck;
     639                    prob->numberNextColsToDo_ = 0;
     640                    if (paction_ == paction1 && fill_level > 0)
     641                         break;
     642               }
     643               // say look at all
     644               int i;
     645               if (!prob->anyProhibited()) {
     646                    for (i = 0; i < nrows_; i++)
     647                         prob->rowsToDo_[i] = i;
     648                    prob->numberRowsToDo_ = nrows_;
     649                    for (i = 0; i < ncols_; i++)
     650                         prob->colsToDo_[i] = i;
     651                    prob->numberColsToDo_ = ncols_;
     652               } else {
     653                    // some stuff must be left alone
     654                    prob->numberRowsToDo_ = 0;
     655                    for (i = 0; i < nrows_; i++)
     656                         if (!prob->rowProhibited(i))
     657                              prob->rowsToDo_[prob->numberRowsToDo_++] = i;
     658                    prob->numberColsToDo_ = 0;
     659                    for (i = 0; i < ncols_; i++)
     660                         if (!prob->colProhibited(i))
     661                              prob->colsToDo_[prob->numberColsToDo_++] = i;
     662               }
     663               // now expensive things
     664               // this caused world.mps to run into numerical difficulties
    665665#ifdef PRESOLVE_SUMMARY
    666       printf("Starting expensive\n");
    667 #endif
    668 
    669       if (dual) {
    670         int itry;
    671         for (itry=0;itry<5;itry++) {
    672           paction_ = remove_dual_action::presolve(prob, paction_);
    673           if (prob->status_)
    674             break;
    675           const CoinPresolveAction * const paction2 = paction_;
    676           if (ifree) {
     666               printf("Starting expensive\n");
     667#endif
     668
     669               if (dual) {
     670                    int itry;
     671                    for (itry = 0; itry < 5; itry++) {
     672                         paction_ = remove_dual_action::presolve(prob, paction_);
     673                         if (prob->status_)
     674                              break;
     675                         const CoinPresolveAction * const paction2 = paction_;
     676                         if (ifree) {
    677677#ifdef IMPLIED
    678678#if IMPLIED2 ==0
    679             int fill_level=0; // switches off substitution
     679                              int fill_level = 0; // switches off substitution
    680680#elif IMPLIED2!=99
    681             int fill_level=IMPLIED2;
    682 #endif
    683 #endif
    684             if ((itry&1)==0)
    685               paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    686             if (prob->status_)
    687               break;
    688           }
    689           if (paction_ == paction2)
    690             break;
    691         }
    692       } else if (ifree) {
    693         // just free
     681                              int fill_level = IMPLIED2;
     682#endif
     683#endif
     684                              if ((itry & 1) == 0)
     685                                   paction_ = implied_free_action::presolve(prob, paction_, fill_level);
     686                              if (prob->status_)
     687                                   break;
     688                         }
     689                         if (paction_ == paction2)
     690                              break;
     691                    }
     692               } else if (ifree) {
     693                    // just free
    694694#ifdef IMPLIED
    695695#if IMPLIED2 ==0
    696         int fill_level=0; // switches off substitution
     696                    int fill_level = 0; // switches off substitution
    697697#elif IMPLIED2!=99
    698         int fill_level=IMPLIED2;
    699 #endif
    700 #endif
    701         paction_ = implied_free_action::presolve(prob, paction_,fill_level);
    702         if (prob->status_)
    703           break;
    704       }
     698                    int fill_level = IMPLIED2;
     699#endif
     700#endif
     701                    paction_ = implied_free_action::presolve(prob, paction_, fill_level);
     702                    if (prob->status_)
     703                         break;
     704               }
    705705#if     PRESOLVE_DEBUG
    706       check_sol(prob,1.0e0);
    707 #endif
    708       if (dupcol) {
    709         // maybe allow integer columns to be checked
    710         if ((presolveActions_&512)!=0)
    711           prob->setPresolveOptions(prob->presolveOptions()|1);
    712         paction_ = dupcol_action::presolve(prob, paction_);
    713         if (prob->status_)
    714           break;
    715       }
     706               check_sol(prob, 1.0e0);
     707#endif
     708               if (dupcol) {
     709                    // maybe allow integer columns to be checked
     710                    if ((presolveActions_ & 512) != 0)
     711                         prob->setPresolveOptions(prob->presolveOptions() | 1);
     712                    paction_ = dupcol_action::presolve(prob, paction_);
     713                    if (prob->status_)
     714                         break;
     715               }
    716716#if     PRESOLVE_DEBUG
    717         check_sol(prob,1.0e0);
    718 #endif
    719      
    720       if (duprow) {
    721         paction_ = duprow_action::presolve(prob, paction_);
    722         if (prob->status_)
    723           break;
    724       }
     717               check_sol(prob, 1.0e0);
     718#endif
     719
     720               if (duprow) {
     721                    paction_ = duprow_action::presolve(prob, paction_);
     722                    if (prob->status_)
     723                         break;
     724               }
    725725#if     PRESOLVE_DEBUG
    726       check_sol(prob,1.0e0);
    727 #endif
    728       bool stopLoop=false;
    729       {
    730         int * hinrow = prob->hinrow_;
    731         int numberDropped=0;
    732         for (i=0;i<nrows_;i++)
    733           if (!hinrow[i])
    734             numberDropped++;
    735 
    736         prob->messageHandler()->message(COIN_PRESOLVE_PASS,
    737                                         messages)
    738                                           <<numberDropped<<iLoop+1
    739                                           <<CoinMessageEol;
    740         //printf("%d rows dropped after pass %d\n",numberDropped,
    741         //     iLoop+1);
    742         if (numberDropped==lastDropped)
    743           stopLoop=true;
    744         else
    745           lastDropped = numberDropped;
    746       }
    747       // Do this here as not very loopy
    748       if (slackSingleton) {
    749         // On most passes do not touch costed slacks
    750         if (paction_ != paction0&&!stopLoop) {
    751           paction_ = slack_singleton_action::presolve(prob, paction_,NULL);
    752         } else {
    753           // do costed if Clp (at end as ruins rest of presolve)
    754           paction_ = slack_singleton_action::presolve(prob, paction_,rowObjective_);
    755           stopLoop=true;
    756         }
    757       }
     726               check_sol(prob, 1.0e0);
     727#endif
     728               bool stopLoop = false;
     729               {
     730                    int * hinrow = prob->hinrow_;
     731                    int numberDropped = 0;
     732                    for (i = 0; i < nrows_; i++)
     733                         if (!hinrow[i])
     734                              numberDropped++;
     735
     736                    prob->messageHandler()->message(COIN_PRESOLVE_PASS,
     737                                                    messages)
     738                              << numberDropped << iLoop + 1
     739                              << CoinMessageEol;
     740                    //printf("%d rows dropped after pass %d\n",numberDropped,
     741                    //     iLoop+1);
     742                    if (numberDropped == lastDropped)
     743                         stopLoop = true;
     744                    else
     745                         lastDropped = numberDropped;
     746               }
     747               // Do this here as not very loopy
     748               if (slackSingleton) {
     749                    // On most passes do not touch costed slacks
     750                    if (paction_ != paction0 && !stopLoop) {
     751                         paction_ = slack_singleton_action::presolve(prob, paction_, NULL);
     752                    } else {
     753                         // do costed if Clp (at end as ruins rest of presolve)
     754                         paction_ = slack_singleton_action::presolve(prob, paction_, rowObjective_);
     755                         stopLoop = true;
     756                    }
     757               }
    758758#if     PRESOLVE_DEBUG
    759       check_sol(prob,1.0e0);
    760 #endif
    761       if (paction_ == paction0||stopLoop)
    762         break;
    763          
    764     }
    765   }
    766   if (!prob->status_) {
    767     paction_ = drop_zero_coefficients(prob, paction_);
     759               check_sol(prob, 1.0e0);
     760#endif
     761               if (paction_ == paction0 || stopLoop)
     762                    break;
     763
     764          }
     765     }
     766     if (!prob->status_) {
     767          paction_ = drop_zero_coefficients(prob, paction_);
    768768#if     PRESOLVE_DEBUG
    769         check_sol(prob,1.0e0);
    770 #endif
    771 
    772     paction_ = drop_empty_cols_action::presolve(prob, paction_);
    773     paction_ = drop_empty_rows_action::presolve(prob, paction_);
     769          check_sol(prob, 1.0e0);
     770#endif
     771
     772          paction_ = drop_empty_cols_action::presolve(prob, paction_);
     773          paction_ = drop_empty_rows_action::presolve(prob, paction_);
    774774#if     PRESOLVE_DEBUG
    775         check_sol(prob,1.0e0);
    776 #endif
    777   }
    778  
    779   if (prob->status_) {
    780     if (prob->status_==1)
    781           prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
    782                                              messages)
    783                                                <<prob->feasibilityTolerance_
    784                                                <<CoinMessageEol;
    785     else if (prob->status_==2)
    786           prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
    787                                              messages)
    788                                                <<CoinMessageEol;
    789     else
    790           prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
    791                                              messages)
    792                                                <<CoinMessageEol;
    793     // get rid of data
    794     destroyPresolve();
    795   }
    796   return (paction_);
     775          check_sol(prob, 1.0e0);
     776#endif
     777     }
     778
     779     if (prob->status_) {
     780          if (prob->status_ == 1)
     781               prob->messageHandler()->message(COIN_PRESOLVE_INFEAS,
     782                                               messages)
     783                         << prob->feasibilityTolerance_
     784                         << CoinMessageEol;
     785          else if (prob->status_ == 2)
     786               prob->messageHandler()->message(COIN_PRESOLVE_UNBOUND,
     787                                               messages)
     788                         << CoinMessageEol;
     789          else
     790               prob->messageHandler()->message(COIN_PRESOLVE_INFEASUNBOUND,
     791                                               messages)
     792                         << CoinMessageEol;
     793          // get rid of data
     794          destroyPresolve();
     795     }
     796     return (paction_);
    797797}
    798798
     
    804804void ClpPresolve::postsolve(CoinPostsolveMatrix &prob)
    805805{
    806   {
    807     // Check activities
    808     double *colels      = prob.colels_;
    809     int *hrow           = prob.hrow_;
    810     CoinBigIndex *mcstrt                = prob.mcstrt_;
    811     int *hincol         = prob.hincol_;
    812     int *link           = prob.link_;
    813     int ncols           = prob.ncols_;
    814 
    815     char *cdone = prob.cdone_;
    816 
    817     double * csol = prob.sol_;
    818     int nrows = prob.nrows_;
    819 
    820     int colx;
    821 
    822     double * rsol = prob.acts_;
    823     memset(rsol,0,nrows*sizeof(double));
    824 
    825     for (colx = 0; colx < ncols; ++colx) {
    826       if (cdone[colx]) {
    827         CoinBigIndex k = mcstrt[colx];
    828         int nx = hincol[colx];
    829         double solutionValue = csol[colx];
    830         for (int i=0; i<nx; ++i) {
    831           int row = hrow[k];
    832           double coeff = colels[k];
    833           k = link[k];
    834           rsol[row] += solutionValue*coeff;
    835         }
    836       }
    837     }
    838   }
    839   const CoinPresolveAction *paction = paction_;
    840   //#define PRESOLVE_DEBUG 1
     806     {
     807          // Check activities
     808          double *colels        = prob.colels_;
     809          int *hrow             = prob.hrow_;
     810          CoinBigIndex *mcstrt          = prob.mcstrt_;
     811          int *hincol           = prob.hincol_;
     812          int *link             = prob.link_;
     813          int ncols             = prob.ncols_;
     814
     815          char *cdone   = prob.cdone_;
     816
     817          double * csol = prob.sol_;
     818          int nrows = prob.nrows_;
     819
     820          int colx;
     821
     822          double * rsol = prob.acts_;
     823          memset(rsol, 0, nrows * sizeof(double));
     824
     825          for (colx = 0; colx < ncols; ++colx) {
     826               if (cdone[colx]) {
     827                    CoinBigIndex k = mcstrt[colx];
     828                    int nx = hincol[colx];
     829                    double solutionValue = csol[colx];
     830                    for (int i = 0; i < nx; ++i) {
     831                         int row = hrow[k];
     832                         double coeff = colels[k];
     833                         k = link[k];
     834                         rsol[row] += solutionValue * coeff;
     835                    }
     836               }
     837          }
     838     }
     839     const CoinPresolveAction *paction = paction_;
     840     //#define PRESOLVE_DEBUG 1
    841841#if     PRESOLVE_DEBUG
    842   // Check only works after first one
    843   int checkit=-1;
    844 #endif
    845  
    846   while (paction) {
     842     // Check only works after first one
     843     int checkit = -1;
     844#endif
     845
     846     while (paction) {
    847847#if PRESOLVE_DEBUG
    848     printf("POSTSOLVING %s\n", paction->name());
    849 #endif
    850 
    851     paction->postsolve(&prob);
    852    
     848          printf("POSTSOLVING %s\n", paction->name());
     849#endif
     850
     851          paction->postsolve(&prob);
     852
    853853#if     PRESOLVE_DEBUG
    854   {
    855     int nr=0;
    856     int i;
    857     for (i=0;i<prob.nrows_;i++) {
    858       if ((prob.rowstat_[i]&7)==1) {
    859         nr++;
    860       } else if ((prob.rowstat_[i]&7)==2) {
    861         // at ub
    862         assert (prob.acts_[i]>prob.rup_[i]-1.0e-6);
    863       } else if ((prob.rowstat_[i]&7)==3) {
    864         // at lb
    865         assert (prob.acts_[i]<prob.rlo_[i]+1.0e-6);
    866       }
    867     }
    868     int nc=0;
    869     for (i=0;i<prob.ncols_;i++) {
    870       if ((prob.colstat_[i]&7)==1)
    871         nc++;
    872     }
    873     printf("%d rows (%d basic), %d cols (%d basic)\n",prob.nrows_,nr,prob.ncols_,nc);
    874   }
    875     checkit++;
    876     if (prob.colstat_&&checkit>0) {
    877       presolve_check_nbasic(&prob) ;
    878       presolve_check_sol(&prob,2,2,1) ;
    879     }
    880 #endif
    881     paction = paction->next;
     854          {
     855               int nr = 0;
     856               int i;
     857               for (i = 0; i < prob.nrows_; i++) {
     858                    if ((prob.rowstat_[i] & 7) == 1) {
     859                         nr++;
     860                    } else if ((prob.rowstat_[i] & 7) == 2) {
     861                         // at ub
     862                         assert (prob.acts_[i] > prob.rup_[i] - 1.0e-6);
     863                    } else if ((prob.rowstat_[i] & 7) == 3) {
     864                         // at lb
     865                         assert (prob.acts_[i] < prob.rlo_[i] + 1.0e-6);
     866                    }
     867               }
     868               int nc = 0;
     869               for (i = 0; i < prob.ncols_; i++) {
     870                    if ((prob.colstat_[i] & 7) == 1)
     871                         nc++;
     872               }
     873               printf("%d rows (%d basic), %d cols (%d basic)\n", prob.nrows_, nr, prob.ncols_, nc);
     874          }
     875          checkit++;
     876          if (prob.colstat_ && checkit > 0) {
     877               presolve_check_nbasic(&prob) ;
     878               presolve_check_sol(&prob, 2, 2, 1) ;
     879          }
     880#endif
     881          paction = paction->next;
    882882#if     PRESOLVE_DEBUG
    883883//  check_djs(&prob);
    884     if (checkit>0)
    885       presolve_check_reduced_costs(&prob) ;
    886 #endif
    887   }   
     884          if (checkit > 0)
     885               presolve_check_reduced_costs(&prob) ;
     886#endif
     887     }
    888888#if     PRESOLVE_DEBUG
    889   if (prob.colstat_) {
    890     presolve_check_nbasic(&prob) ;
    891     presolve_check_sol(&prob,2,2,1) ;
    892   }
     889     if (prob.colstat_) {
     890          presolve_check_nbasic(&prob) ;
     891          presolve_check_sol(&prob, 2, 2, 1) ;
     892     }
    893893#endif
    894894#undef PRESOLVE_DEBUG
    895  
     895
    896896#if     0 && PRESOLVE_DEBUG
    897   for (i=0; i<ncols0; i++) {
    898     if (!cdone[i]) {
    899       printf("!cdone[%d]\n", i);
    900       abort();
    901     }
    902   }
    903  
    904   for (i=0; i<nrows0; i++) {
    905     if (!rdone[i]) {
    906       printf("!rdone[%d]\n", i);
    907       abort();
    908     }
    909   }
    910  
    911  
    912   for (i=0; i<ncols0; i++) {
    913     if (sol[i] < -1e10 || sol[i] > 1e10)
    914       printf("!!!%d %g\n", i, sol[i]);
    915    
    916   }
    917  
    918  
    919 #endif
    920  
     897     for (i = 0; i < ncols0; i++) {
     898          if (!cdone[i]) {
     899               printf("!cdone[%d]\n", i);
     900               abort();
     901          }
     902     }
     903
     904     for (i = 0; i < nrows0; i++) {
     905          if (!rdone[i]) {
     906               printf("!rdone[%d]\n", i);
     907               abort();
     908          }
     909     }
     910
     911
     912     for (i = 0; i < ncols0; i++) {
     913          if (sol[i] < -1e10 || sol[i] > 1e10)
     914               printf("!!!%d %g\n", i, sol[i]);
     915
     916     }
     917
     918
     919#endif
     920
    921921#if     0 && PRESOLVE_DEBUG
    922   // debug check:  make sure we ended up with same original matrix
    923   {
    924     int identical = 1;
    925    
    926     for (int i=0; i<ncols0; i++) {
    927       PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
    928       CoinBigIndex kcs0 = &prob->mcstrt0[i];
    929       CoinBigIndex kcs = mcstrt[i];
    930       int n = hincol[i];
    931       for (int k=0; k<n; k++) {
    932         CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs+n, hrow);
    933 
    934         if (k1 == kcs+n) {
    935           printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
    936           abort();
    937         }
    938 
    939         if (colels[k1] != &prob->dels0[kcs0+k])
    940           printf("BAD COLEL[%d %d %d]:  %g\n",
    941                 k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
    942 
    943         if (kcs0+k != k1)
    944           identical=0;
    945       }
    946     }
    947     printf("identical? %d\n", identical);
    948   }
     922     // debug check:  make sure we ended up with same original matrix
     923     {
     924          int identical = 1;
     925
     926          for (int i = 0; i < ncols0; i++) {
     927               PRESOLVEASSERT(hincol[i] == &prob->mcstrt0[i+1] - &prob->mcstrt0[i]);
     928               CoinBigIndex kcs0 = &prob->mcstrt0[i];
     929               CoinBigIndex kcs = mcstrt[i];
     930               int n = hincol[i];
     931               for (int k = 0; k < n; k++) {
     932                    CoinBigIndex k1 = presolve_find_row1(&prob->hrow0[kcs0+k], kcs, kcs + n, hrow);
     933
     934                    if (k1 == kcs + n) {
     935                         printf("ROW %d NOT IN COL %d\n", &prob->hrow0[kcs0+k], i);
     936                         abort();
     937                    }
     938
     939                    if (colels[k1] != &prob->dels0[kcs0+k])
     940                         printf("BAD COLEL[%d %d %d]:  %g\n",
     941                                k1, i, &prob->hrow0[kcs0+k], colels[k1] - &prob->dels0[kcs0+k]);
     942
     943                    if (kcs0 + k != k1)
     944                         identical = 0;
     945               }
     946          }
     947          printf("identical? %d\n", identical);
     948     }
    949949#endif
    950950}
     
    959959static inline double getTolerance(const ClpSimplex  *si, ClpDblParam key)
    960960{
    961   double tol;
    962   if (! si->getDblParam(key, tol)) {
    963     CoinPresolveAction::throwCoinError("getDblParam failed",
    964                                       "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
    965   }
    966   return (tol);
     961     double tol;
     962     if (! si->getDblParam(key, tol)) {
     963          CoinPresolveAction::throwCoinError("getDblParam failed",
     964                                             "CoinPrePostsolveMatrix::CoinPrePostsolveMatrix");
     965     }
     966     return (tol);
    967967}
    968968
     
    979979// but we need to reserve enough space for the original problem.
    980980CoinPrePostsolveMatrix::CoinPrePostsolveMatrix(const ClpSimplex * si,
    981                                              int ncols_in,
    982                                              int nrows_in,
    983                                              CoinBigIndex nelems_in,
    984                                                double bulkRatio)
    985   : ncols_(si->getNumCols()),
    986     nrows_(si->getNumRows()),
    987     nelems_(si->getNumElements()),
    988     ncols0_(ncols_in),
    989     nrows0_(nrows_in),
    990     bulkRatio_(bulkRatio),
    991     mcstrt_(new CoinBigIndex[ncols_in+1]),
    992     hincol_(new int[ncols_in+1]),
    993     cost_(new double[ncols_in]),
    994     clo_(new double[ncols_in]),
    995     cup_(new double[ncols_in]),
    996     rlo_(new double[nrows_in]),
    997     rup_(new double[nrows_in]),
    998     originalColumn_(new int[ncols_in]),
    999     originalRow_(new int[nrows_in]),
    1000     ztolzb_(getTolerance(si, ClpPrimalTolerance)),
    1001     ztoldj_(getTolerance(si, ClpDualTolerance)),
    1002     maxmin_(si->getObjSense()),
    1003     sol_(NULL),
    1004     rowduals_(NULL),
    1005     acts_(NULL),
    1006     rcosts_(NULL),
    1007     colstat_(NULL),
    1008     rowstat_(NULL),
    1009     handler_(NULL),
    1010     defaultHandler_(false),
    1011     messages_()
    1012 
    1013 {
    1014   bulk0_ = static_cast<CoinBigIndex> (bulkRatio_*nelems_in);
    1015   hrow_  = new int   [bulk0_];
    1016   colels_ = new double[bulk0_];
    1017   si->getDblParam(ClpObjOffset,originalOffset_);
    1018   int ncols = si->getNumCols();
    1019   int nrows = si->getNumRows();
    1020 
    1021   setMessageHandler(si->messageHandler()) ;
    1022 
    1023   ClpDisjointCopyN(si->getColLower(), ncols, clo_);
    1024   ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
    1025   //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
    1026   double offset;
    1027   ClpDisjointCopyN(si->objectiveAsObject()->gradient(si,si->getColSolution(),offset,true), ncols, cost_);
    1028   ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
    1029   ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
    1030   int i;
    1031   for (i=0;i<ncols_in;i++)
    1032     originalColumn_[i]=i;
    1033   for (i=0;i<nrows_in;i++)
    1034     originalRow_[i]=i;
    1035   sol_=NULL;
    1036   rowduals_=NULL;
    1037   acts_=NULL;
    1038 
    1039   rcosts_=NULL;
    1040   colstat_=NULL;
    1041   rowstat_=NULL;
     981          int ncols_in,
     982          int nrows_in,
     983          CoinBigIndex nelems_in,
     984          double bulkRatio)
     985     : ncols_(si->getNumCols()),
     986       nrows_(si->getNumRows()),
     987       nelems_(si->getNumElements()),
     988       ncols0_(ncols_in),
     989       nrows0_(nrows_in),
     990       bulkRatio_(bulkRatio),
     991       mcstrt_(new CoinBigIndex[ncols_in+1]),
     992       hincol_(new int[ncols_in+1]),
     993       cost_(new double[ncols_in]),
     994       clo_(new double[ncols_in]),
     995       cup_(new double[ncols_in]),
     996       rlo_(new double[nrows_in]),
     997       rup_(new double[nrows_in]),
     998       originalColumn_(new int[ncols_in]),
     999       originalRow_(new int[nrows_in]),
     1000       ztolzb_(getTolerance(si, ClpPrimalTolerance)),
     1001       ztoldj_(getTolerance(si, ClpDualTolerance)),
     1002       maxmin_(si->getObjSense()),
     1003       sol_(NULL),
     1004       rowduals_(NULL),
     1005       acts_(NULL),
     1006       rcosts_(NULL),
     1007       colstat_(NULL),
     1008       rowstat_(NULL),
     1009       handler_(NULL),
     1010       defaultHandler_(false),
     1011       messages_()
     1012
     1013{
     1014     bulk0_ = static_cast<CoinBigIndex> (bulkRatio_ * nelems_in);
     1015     hrow_  = new int   [bulk0_];
     1016     colels_ = new double[bulk0_];
     1017     si->getDblParam(ClpObjOffset, originalOffset_);
     1018     int ncols = si->getNumCols();
     1019     int nrows = si->getNumRows();
     1020
     1021     setMessageHandler(si->messageHandler()) ;
     1022
     1023     ClpDisjointCopyN(si->getColLower(), ncols, clo_);
     1024     ClpDisjointCopyN(si->getColUpper(), ncols, cup_);
     1025     //ClpDisjointCopyN(si->getObjCoefficients(), ncols, cost_);
     1026     double offset;
     1027     ClpDisjointCopyN(si->objectiveAsObject()->gradient(si, si->getColSolution(), offset, true), ncols, cost_);
     1028     ClpDisjointCopyN(si->getRowLower(), nrows,  rlo_);
     1029     ClpDisjointCopyN(si->getRowUpper(), nrows,  rup_);
     1030     int i;
     1031     for (i = 0; i < ncols_in; i++)
     1032          originalColumn_[i] = i;
     1033     for (i = 0; i < nrows_in; i++)
     1034          originalRow_[i] = i;
     1035     sol_ = NULL;
     1036     rowduals_ = NULL;
     1037     acts_ = NULL;
     1038
     1039     rcosts_ = NULL;
     1040     colstat_ = NULL;
     1041     rowstat_ = NULL;
    10421042}
    10431043
     
    10471047static bool isGapFree(const CoinPackedMatrix& matrix)
    10481048{
    1049   const CoinBigIndex * start = matrix.getVectorStarts();
    1050   const int * length = matrix.getVectorLengths();
    1051   int i = matrix.getSizeVectorLengths() - 1;
    1052   // Quick check
    1053   if (matrix.getNumElements()==start[i]) {
    1054     return true;
    1055   } else {
    1056     for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
    1057       if (start[i+1] - start[i] != length[i])
    1058         break;
    1059     }
    1060     return (! (i >= 0));
    1061   }
     1049     const CoinBigIndex * start = matrix.getVectorStarts();
     1050     const int * length = matrix.getVectorLengths();
     1051     int i = matrix.getSizeVectorLengths() - 1;
     1052     // Quick check
     1053     if (matrix.getNumElements() == start[i]) {
     1054          return true;
     1055     } else {
     1056          for (i = matrix.getSizeVectorLengths() - 1; i >= 0; --i) {
     1057               if (start[i+1] - start[i] != length[i])
     1058                    break;
     1059          }
     1060          return (! (i >= 0));
     1061     }
    10621062}
    10631063#if     PRESOLVE_DEBUG
    10641064static void matrix_bounds_ok(const double *lo, const double *up, int n)
    10651065{
    1066   int i;
    1067   for (i=0; i<n; i++) {
    1068     PRESOLVEASSERT(lo[i] <= up[i]);
    1069     PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
    1070     PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
    1071   }
     1066     int i;
     1067     for (i = 0; i < n; i++) {
     1068          PRESOLVEASSERT(lo[i] <= up[i]);
     1069          PRESOLVEASSERT(lo[i] < PRESOLVE_INF);
     1070          PRESOLVEASSERT(-PRESOLVE_INF < up[i]);
     1071     }
    10721072}
    10731073#endif
    10741074CoinPresolveMatrix::CoinPresolveMatrix(int ncols0_in,
    1075                                        double /*maxmin*/,
    1076                                      // end prepost members
    1077 
    1078                                      ClpSimplex * si,
    1079 
    1080                                      // rowrep
    1081                                      int nrows_in,
    1082                                      CoinBigIndex nelems_in,
     1075                                       double /*maxmin*/,
     1076                                       // end prepost members
     1077
     1078                                       ClpSimplex * si,
     1079
     1080                                       // rowrep
     1081                                       int nrows_in,
     1082                                       CoinBigIndex nelems_in,
    10831083                                       bool doStatus,
    10841084                                       double nonLinearValue,
    10851085                                       double bulkRatio) :
    10861086
    1087   CoinPrePostsolveMatrix(si,
    1088                         ncols0_in, nrows_in, nelems_in,bulkRatio),
    1089   clink_(new presolvehlink[ncols0_in+1]),
    1090   rlink_(new presolvehlink[nrows_in+1]),
    1091 
    1092   dobias_(0.0),
    1093 
    1094 
    1095   // temporary init
    1096   integerType_(new unsigned char[ncols0_in]),
    1097   tuning_(false),
    1098   startTime_(0.0),
    1099   feasibilityTolerance_(0.0),
    1100   status_(-1),
    1101   colsToDo_(new int [ncols0_in]),
    1102   numberColsToDo_(0),
    1103   nextColsToDo_(new int[ncols0_in]),
    1104   numberNextColsToDo_(0),
    1105   rowsToDo_(new int [nrows_in]),
    1106   numberRowsToDo_(0),
    1107   nextRowsToDo_(new int[nrows_in]),
    1108   numberNextRowsToDo_(0),
    1109   presolveOptions_(0)
    1110 {
    1111   const int bufsize = bulk0_;
    1112 
    1113   nrows_ = si->getNumRows() ;
    1114 
    1115   // Set up change bits etc
    1116   rowChanged_ = new unsigned char[nrows_];
    1117   memset(rowChanged_,0,nrows_);
    1118   colChanged_ = new unsigned char[ncols_];
    1119   memset(colChanged_,0,ncols_);
    1120   CoinPackedMatrix * m = si->matrix();
    1121 
    1122   // The coefficient matrix is a big hunk of stuff.
    1123   // Do the copy here to try to avoid running out of memory.
    1124 
    1125   const CoinBigIndex * start = m->getVectorStarts();
    1126   const int * row = m->getIndices();
    1127   const double * element = m->getElements();
    1128   int icol,nel=0;
    1129   mcstrt_[0]=0;
    1130   ClpDisjointCopyN(m->getVectorLengths(),ncols_,  hincol_);
    1131   for (icol=0;icol<ncols_;icol++) {
    1132     CoinBigIndex j;
    1133     for (j=start[icol];j<start[icol]+hincol_[icol];j++) {
    1134       hrow_[nel]=row[j];
    1135       if (fabs(element[j])>ZTOLDP)
    1136         colels_[nel++]=element[j];
    1137     }
    1138     mcstrt_[icol+1]=nel;
    1139     hincol_[icol]=nel-mcstrt_[icol];
    1140   }
    1141 
    1142   // same thing for row rep
    1143   CoinPackedMatrix * mRow = new CoinPackedMatrix();
    1144   mRow->setExtraGap(0.0);
    1145   mRow->setExtraMajor(0.0);
    1146   mRow->reverseOrderedCopyOf(*m);
    1147   //mRow->removeGaps();
    1148   //mRow->setExtraGap(0.0);
    1149 
    1150   // Now get rid of matrix
    1151   si->createEmptyMatrix();
    1152 
    1153   double * el = mRow->getMutableElements();
    1154   int * ind = mRow->getMutableIndices();
    1155   CoinBigIndex * strt = mRow->getMutableVectorStarts();
    1156   int * len = mRow->getMutableVectorLengths();
    1157   // Do carefully to save memory
    1158   rowels_ = new double[bulk0_];
    1159   ClpDisjointCopyN(el,      nelems_, rowels_);
    1160   mRow->nullElementArray();
    1161   delete [] el;
    1162   hcol_ = new int[bulk0_];
    1163   ClpDisjointCopyN(ind,       nelems_, hcol_);
    1164   mRow->nullIndexArray();
    1165   delete [] ind;
    1166   mrstrt_ = new CoinBigIndex[nrows_in+1];
    1167   ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
    1168   mRow->nullStartArray();
    1169   mrstrt_[nrows_] = nelems_;
    1170   delete [] strt;
    1171   hinrow_ = new int[nrows_in+1];
    1172   ClpDisjointCopyN(len, nrows_,  hinrow_);
    1173   if (nelems_>nel) {
    1174     nelems_ = nel;
    1175     // Clean any small elements
    1176     int irow;
    1177     nel=0;
    1178     CoinBigIndex start=0;
    1179     for (irow=0;irow<nrows_;irow++) {
    1180       CoinBigIndex j;
    1181       for (j=start;j<start+hinrow_[irow];j++) {
    1182         hcol_[nel]=hcol_[j];
    1183         if (fabs(rowels_[j])>ZTOLDP)
    1184           rowels_[nel++]=rowels_[j];
    1185       }
    1186       start=mrstrt_[irow+1];
    1187       mrstrt_[irow+1]=nel;
    1188       hinrow_[irow]=nel-mrstrt_[irow];
    1189     }
    1190   }
    1191 
    1192   delete mRow;
    1193   if (si->integerInformation()) {
    1194     CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()),ncols_,integerType_);
    1195   } else {
    1196     ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
    1197   }
     1087     CoinPrePostsolveMatrix(si,
     1088                            ncols0_in, nrows_in, nelems_in, bulkRatio),
     1089     clink_(new presolvehlink[ncols0_in+1]),
     1090     rlink_(new presolvehlink[nrows_in+1]),
     1091
     1092     dobias_(0.0),
     1093
     1094
     1095     // temporary init
     1096     integerType_(new unsigned char[ncols0_in]),
     1097     tuning_(false),
     1098     startTime_(0.0),
     1099     feasibilityTolerance_(0.0),
     1100     status_(-1),
     1101     colsToDo_(new int [ncols0_in]),
     1102     numberColsToDo_(0),
     1103     nextColsToDo_(new int[ncols0_in]),
     1104     numberNextColsToDo_(0),
     1105     rowsToDo_(new int [nrows_in]),
     1106     numberRowsToDo_(0),
     1107     nextRowsToDo_(new int[nrows_in]),
     1108     numberNextRowsToDo_(0),
     1109     presolveOptions_(0)
     1110{
     1111     const int bufsize = bulk0_;
     1112
     1113     nrows_ = si->getNumRows() ;
     1114
     1115     // Set up change bits etc
     1116     rowChanged_ = new unsigned char[nrows_];
     1117     memset(rowChanged_, 0, nrows_);
     1118     colChanged_ = new unsigned char[ncols_];
     1119     memset(colChanged_, 0, ncols_);
     1120     CoinPackedMatrix * m = si->matrix();
     1121
     1122     // The coefficient matrix is a big hunk of stuff.
     1123     // Do the copy here to try to avoid running out of memory.
     1124
     1125     const CoinBigIndex * start = m->getVectorStarts();
     1126     const int * row = m->getIndices();
     1127     const double * element = m->getElements();
     1128     int icol, nel = 0;
     1129     mcstrt_[0] = 0;
     1130     ClpDisjointCopyN(m->getVectorLengths(), ncols_,  hincol_);
     1131     for (icol = 0; icol < ncols_; icol++) {
     1132          CoinBigIndex j;
     1133          for (j = start[icol]; j < start[icol] + hincol_[icol]; j++) {
     1134               hrow_[nel] = row[j];
     1135               if (fabs(element[j]) > ZTOLDP)
     1136                    colels_[nel++] = element[j];
     1137          }
     1138          mcstrt_[icol+1] = nel;
     1139          hincol_[icol] = nel - mcstrt_[icol];
     1140     }
     1141
     1142     // same thing for row rep
     1143     CoinPackedMatrix * mRow = new CoinPackedMatrix();
     1144     mRow->setExtraGap(0.0);
     1145     mRow->setExtraMajor(0.0);
     1146     mRow->reverseOrderedCopyOf(*m);
     1147     //mRow->removeGaps();
     1148     //mRow->setExtraGap(0.0);
     1149
     1150     // Now get rid of matrix
     1151     si->createEmptyMatrix();
     1152
     1153     double * el = mRow->getMutableElements();
     1154     int * ind = mRow->getMutableIndices();
     1155     CoinBigIndex * strt = mRow->getMutableVectorStarts();
     1156     int * len = mRow->getMutableVectorLengths();
     1157     // Do carefully to save memory
     1158     rowels_ = new double[bulk0_];
     1159     ClpDisjointCopyN(el,      nelems_, rowels_);
     1160     mRow->nullElementArray();
     1161     delete [] el;
     1162     hcol_ = new int[bulk0_];
     1163     ClpDisjointCopyN(ind,       nelems_, hcol_);
     1164     mRow->nullIndexArray();
     1165     delete [] ind;
     1166     mrstrt_ = new CoinBigIndex[nrows_in+1];
     1167     ClpDisjointCopyN(strt,  nrows_,  mrstrt_);
     1168     mRow->nullStartArray();
     1169     mrstrt_[nrows_] = nelems_;
     1170     delete [] strt;
     1171     hinrow_ = new int[nrows_in+1];
     1172     ClpDisjointCopyN(len, nrows_,  hinrow_);
     1173     if (nelems_ > nel) {
     1174          nelems_ = nel;
     1175          // Clean any small elements
     1176          int irow;
     1177          nel = 0;
     1178          CoinBigIndex start = 0;
     1179          for (irow = 0; irow < nrows_; irow++) {
     1180               CoinBigIndex j;
     1181               for (j = start; j < start + hinrow_[irow]; j++) {
     1182                    hcol_[nel] = hcol_[j];
     1183                    if (fabs(rowels_[j]) > ZTOLDP)
     1184                         rowels_[nel++] = rowels_[j];
     1185               }
     1186               start = mrstrt_[irow+1];
     1187               mrstrt_[irow+1] = nel;
     1188               hinrow_[irow] = nel - mrstrt_[irow];
     1189          }
     1190     }
     1191
     1192     delete mRow;
     1193     if (si->integerInformation()) {
     1194          CoinMemcpyN(reinterpret_cast<unsigned char *> (si->integerInformation()), ncols_, integerType_);
     1195     } else {
     1196          ClpFillN<unsigned char>(integerType_, ncols_, static_cast<unsigned char> (0));
     1197     }
    11981198
    11991199#ifndef SLIM_CLP
    12001200#ifndef NO_RTTI
    1201     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
     1201     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
    12021202#else
    1203     ClpQuadraticObjective * quadraticObj = NULL;
    1204     if (si->objectiveAsObject()->type()==2)
    1205       quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
    1206 #endif
    1207 #endif
    1208   // Set up prohibited bits if needed
    1209   if (nonLinearValue) {
    1210     anyProhibited_ = true;
    1211     for (icol=0;icol<ncols_;icol++) {
    1212       int j;
    1213       bool nonLinearColumn = false;
    1214       if (cost_[icol]==nonLinearValue)
    1215         nonLinearColumn=true;
    1216       for (j=mcstrt_[icol];j<mcstrt_[icol+1];j++) {
    1217         if (colels_[j]==nonLinearValue) {
    1218           nonLinearColumn=true;
    1219           setRowProhibited(hrow_[j]);
    1220         }
    1221       }
    1222       if (nonLinearColumn)
    1223         setColProhibited(icol);
    1224     }
     1203     ClpQuadraticObjective * quadraticObj = NULL;
     1204     if (si->objectiveAsObject()->type() == 2)
     1205          quadraticObj = (static_cast< ClpQuadraticObjective*>(si->objectiveAsObject()));
     1206#endif
     1207#endif
     1208     // Set up prohibited bits if needed
     1209     if (nonLinearValue) {
     1210          anyProhibited_ = true;
     1211          for (icol = 0; icol < ncols_; icol++) {
     1212               int j;
     1213               bool nonLinearColumn = false;
     1214               if (cost_[icol] == nonLinearValue)
     1215                    nonLinearColumn = true;
     1216               for (j = mcstrt_[icol]; j < mcstrt_[icol+1]; j++) {
     1217                    if (colels_[j] == nonLinearValue) {
     1218                         nonLinearColumn = true;
     1219                         setRowProhibited(hrow_[j]);
     1220                    }
     1221               }
     1222               if (nonLinearColumn)
     1223                    setColProhibited(icol);
     1224          }
    12251225#ifndef SLIM_CLP
    1226   } else if (quadraticObj) {
    1227     CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    1228     //const int * columnQuadratic = quadratic->getIndices();
    1229     //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    1230     const int * columnQuadraticLength = quadratic->getVectorLengths();
    1231     //double * quadraticElement = quadratic->getMutableElements();
    1232     int numberColumns = quadratic->getNumCols();
    1233     anyProhibited_ = true;
    1234     for (int iColumn=0;iColumn<numberColumns;iColumn++) {
    1235       if (columnQuadraticLength[iColumn]) {
    1236         setColProhibited(iColumn);
    1237         //printf("%d prohib\n",iColumn);
    1238       }
    1239     }
    1240 #endif
    1241   } else {
    1242     anyProhibited_ = false;
    1243   }
    1244 
    1245   if (doStatus) {
    1246     // allow for status and solution
    1247     sol_ = new double[ncols_];
    1248     CoinMemcpyN(si->primalColumnSolution(),ncols_,sol_);;
    1249     acts_ = new double [nrows_];
    1250     CoinMemcpyN(si->primalRowSolution(),nrows_,acts_);
    1251     if (!si->statusArray())
    1252       si->createStatus();
    1253     colstat_ = new unsigned char [nrows_+ncols_];
    1254     CoinMemcpyN(si->statusArray(),      (nrows_+ncols_),colstat_);
    1255     rowstat_ = colstat_+ncols_;
    1256   }
    1257 
    1258   // the original model's fields are now unneeded - free them
    1259  
    1260   si->resize(0,0);
     1226     } else if (quadraticObj) {
     1227          CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1228          //const int * columnQuadratic = quadratic->getIndices();
     1229          //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1230          const int * columnQuadraticLength = quadratic->getVectorLengths();
     1231          //double * quadraticElement = quadratic->getMutableElements();
     1232          int numberColumns = quadratic->getNumCols();
     1233          anyProhibited_ = true;
     1234          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
     1235               if (columnQuadraticLength[iColumn]) {
     1236                    setColProhibited(iColumn);
     1237                    //printf("%d prohib\n",iColumn);
     1238               }
     1239          }
     1240#endif
     1241     } else {
     1242          anyProhibited_ = false;
     1243     }
     1244
     1245     if (doStatus) {
     1246          // allow for status and solution
     1247          sol_ = new double[ncols_];
     1248          CoinMemcpyN(si->primalColumnSolution(), ncols_, sol_);;
     1249          acts_ = new double [nrows_];
     1250          CoinMemcpyN(si->primalRowSolution(), nrows_, acts_);
     1251          if (!si->statusArray())
     1252               si->createStatus();
     1253          colstat_ = new unsigned char [nrows_+ncols_];
     1254          CoinMemcpyN(si->statusArray(),        (nrows_ + ncols_), colstat_);
     1255          rowstat_ = colstat_ + ncols_;
     1256     }
     1257
     1258     // the original model's fields are now unneeded - free them
     1259
     1260     si->resize(0, 0);
    12611261
    12621262#if     PRESOLVE_DEBUG
    1263   matrix_bounds_ok(rlo_, rup_, nrows_);
    1264   matrix_bounds_ok(clo_, cup_, ncols_);
     1263     matrix_bounds_ok(rlo_, rup_, nrows_);
     1264     matrix_bounds_ok(clo_, cup_, ncols_);
    12651265#endif
    12661266
    12671267#if 0
    1268   for (i=0; i<nrows; ++i)
    1269     printf("NR: %6d\n", hinrow[i]);
    1270   for (int i=0; i<ncols; ++i)
    1271     printf("NC: %6d\n", hincol[i]);
    1272 #endif
    1273 
    1274   presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
    1275   presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
    1276 
    1277   // this allows last col/row to expand up to bufsize-1 (22);
    1278   // this must come after the calls to presolve_prefix
    1279   mcstrt_[ncols_] = bufsize-1;
    1280   mrstrt_[nrows_] = bufsize-1;
    1281   // Allocate useful arrays
    1282   initializeStuff();
     1268     for (i = 0; i < nrows; ++i)
     1269          printf("NR: %6d\n", hinrow[i]);
     1270     for (int i = 0; i < ncols; ++i)
     1271          printf("NC: %6d\n", hincol[i]);
     1272#endif
     1273
     1274     presolve_make_memlists(/*mcstrt_,*/ hincol_, clink_, ncols_);
     1275     presolve_make_memlists(/*mrstrt_,*/ hinrow_, rlink_, nrows_);
     1276
     1277     // this allows last col/row to expand up to bufsize-1 (22);
     1278     // this must come after the calls to presolve_prefix
     1279     mcstrt_[ncols_] = bufsize - 1;
     1280     mrstrt_[nrows_] = bufsize - 1;
     1281     // Allocate useful arrays
     1282     initializeStuff();
    12831283
    12841284#if     PRESOLVE_CONSISTENCY
    12851285//consistent(false);
    1286   presolve_consistent(this,false) ;
     1286     presolve_consistent(this, false) ;
    12871287#endif
    12881288}
    12891289
    12901290void CoinPresolveMatrix::update_model(ClpSimplex * si,
    1291                                       int /*nrows0*/,
    1292                                       int /*ncols0*/,
    1293                                       CoinBigIndex /*nelems0*/)
    1294 {
    1295   si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
    1296                 clo_, cup_, cost_, rlo_, rup_);
    1297   //delete [] si->integerInformation();
    1298   int numberIntegers=0;
    1299   for (int i=0; i<ncols_; i++) {
    1300     if (integerType_[i])
    1301       numberIntegers++;
    1302   }
    1303   if (numberIntegers)
    1304     si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
    1305   else
    1306     si->copyInIntegerInformation(NULL);
     1291                                      int /*nrows0*/,
     1292                                      int /*ncols0*/,
     1293                                      CoinBigIndex /*nelems0*/)
     1294{
     1295     si->loadProblem(ncols_, nrows_, mcstrt_, hrow_, colels_, hincol_,
     1296                    clo_, cup_, cost_, rlo_, rup_);
     1297     //delete [] si->integerInformation();
     1298     int numberIntegers = 0;
     1299     for (int i = 0; i < ncols_; i++) {
     1300          if (integerType_[i])
     1301               numberIntegers++;
     1302     }
     1303     if (numberIntegers)
     1304          si->copyInIntegerInformation(reinterpret_cast<const char *> (integerType_));
     1305     else
     1306          si->copyInIntegerInformation(NULL);
    13071307
    13081308#if     PRESOLVE_SUMMARY
    1309   printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
    1310          ncols_, ncols0-ncols_,
    1311          nrows_, nrows0-nrows_,
    1312          si->getNumElements(), nelems0-si->getNumElements());
    1313 #endif
    1314   si->setDblParam(ClpObjOffset,originalOffset_-dobias_);
     1309     printf("NEW NCOL/NROW/NELS:  %d(-%d) %d(-%d) %d(-%d)\n",
     1310            ncols_, ncols0 - ncols_,
     1311            nrows_, nrows0 - nrows_,
     1312            si->getNumElements(), nelems0 - si->getNumElements());
     1313#endif
     1314     si->setDblParam(ClpObjOffset, originalOffset_ - dobias_);
    13151315
    13161316}
     
    13291329
    13301330CoinPostsolveMatrix::CoinPostsolveMatrix(ClpSimplex*  si,
    1331                                        int ncols0_in,
    1332                                        int nrows0_in,
    1333                                        CoinBigIndex nelems0,
    1334                                    
    1335                                        double maxmin,
    1336                                        // end prepost members
    1337 
    1338                                        double *sol_in,
    1339                                        double *acts_in,
    1340 
    1341                                        unsigned char *colstat_in,
    1342                                        unsigned char *rowstat_in) :
    1343   CoinPrePostsolveMatrix(si,
    1344                         ncols0_in, nrows0_in, nelems0,2.0),
    1345 
    1346   free_list_(0),
    1347   // link, free_list, maxlink
    1348   maxlink_(bulk0_),
    1349   link_(new int[/*maxlink*/ bulk0_]),
    1350      
    1351   cdone_(new char[ncols0_]),
    1352   rdone_(new char[nrows0_in])
    1353 
    1354 {
    1355   bulk0_ = maxlink_ ;
    1356   nrows_ = si->getNumRows() ;
    1357   ncols_ = si->getNumCols() ;
    1358 
    1359   sol_=sol_in;
    1360   rowduals_=NULL;
    1361   acts_=acts_in;
    1362 
    1363   rcosts_=NULL;
    1364   colstat_=colstat_in;
    1365   rowstat_=rowstat_in;
    1366 
    1367   // this is the *reduced* model, which is probably smaller
    1368   int ncols1 = ncols_ ;
    1369   int nrows1 = nrows_ ;
    1370 
    1371   const CoinPackedMatrix * m = si->matrix();
    1372 
    1373   const CoinBigIndex nelemsr = m->getNumElements();
    1374   if (m->getNumElements()&&!isGapFree(*m)) {
    1375     // Odd - gaps
    1376     CoinPackedMatrix mm(*m);
    1377     mm.removeGaps();
    1378     mm.setExtraGap(0.0);
    1379    
    1380     ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
    1381     CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    1382     mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
    1383     ClpDisjointCopyN(mm.getVectorLengths(),ncols1,  hincol_);
    1384     ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
    1385     ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
    1386   } else {
    1387     // No gaps
    1388    
    1389     ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
    1390     CoinZeroN(mcstrt_+ncols1,ncols0_-ncols1);
    1391     mcstrt_[ncols1] = nelems0;  // ??    (should point to end of bulk store   -- lh --)
    1392     ClpDisjointCopyN(m->getVectorLengths(),ncols1,  hincol_);
    1393     ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
    1394     ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
    1395   }
     1331          int ncols0_in,
     1332          int nrows0_in,
     1333          CoinBigIndex nelems0,
     1334
     1335          double maxmin,
     1336          // end prepost members
     1337
     1338          double *sol_in,
     1339          double *acts_in,
     1340
     1341          unsigned char *colstat_in,
     1342          unsigned char *rowstat_in) :
     1343     CoinPrePostsolveMatrix(si,
     1344                            ncols0_in, nrows0_in, nelems0, 2.0),
     1345
     1346     free_list_(0),
     1347     // link, free_list, maxlink
     1348     maxlink_(bulk0_),
     1349     link_(new int[/*maxlink*/ bulk0_]),
     1350
     1351     cdone_(new char[ncols0_]),
     1352     rdone_(new char[nrows0_in])
     1353
     1354{
     1355     bulk0_ = maxlink_ ;
     1356     nrows_ = si->getNumRows() ;
     1357     ncols_ = si->getNumCols() ;
     1358
     1359     sol_ = sol_in;
     1360     rowduals_ = NULL;
     1361     acts_ = acts_in;
     1362
     1363     rcosts_ = NULL;
     1364     colstat_ = colstat_in;
     1365     rowstat_ = rowstat_in;
     1366
     1367     // this is the *reduced* model, which is probably smaller
     1368     int ncols1 = ncols_ ;
     1369     int nrows1 = nrows_ ;
     1370
     1371     const CoinPackedMatrix * m = si->matrix();
     1372
     1373     const CoinBigIndex nelemsr = m->getNumElements();
     1374     if (m->getNumElements() && !isGapFree(*m)) {
     1375          // Odd - gaps
     1376          CoinPackedMatrix mm(*m);
     1377          mm.removeGaps();
     1378          mm.setExtraGap(0.0);
     1379
     1380          ClpDisjointCopyN(mm.getVectorStarts(), ncols1, mcstrt_);
     1381          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
     1382          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
     1383          ClpDisjointCopyN(mm.getVectorLengths(), ncols1,  hincol_);
     1384          ClpDisjointCopyN(mm.getIndices(),      nelemsr, hrow_);
     1385          ClpDisjointCopyN(mm.getElements(),     nelemsr, colels_);
     1386     } else {
     1387          // No gaps
     1388
     1389          ClpDisjointCopyN(m->getVectorStarts(), ncols1, mcstrt_);
     1390          CoinZeroN(mcstrt_ + ncols1, ncols0_ - ncols1);
     1391          mcstrt_[ncols1] = nelems0;    // ??    (should point to end of bulk store   -- lh --)
     1392          ClpDisjointCopyN(m->getVectorLengths(), ncols1,  hincol_);
     1393          ClpDisjointCopyN(m->getIndices(),      nelemsr, hrow_);
     1394          ClpDisjointCopyN(m->getElements(),     nelemsr, colels_);
     1395     }
    13961396
    13971397
    13981398
    13991399#if     0 && PRESOLVE_DEBUG
    1400   presolve_check_costs(model, &colcopy);
    1401 #endif
    1402 
    1403   // This determines the size of the data structure that contains
    1404   // the matrix being postsolved.  Links are taken from the free_list
    1405   // to recreate matrix entries that were presolved away,
    1406   // and links are added to the free_list when entries created during
    1407   // presolve are discarded.  There is never a need to gc this list.
    1408   // Naturally, it should contain
    1409   // exactly nelems0 entries "in use" when postsolving is done,
    1410   // but I don't know whether the matrix could temporarily get
    1411   // larger during postsolving.  Substitution into more than two
    1412   // rows could do that, in principle.  I am being very conservative
    1413   // here by reserving much more than the amount of space I probably need.
    1414   // If this guess is wrong, check_free_list may be called.
    1415   //  int bufsize = 2*nelems0;
    1416 
    1417   memset(cdone_, -1, ncols0_);
    1418   memset(rdone_, -1, nrows0_);
    1419 
    1420   rowduals_ = new double[nrows0_];
    1421   ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
    1422 
    1423   rcosts_ = new double[ncols0_];
    1424   ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
    1425   if (maxmin<0.0) {
    1426     // change so will look as if minimize
    1427     int i;
    1428     for (i=0;i<nrows1;i++)
    1429       rowduals_[i] = - rowduals_[i];
    1430     for (i=0;i<ncols1;i++) {
    1431       rcosts_[i] = - rcosts_[i];
    1432     }
    1433   }
    1434 
    1435   //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
    1436   //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
    1437 
    1438   ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
    1439   si->setDblParam(ClpObjOffset,originalOffset_);
    1440 
    1441   for (int j=0; j<ncols1; j++) {
    1442     CoinBigIndex kcs = mcstrt_[j];
    1443     CoinBigIndex kce = kcs + hincol_[j];
    1444     for (CoinBigIndex k=kcs; k<kce; ++k) {
    1445       link_[k] = k+1;
    1446     }
    1447     link_[kce-1] = NO_LINK ;
    1448   }
    1449   {
    1450     int ml = maxlink_;
    1451     for (CoinBigIndex k=nelemsr; k<ml; ++k)
    1452       link_[k] = k+1;
    1453     if (ml)
    1454       link_[ml-1] = NO_LINK;
    1455   }
    1456   free_list_ = nelemsr;
     1400     presolve_check_costs(model, &colcopy);
     1401#endif
     1402
     1403     // This determines the size of the data structure that contains
     1404     // the matrix being postsolved.  Links are taken from the free_list
     1405     // to recreate matrix entries that were presolved away,
     1406     // and links are added to the free_list when entries created during
     1407     // presolve are discarded.  There is never a need to gc this list.
     1408     // Naturally, it should contain
     1409     // exactly nelems0 entries "in use" when postsolving is done,
     1410     // but I don't know whether the matrix could temporarily get
     1411     // larger during postsolving.  Substitution into more than two
     1412     // rows could do that, in principle.  I am being very conservative
     1413     // here by reserving much more than the amount of space I probably need.
     1414     // If this guess is wrong, check_free_list may be called.
     1415     //  int bufsize = 2*nelems0;
     1416
     1417     memset(cdone_, -1, ncols0_);
     1418     memset(rdone_, -1, nrows0_);
     1419
     1420     rowduals_ = new double[nrows0_];
     1421     ClpDisjointCopyN(si->getRowPrice(), nrows1, rowduals_);
     1422
     1423     rcosts_ = new double[ncols0_];
     1424     ClpDisjointCopyN(si->getReducedCost(), ncols1, rcosts_);
     1425     if (maxmin < 0.0) {
     1426          // change so will look as if minimize
     1427          int i;
     1428          for (i = 0; i < nrows1; i++)
     1429               rowduals_[i] = - rowduals_[i];
     1430          for (i = 0; i < ncols1; i++) {
     1431               rcosts_[i] = - rcosts_[i];
     1432          }
     1433     }
     1434
     1435     //ClpDisjointCopyN(si->getRowUpper(), nrows1, rup_);
     1436     //ClpDisjointCopyN(si->getRowLower(), nrows1, rlo_);
     1437
     1438     ClpDisjointCopyN(si->getColSolution(), ncols1, sol_);
     1439     si->setDblParam(ClpObjOffset, originalOffset_);
     1440
     1441     for (int j = 0; j < ncols1; j++) {
     1442          CoinBigIndex kcs = mcstrt_[j];
     1443          CoinBigIndex kce = kcs + hincol_[j];
     1444          for (CoinBigIndex k = kcs; k < kce; ++k) {
     1445               link_[k] = k + 1;
     1446          }
     1447          link_[kce-1] = NO_LINK ;
     1448     }
     1449     {
     1450          int ml = maxlink_;
     1451          for (CoinBigIndex k = nelemsr; k < ml; ++k)
     1452               link_[k] = k + 1;
     1453          if (ml)
     1454               link_[ml-1] = NO_LINK;
     1455     }
     1456     free_list_ = nelemsr;
    14571457# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
    1458 /*
    1459   These are used to track the action of postsolve transforms during debugging.
    1460 */
    1461   CoinFillN(cdone_,ncols1,PRESENT_IN_REDUCED) ;
    1462   CoinZeroN(cdone_+ncols1,ncols0_in-ncols1) ;
    1463   CoinFillN(rdone_,nrows1,PRESENT_IN_REDUCED) ;
    1464   CoinZeroN(rdone_+nrows1,nrows0_in-nrows1) ;
     1458     /*
     1459       These are used to track the action of postsolve transforms during debugging.
     1460     */
     1461     CoinFillN(cdone_, ncols1, PRESENT_IN_REDUCED) ;
     1462     CoinZeroN(cdone_ + ncols1, ncols0_in - ncols1) ;
     1463     CoinFillN(rdone_, nrows1, PRESENT_IN_REDUCED) ;
     1464     CoinZeroN(rdone_ + nrows1, nrows0_in - nrows1) ;
    14651465# endif
    14661466}
    14671467/* This is main part of Presolve */
    1468 ClpSimplex * 
     1468ClpSimplex *
    14691469ClpPresolve::gutsOfPresolvedModel(ClpSimplex * originalModel,
    1470                                   double feasibilityTolerance,
    1471                                   bool keepIntegers,
    1472                                   int numberPasses,
    1473                                   bool dropNames,
     1470                                  double feasibilityTolerance,
     1471                                  bool keepIntegers,
     1472                                  int numberPasses,
     1473                                  bool dropNames,
    14741474                                  bool doRowObjective)
    14751475{
    1476   ncols_ = originalModel->getNumCols();
    1477   nrows_ = originalModel->getNumRows();
    1478   nelems_ = originalModel->getNumElements();
    1479   numberPasses_ = numberPasses;
    1480 
    1481   double maxmin = originalModel->getObjSense();
    1482   originalModel_ = originalModel;
    1483   delete [] originalColumn_;
    1484   originalColumn_ = new int[ncols_];
    1485   delete [] originalRow_;
    1486   originalRow_ = new int[nrows_];
    1487   // and fill in case returns early
    1488   int i;
    1489   for (i=0;i<ncols_;i++)
    1490     originalColumn_[i]=i;
    1491   for (i=0;i<nrows_;i++)
    1492     originalRow_[i]=i;
    1493   delete [] rowObjective_;
    1494   if (doRowObjective) {
    1495     rowObjective_ = new double [nrows_];
    1496     memset(rowObjective_,0,nrows_*sizeof(double));
    1497   } else {
    1498     rowObjective_=NULL;
    1499   }
    1500 
    1501   // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
    1502   int result = -1;
    1503  
    1504   // User may have deleted - its their responsibility
    1505   presolvedModel_=NULL;
    1506   // Messages
    1507   CoinMessages messages = originalModel->coinMessages();
    1508   // Only go round 100 times even if integer preprocessing
    1509   int totalPasses=100;
    1510   while (result==-1) {
     1476     ncols_ = originalModel->getNumCols();
     1477     nrows_ = originalModel->getNumRows();
     1478     nelems_ = originalModel->getNumElements();
     1479     numberPasses_ = numberPasses;
     1480
     1481     double maxmin = originalModel->getObjSense();
     1482     originalModel_ = originalModel;
     1483     delete [] originalColumn_;
     1484     originalColumn_ = new int[ncols_];
     1485     delete [] originalRow_;
     1486     originalRow_ = new int[nrows_];
     1487     // and fill in case returns early
     1488     int i;
     1489     for (i = 0; i < ncols_; i++)
     1490          originalColumn_[i] = i;
     1491     for (i = 0; i < nrows_; i++)
     1492          originalRow_[i] = i;
     1493     delete [] rowObjective_;
     1494     if (doRowObjective) {
     1495          rowObjective_ = new double [nrows_];
     1496          memset(rowObjective_, 0, nrows_ * sizeof(double));
     1497     } else {
     1498          rowObjective_ = NULL;
     1499     }
     1500
     1501     // result is 0 - okay, 1 infeasible, -1 go round again, 2 - original model
     1502     int result = -1;
     1503
     1504     // User may have deleted - its their responsibility
     1505     presolvedModel_ = NULL;
     1506     // Messages
     1507     CoinMessages messages = originalModel->coinMessages();
     1508     // Only go round 100 times even if integer preprocessing
     1509     int totalPasses = 100;
     1510     while (result == -1) {
    15111511
    15121512#ifndef CLP_NO_STD
    1513     // make new copy
    1514     if (saveFile_=="") {
    1515 #endif
    1516       delete presolvedModel_;
     1513          // make new copy
     1514          if (saveFile_ == "") {
     1515#endif
     1516               delete presolvedModel_;
    15171517#ifndef CLP_NO_STD
    1518       // So won't get names
    1519       int lengthNames = originalModel->lengthNames();
    1520       originalModel->setLengthNames(0);
    1521 #endif
    1522       presolvedModel_ = new ClpSimplex(*originalModel);
     1518               // So won't get names
     1519               int lengthNames = originalModel->lengthNames();
     1520               originalModel->setLengthNames(0);
     1521#endif
     1522               presolvedModel_ = new ClpSimplex(*originalModel);
    15231523#ifndef CLP_NO_STD
    1524       originalModel->setLengthNames(lengthNames);
    1525     } else {
    1526       presolvedModel_=originalModel;
    1527     }
    1528     presolvedModel_->dropNames();
    1529 #endif
    1530 
    1531     // drop integer information if wanted
    1532     if (!keepIntegers)
    1533       presolvedModel_->deleteIntegerInformation();
    1534     totalPasses--;
    1535 
    1536     double ratio=2.0;
    1537     if (substitution_>3)
    1538       ratio=substitution_;
    1539     else if (substitution_==2)
    1540       ratio=1.5;
    1541     CoinPresolveMatrix prob(ncols_,
    1542                         maxmin,
    1543                         presolvedModel_,
    1544                         nrows_, nelems_,true,nonLinearValue_,ratio);
    1545     prob.setMaximumSubstitutionLevel(substitution_);
    1546     if (doRowObjective)
    1547       memset(rowObjective_,0,nrows_*sizeof(double));
    1548     // See if we want statistics
    1549     if ((presolveActions_&0x80000000)!=0)
    1550       prob.statistics();
    1551     // make sure row solution correct
    1552     {
    1553       double *colels    = prob.colels_;
    1554       int *hrow         = prob.hrow_;
    1555       CoinBigIndex *mcstrt              = prob.mcstrt_;
    1556       int *hincol               = prob.hincol_;
    1557       int ncols         = prob.ncols_;
    1558      
    1559      
    1560       double * csol = prob.sol_;
    1561       double * acts = prob.acts_;
    1562       int nrows = prob.nrows_;
    1563 
    1564       int colx;
    1565 
    1566       memset(acts,0,nrows*sizeof(double));
    1567      
    1568       for (colx = 0; colx < ncols; ++colx) {
    1569         double solutionValue = csol[colx];
    1570         for (int i=mcstrt[colx]; i<mcstrt[colx]+hincol[colx]; ++i) {
    1571           int row = hrow[i];
    1572           double coeff = colels[i];
    1573           acts[row] += solutionValue*coeff;
    1574         }
    1575       }
    1576     }
    1577 
    1578     // move across feasibility tolerance
    1579     prob.feasibilityTolerance_ = feasibilityTolerance;
    1580 
    1581     // Do presolve
    1582     paction_ = presolve(&prob);
    1583     // Get rid of useful arrays
    1584     prob.deleteStuff();
    1585 
    1586     result =0;
    1587 
    1588     if (prob.status_==0&&paction_) {
    1589       // Looks feasible but double check to see if anything slipped through
    1590       int n             = prob.ncols_;
    1591       double * lo = prob.clo_;
    1592       double * up = prob.cup_;
    1593       int i;
    1594      
    1595       for (i=0;i<n;i++) {
    1596         if (up[i]<lo[i]) {
    1597           if (up[i]<lo[i]-1.0e-8) {
    1598             // infeasible
    1599             prob.status_=1;
    1600           } else {
    1601             up[i]=lo[i];
    1602           }
    1603         }
    1604       }
    1605      
    1606       n = prob.nrows_;
    1607       lo = prob.rlo_;
    1608       up = prob.rup_;
    1609 
    1610       for (i=0;i<n;i++) {
    1611         if (up[i]<lo[i]) {
    1612           if (up[i]<lo[i]-1.0e-8) {
    1613             // infeasible
    1614             prob.status_=1;
    1615           } else {
    1616             up[i]=lo[i];
    1617           }
    1618         }
    1619       }
    1620     }
    1621     if (prob.status_==0&&paction_) {
    1622       // feasible
    1623    
    1624       prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
    1625       // copy status and solution
    1626       CoinMemcpyN(           prob.sol_,prob.ncols_,presolvedModel_->primalColumnSolution());
    1627       CoinMemcpyN(           prob.acts_,prob.nrows_,presolvedModel_->primalRowSolution());
    1628       CoinMemcpyN(           prob.colstat_,prob.ncols_,presolvedModel_->statusArray());
    1629       CoinMemcpyN(           prob.rowstat_,prob.nrows_,presolvedModel_->statusArray()+prob.ncols_);
    1630       delete [] prob.sol_;
    1631       delete [] prob.acts_;
    1632       delete [] prob.colstat_;
    1633       prob.sol_=NULL;
    1634       prob.acts_=NULL;
    1635       prob.colstat_=NULL;
    1636      
    1637       int ncolsNow = presolvedModel_->getNumCols();
    1638       CoinMemcpyN(prob.originalColumn_,ncolsNow,originalColumn_);
     1524               originalModel->setLengthNames(lengthNames);
     1525          } else {
     1526               presolvedModel_ = originalModel;
     1527          }
     1528          presolvedModel_->dropNames();
     1529#endif
     1530
     1531          // drop integer information if wanted
     1532          if (!keepIntegers)
     1533               presolvedModel_->deleteIntegerInformation();
     1534          totalPasses--;
     1535
     1536          double ratio = 2.0;
     1537          if (substitution_ > 3)
     1538               ratio = substitution_;
     1539          else if (substitution_ == 2)
     1540               ratio = 1.5;
     1541          CoinPresolveMatrix prob(ncols_,
     1542                                  maxmin,
     1543                                  presolvedModel_,
     1544                                  nrows_, nelems_, true, nonLinearValue_, ratio);
     1545          prob.setMaximumSubstitutionLevel(substitution_);
     1546          if (doRowObjective)
     1547               memset(rowObjective_, 0, nrows_ * sizeof(double));
     1548          // See if we want statistics
     1549          if ((presolveActions_ & 0x80000000) != 0)
     1550               prob.statistics();
     1551          // make sure row solution correct
     1552          {
     1553               double *colels   = prob.colels_;
     1554               int *hrow                = prob.hrow_;
     1555               CoinBigIndex *mcstrt             = prob.mcstrt_;
     1556               int *hincol              = prob.hincol_;
     1557               int ncols                = prob.ncols_;
     1558
     1559
     1560               double * csol = prob.sol_;
     1561               double * acts = prob.acts_;
     1562               int nrows = prob.nrows_;
     1563
     1564               int colx;
     1565
     1566               memset(acts, 0, nrows * sizeof(double));
     1567
     1568               for (colx = 0; colx < ncols; ++colx) {
     1569                    double solutionValue = csol[colx];
     1570                    for (int i = mcstrt[colx]; i < mcstrt[colx] + hincol[colx]; ++i) {
     1571                         int row = hrow[i];
     1572                         double coeff = colels[i];
     1573                         acts[row] += solutionValue * coeff;
     1574                    }
     1575               }
     1576          }
     1577
     1578          // move across feasibility tolerance
     1579          prob.feasibilityTolerance_ = feasibilityTolerance;
     1580
     1581          // Do presolve
     1582          paction_ = presolve(&prob);
     1583          // Get rid of useful arrays
     1584          prob.deleteStuff();
     1585
     1586          result = 0;
     1587
     1588          if (prob.status_ == 0 && paction_) {
     1589               // Looks feasible but double check to see if anything slipped through
     1590               int n            = prob.ncols_;
     1591               double * lo = prob.clo_;
     1592               double * up = prob.cup_;
     1593               int i;
     1594
     1595               for (i = 0; i < n; i++) {
     1596                    if (up[i] < lo[i]) {
     1597                         if (up[i] < lo[i] - 1.0e-8) {
     1598                              // infeasible
     1599                              prob.status_ = 1;
     1600                         } else {
     1601                              up[i] = lo[i];
     1602                         }
     1603                    }
     1604               }
     1605
     1606               n = prob.nrows_;
     1607               lo = prob.rlo_;
     1608               up = prob.rup_;
     1609
     1610               for (i = 0; i < n; i++) {
     1611                    if (up[i] < lo[i]) {
     1612                         if (up[i] < lo[i] - 1.0e-8) {
     1613                              // infeasible
     1614                              prob.status_ = 1;
     1615                         } else {
     1616                              up[i] = lo[i];
     1617                         }
     1618                    }
     1619               }
     1620          }
     1621          if (prob.status_ == 0 && paction_) {
     1622               // feasible
     1623
     1624               prob.update_model(presolvedModel_, nrows_, ncols_, nelems_);
     1625               // copy status and solution
     1626               CoinMemcpyN(          prob.sol_, prob.ncols_, presolvedModel_->primalColumnSolution());
     1627               CoinMemcpyN(          prob.acts_, prob.nrows_, presolvedModel_->primalRowSolution());
     1628               CoinMemcpyN(          prob.colstat_, prob.ncols_, presolvedModel_->statusArray());
     1629               CoinMemcpyN(          prob.rowstat_, prob.nrows_, presolvedModel_->statusArray() + prob.ncols_);
     1630               delete [] prob.sol_;
     1631               delete [] prob.acts_;
     1632               delete [] prob.colstat_;
     1633               prob.sol_ = NULL;
     1634               prob.acts_ = NULL;
     1635               prob.colstat_ = NULL;
     1636
     1637               int ncolsNow = presolvedModel_->getNumCols();
     1638               CoinMemcpyN(prob.originalColumn_, ncolsNow, originalColumn_);
    16391639#ifndef SLIM_CLP
    16401640#ifndef NO_RTTI
    1641       ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
     1641               ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
    16421642#else
    1643       ClpQuadraticObjective * quadraticObj = NULL;
    1644       if (originalModel->objectiveAsObject()->type()==2)
    1645         quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
    1646 #endif
    1647       if (quadraticObj) {
    1648         // set up for subset
    1649         char * mark = new char [ncols_];
    1650         memset(mark,0,ncols_);
    1651         CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
    1652         //const int * columnQuadratic = quadratic->getIndices();
    1653         //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
    1654         const int * columnQuadraticLength = quadratic->getVectorLengths();
    1655         //double * quadraticElement = quadratic->getMutableElements();
    1656         int numberColumns = quadratic->getNumCols();
    1657         ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
    1658                                                                    ncolsNow,
    1659                                                                    originalColumn_);
    1660         // and modify linear and check
    1661         double * linear = newObj->linearObjective();
    1662  CoinMemcpyN(presolvedModel_->objective(),ncolsNow,linear);
    1663         int iColumn;
    1664         for ( iColumn=0;iColumn<numberColumns;iColumn++) {
    1665           if (columnQuadraticLength[iColumn])
    1666             mark[iColumn]=1;
    1667         }
    1668         // and new
    1669         quadratic = newObj->quadraticObjective();
    1670         columnQuadraticLength = quadratic->getVectorLengths();
    1671         int numberColumns2 = quadratic->getNumCols();
    1672         for ( iColumn=0;iColumn<numberColumns2;iColumn++) {
    1673           if (columnQuadraticLength[iColumn])
    1674             mark[originalColumn_[iColumn]]=0;
    1675         }
    1676         presolvedModel_->setObjective(newObj);
    1677         delete newObj;
    1678         // final check
    1679         for ( iColumn=0;iColumn<numberColumns;iColumn++)
    1680           if (mark[iColumn])
    1681             printf("Quadratic column %d modified - may be okay\n",iColumn);
    1682         delete [] mark;
    1683       }
    1684 #endif
    1685       delete [] prob.originalColumn_;
    1686       prob.originalColumn_=NULL;
    1687       int nrowsNow = presolvedModel_->getNumRows();
    1688       CoinMemcpyN(prob.originalRow_,nrowsNow,originalRow_);
    1689       delete [] prob.originalRow_;
    1690       prob.originalRow_=NULL;
     1643               ClpQuadraticObjective * quadraticObj = NULL;
     1644               if (originalModel->objectiveAsObject()->type() == 2)
     1645                    quadraticObj = (static_cast< ClpQuadraticObjective*>(originalModel->objectiveAsObject()));
     1646#endif
     1647               if (quadraticObj) {
     1648                    // set up for subset
     1649                    char * mark = new char [ncols_];
     1650                    memset(mark, 0, ncols_);
     1651                    CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
     1652                    //const int * columnQuadratic = quadratic->getIndices();
     1653                    //const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
     1654                    const int * columnQuadraticLength = quadratic->getVectorLengths();
     1655                    //double * quadraticElement = quadratic->getMutableElements();
     1656                    int numberColumns = quadratic->getNumCols();
     1657                    ClpQuadraticObjective * newObj = new ClpQuadraticObjective(*quadraticObj,
     1658                              ncolsNow,
     1659                              originalColumn_);
     1660                    // and modify linear and check
     1661                    double * linear = newObj->linearObjective();
     1662                    CoinMemcpyN(presolvedModel_->objective(), ncolsNow, linear);
     1663                    int iColumn;
     1664                    for ( iColumn = 0; iColumn < numberColumns; iColumn++) {
     1665                         if (columnQuadraticLength[iColumn])
     1666                              mark[iColumn] = 1;
     1667                    }
     1668                    // and new
     1669                    quadratic = newObj->quadraticObjective();
     1670                    columnQuadraticLength = quadratic->getVectorLengths();
     1671                    int numberColumns2 = quadratic->getNumCols();
     1672                    for ( iColumn = 0; iColumn < numberColumns2; iColumn++) {
     1673                         if (columnQuadraticLength[iColumn])
     1674                              mark[originalColumn_[iColumn]] = 0;
     1675                    }
     1676                    presolvedModel_->setObjective(newObj);
     1677                    delete newObj;
     1678                    // final check
     1679                    for ( iColumn = 0; iColumn < numberColumns; iColumn++)
     1680                         if (mark[iColumn])
     1681                              printf("Quadratic column %d modified - may be okay\n", iColumn);
     1682                    delete [] mark;
     1683               }
     1684#endif
     1685               delete [] prob.originalColumn_;
     1686               prob.originalColumn_ = NULL;
     1687               int nrowsNow = presolvedModel_->getNumRows();
     1688               CoinMemcpyN(prob.originalRow_, nrowsNow, originalRow_);
     1689               delete [] prob.originalRow_;
     1690               prob.originalRow_ = NULL;
    16911691#ifndef CLP_NO_STD
    1692       if (!dropNames&&originalModel->lengthNames()) {
    1693         // Redo names
    1694         int iRow;
    1695         std::vector<std::string> rowNames;
    1696         rowNames.reserve(nrowsNow);
    1697         for (iRow=0;iRow<nrowsNow;iRow++) {
    1698           int kRow = originalRow_[iRow];
    1699           rowNames.push_back(originalModel->rowName(kRow));
    1700         }
    1701      
    1702         int iColumn;
    1703         std::vector<std::string> columnNames;
    1704         columnNames.reserve(ncolsNow);
    1705         for (iColumn=0;iColumn<ncolsNow;iColumn++) {
    1706           int kColumn = originalColumn_[iColumn];
    1707           columnNames.push_back(originalModel->columnName(kColumn));
    1708         }
    1709         presolvedModel_->copyNames(rowNames,columnNames);
    1710       } else {
    1711         presolvedModel_->setLengthNames(0);
    1712       }
    1713 #endif
    1714       if (rowObjective_) {
    1715         int iRow;
    1716         int k=-1;
    1717         int nObj=0;
    1718         for (iRow=0;iRow<nrowsNow;iRow++) {
    1719           int kRow = originalRow_[iRow];
    1720           assert (kRow>k);
    1721           k=kRow;
    1722           rowObjective_[iRow]=rowObjective_[kRow];
    1723           if (rowObjective_[iRow])
    1724             nObj++;
    1725         }
    1726         if (nObj) {
    1727           printf("%d costed slacks\n",nObj);
    1728           presolvedModel_->setRowObjective(rowObjective_);
    1729         }
    1730       }
    1731       // now clean up integer variables.  This can modify original
    1732       int i;
    1733       const char * information = presolvedModel_->integerInformation();
    1734       if (information) {
    1735         int numberChanges=0;
    1736         double * lower0 = originalModel_->columnLower();
    1737         double * upper0 = originalModel_->columnUpper();
    1738         double * lower = presolvedModel_->columnLower();
    1739         double * upper = presolvedModel_->columnUpper();
    1740         for (i=0;i<ncolsNow;i++) {
    1741           if (!information[i])
    1742             continue;
    1743           int iOriginal = originalColumn_[i];
    1744           double lowerValue0 = lower0[iOriginal];
    1745           double upperValue0 = upper0[iOriginal];
    1746           double lowerValue = ceil(lower[i]-1.0e-5);
    1747           double upperValue = floor(upper[i]+1.0e-5);
    1748           lower[i]=lowerValue;
    1749           upper[i]=upperValue;
    1750           if (lowerValue>upperValue) {
    1751             numberChanges++;
    1752             presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
    1753                                                        messages)
    1754                                                          <<iOriginal
    1755                                                          <<lowerValue
    1756                                                          <<upperValue
    1757                                                          <<CoinMessageEol;
    1758             result=1;
    1759           } else {
    1760             if (lowerValue>lowerValue0+1.0e-8) {
    1761               lower0[iOriginal] = lowerValue;
    1762               numberChanges++;
    1763             }
    1764             if (upperValue<upperValue0-1.0e-8) {
    1765               upper0[iOriginal] = upperValue;
    1766               numberChanges++;
    1767             }
    1768           }       
    1769         }
    1770         if (numberChanges) {
    1771           presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
    1772                                                      messages)
    1773                                                        <<numberChanges
    1774                                                        <<CoinMessageEol;
    1775           if (!result&&totalPasses>0) {
    1776             result = -1; // round again
    1777             const CoinPresolveAction *paction = paction_;
    1778             while (paction) {
    1779               const CoinPresolveAction *next = paction->next;
    1780               delete paction;
    1781               paction = next;
    1782             }
    1783             paction_=NULL;
    1784           }
    1785         }
    1786       }
    1787     } else if (prob.status_) {
    1788       // infeasible or unbounded
    1789       result=1;
    1790       originalModel->setProblemStatus(prob.status_);
    1791     } else {
    1792       // no changes - model needs restoring after Lou's changes
     1692               if (!dropNames && originalModel->lengthNames()) {
     1693                    // Redo names
     1694                    int iRow;
     1695                    std::vector<std::string> rowNames;
     1696                    rowNames.reserve(nrowsNow);
     1697                    for (iRow = 0; iRow < nrowsNow; iRow++) {
     1698                         int kRow = originalRow_[iRow];
     1699                         rowNames.push_back(originalModel->rowName(kRow));
     1700                    }
     1701
     1702                    int iColumn;
     1703                    std::vector<std::string> columnNames;
     1704                    columnNames.reserve(ncolsNow);
     1705                    for (iColumn = 0; iColumn < ncolsNow; iColumn++) {
     1706                         int kColumn = originalColumn_[iColumn];
     1707                         columnNames.push_back(originalModel->columnName(kColumn));
     1708                    }
     1709                    presolvedModel_->copyNames(rowNames, columnNames);
     1710               } else {
     1711                    presolvedModel_->setLengthNames(0);
     1712               }
     1713#endif
     1714               if (rowObjective_) {
     1715                    int iRow;
     1716                    int k = -1;
     1717                    int nObj = 0;
     1718                    for (iRow = 0; iRow < nrowsNow; iRow++) {
     1719                         int kRow = originalRow_[iRow];
     1720                         assert (kRow > k);
     1721                         k = kRow;
     1722                         rowObjective_[iRow] = rowObjective_[kRow];
     1723                         if (rowObjective_[iRow])
     1724                              nObj++;
     1725                    }
     1726                    if (nObj) {
     1727                         printf("%d costed slacks\n", nObj);
     1728                         presolvedModel_->setRowObjective(rowObjective_);
     1729                    }
     1730               }
     1731               // now clean up integer variables.  This can modify original
     1732               int i;
     1733               const char * information = presolvedModel_->integerInformation();
     1734               if (information) {
     1735                    int numberChanges = 0;
     1736                    double * lower0 = originalModel_->columnLower();
     1737                    double * upper0 = originalModel_->columnUpper();
     1738                    double * lower = presolvedModel_->columnLower();
     1739                    double * upper = presolvedModel_->columnUpper();
     1740                    for (i = 0; i < ncolsNow; i++) {
     1741                         if (!information[i])
     1742                              continue;
     1743                         int iOriginal = originalColumn_[i];
     1744                         double lowerValue0 = lower0[iOriginal];
     1745                         double upperValue0 = upper0[iOriginal];
     1746                         double lowerValue = ceil(lower[i] - 1.0e-5);
     1747                         double upperValue = floor(upper[i] + 1.0e-5);
     1748                         lower[i] = lowerValue;
     1749                         upper[i] = upperValue;
     1750                         if (lowerValue > upperValue) {
     1751                              numberChanges++;
     1752                              presolvedModel_->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
     1753                                        messages)
     1754                                        << iOriginal
     1755                                        << lowerValue
     1756                                        << upperValue
     1757                                        << CoinMessageEol;
     1758                              result = 1;
     1759                         } else {
     1760                              if (lowerValue > lowerValue0 + 1.0e-8) {
     1761                                   lower0[iOriginal] = lowerValue;
     1762                                   numberChanges++;
     1763                              }
     1764                              if (upperValue < upperValue0 - 1.0e-8) {
     1765                                   upper0[iOriginal] = upperValue;
     1766                                   numberChanges++;
     1767                              }
     1768                         }
     1769                    }
     1770                    if (numberChanges) {
     1771                         presolvedModel_->messageHandler()->message(COIN_PRESOLVE_INTEGERMODS,
     1772                                   messages)
     1773                                   << numberChanges
     1774                                   << CoinMessageEol;
     1775                         if (!result && totalPasses > 0) {
     1776                              result = -1; // round again
     1777                              const CoinPresolveAction *paction = paction_;
     1778                              while (paction) {
     1779                                   const CoinPresolveAction *next = paction->next;
     1780                                   delete paction;
     1781                                   paction = next;
     1782                              }
     1783                              paction_ = NULL;
     1784                         }
     1785                    }
     1786               }
     1787          } else if (prob.status_) {
     1788               // infeasible or unbounded
     1789               result = 1;
     1790               originalModel->setProblemStatus(prob.status_);
     1791          } else {
     1792               // no changes - model needs restoring after Lou's changes
    17931793#ifndef CLP_NO_STD
    1794       if (saveFile_=="") {
    1795 #endif
    1796         delete presolvedModel_;
    1797         presolvedModel_ = new ClpSimplex(*originalModel);
    1798         // but we need to remove gaps
    1799         ClpPackedMatrix* clpMatrix =
    1800           dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
    1801         if (clpMatrix) {
    1802           clpMatrix->getPackedMatrix()->removeGaps();
    1803         }
     1794               if (saveFile_ == "") {
     1795#endif
     1796                    delete presolvedModel_;
     1797                    presolvedModel_ = new ClpSimplex(*originalModel);
     1798                    // but we need to remove gaps
     1799                    ClpPackedMatrix* clpMatrix =
     1800                         dynamic_cast< ClpPackedMatrix*>(presolvedModel_->clpMatrix());
     1801                    if (clpMatrix) {
     1802                         clpMatrix->getPackedMatrix()->removeGaps();
     1803                    }
    18041804#ifndef CLP_NO_STD
    1805       } else {
    1806         presolvedModel_=originalModel;
    1807       }
    1808       presolvedModel_->dropNames();
    1809 #endif
    1810      
    1811       // drop integer information if wanted
    1812       if (!keepIntegers)
    1813         presolvedModel_->deleteIntegerInformation();
    1814       result=2;
    1815     }
    1816   }
    1817   if (result==0||result==2) {
    1818     int nrowsAfter = presolvedModel_->getNumRows();
    1819     int ncolsAfter = presolvedModel_->getNumCols();
    1820     CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
    1821     presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
    1822                                                messages)
    1823                                                  <<nrowsAfter<< -(nrows_ - nrowsAfter)
    1824                                                  << ncolsAfter<< -(ncols_ - ncolsAfter)
    1825                                                  <<nelsAfter<< -(nelems_ - nelsAfter)
    1826                                                  <<CoinMessageEol;
    1827   } else {
    1828     destroyPresolve();
    1829     if (presolvedModel_!=originalModel_)
    1830       delete presolvedModel_;
    1831     presolvedModel_=NULL;
    1832   }
    1833   return presolvedModel_;
    1834 }
    1835 
    1836 
     1805               } else {
     1806                    presolvedModel_ = originalModel;
     1807               }
     1808               presolvedModel_->dropNames();
     1809#endif
     1810
     1811               // drop integer information if wanted
     1812               if (!keepIntegers)
     1813                    presolvedModel_->deleteIntegerInformation();
     1814               result = 2;
     1815          }
     1816     }
     1817     if (result == 0 || result == 2) {
     1818          int nrowsAfter = presolvedModel_->getNumRows();
     1819          int ncolsAfter = presolvedModel_->getNumCols();
     1820          CoinBigIndex nelsAfter = presolvedModel_->getNumElements();
     1821          presolvedModel_->messageHandler()->message(COIN_PRESOLVE_STATS,
     1822                    messages)
     1823                    << nrowsAfter << -(nrows_ - nrowsAfter)
     1824                    << ncolsAfter << -(ncols_ - ncolsAfter)
     1825                    << nelsAfter << -(nelems_ - nelsAfter)
     1826                    << CoinMessageEol;
     1827     } else {
     1828          destroyPresolve();
     1829          if (presolvedModel_ != originalModel_)
     1830               delete presolvedModel_;
     1831          presolvedModel_ = NULL;
     1832     }
     1833     return presolvedModel_;
     1834}
     1835
     1836
Note: See TracChangeset for help on using the changeset viewer.