Changeset 1525 for trunk/Clp/src/ClpPresolve.cpp
 Timestamp:
 Feb 26, 2010 12:27:59 PM (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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