Changeset 2385 for trunk/Clp/src/ClpPresolve.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (3 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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