00001
00002
00003 #include <stdio.h>
00004 #include <math.h>
00005
00006 #include "CoinPresolveMatrix.hpp"
00007 #include "CoinPresolveSubst.hpp"
00008 #include "CoinPresolveIsolated.hpp"
00009 #include "CoinPresolveImpliedFree.hpp"
00010 #include "CoinPresolveUseless.hpp"
00011 #include "CoinMessage.hpp"
00012 #include "CoinHelperFunctions.hpp"
00013 #include "CoinSort.hpp"
00014 static const CoinPresolveAction * testRedundant(CoinPresolveMatrix *prob,
00015 const CoinPresolveAction *next,
00016 int & numberInfeasible)
00017 {
00018 prob->pass_++;
00019 numberInfeasible=0;
00020 if (prob->pass_%2!=1)
00021 return next;
00022 int numberColumns = prob->ncols_;
00023 double * columnLower = new double[numberColumns];
00024 double * columnUpper = new double[numberColumns];
00025 memcpy(columnLower,prob->clo_,numberColumns*sizeof(double));
00026 memcpy(columnUpper,prob->cup_,numberColumns*sizeof(double));
00027
00028 const double *element = prob->rowels_;
00029 const int *column = prob->hcol_;
00030 const CoinBigIndex *rowStart = prob->mrstrt_;
00031 const int *rowLength = prob->hinrow_;
00032 int numberRows = prob->nrows_;
00033
00034 int *useless_rows = new int[numberRows];
00035 int nuseless_rows = 0;
00036
00037 double *rowLower = prob->rlo_;
00038 double *rowUpper = prob->rup_;
00039
00040 double tolerance = prob->feasibilityTolerance_;
00041 int numberChanged=1,iPass=0;
00042 double large = 1.0e20;
00043 #ifndef NDEBUG
00044 double large2= 1.0e10*large;
00045 #endif
00046 int totalTightened = 0;
00047
00048 int iRow, iColumn;
00049
00050 int * markRow = new int[numberRows];
00051 for (iRow=0;iRow<numberRows;iRow++)
00052 markRow[iRow]=-1;
00053 #define MAXPASS 10
00054
00055
00056
00057
00058 int numberCheck=-1;
00059 while(numberChanged>numberCheck) {
00060
00061 numberChanged = 0;
00062
00063 if (iPass==MAXPASS) break;
00064 iPass++;
00065
00066 for (iRow = 0; iRow < numberRows; iRow++) {
00067
00068 if ((rowLower[iRow]>-large||rowUpper[iRow]<large)&&markRow[iRow]<0&&rowLength[iRow]>0) {
00069
00070
00071 int infiniteUpper = 0;
00072 int infiniteLower = 0;
00073 double maximumUp = 0.0;
00074 double maximumDown = 0.0;
00075 double newBound;
00076 CoinBigIndex rStart = rowStart[iRow];
00077 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
00078 CoinBigIndex j;
00079
00080
00081 for (j = rStart; j < rEnd; ++j) {
00082 double value=element[j];
00083 iColumn = column[j];
00084 if (value > 0.0) {
00085 if (columnUpper[iColumn] >= large) {
00086 ++infiniteUpper;
00087 } else {
00088 maximumUp += columnUpper[iColumn] * value;
00089 }
00090 if (columnLower[iColumn] <= -large) {
00091 ++infiniteLower;
00092 } else {
00093 maximumDown += columnLower[iColumn] * value;
00094 }
00095 } else if (value<0.0) {
00096 if (columnUpper[iColumn] >= large) {
00097 ++infiniteLower;
00098 } else {
00099 maximumDown += columnUpper[iColumn] * value;
00100 }
00101 if (columnLower[iColumn] <= -large) {
00102 ++infiniteUpper;
00103 } else {
00104 maximumUp += columnLower[iColumn] * value;
00105 }
00106 }
00107 }
00108
00109 maximumUp += 1.0e-8*fabs(maximumUp);
00110 maximumDown -= 1.0e-8*fabs(maximumDown);
00111 double maxUp = maximumUp+infiniteUpper*1.0e31;
00112 double maxDown = maximumDown-infiniteLower*1.0e31;
00113 if (maxUp <= rowUpper[iRow] + tolerance &&
00114 maxDown >= rowLower[iRow] - tolerance) {
00115
00116 } else {
00117 if (maxUp < rowLower[iRow] -100.0*tolerance ||
00118 maxDown > rowUpper[iRow]+100.0*tolerance) {
00119
00120 numberInfeasible++;
00121 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00122 prob->messages())
00123 <<iRow
00124 <<rowLower[iRow]
00125 <<rowUpper[iRow]
00126 <<CoinMessageEol;
00127 break;
00128 }
00129 double lower = rowLower[iRow];
00130 double upper = rowUpper[iRow];
00131 for (j = rStart; j < rEnd; ++j) {
00132 double value=element[j];
00133 iColumn = column[j];
00134 double nowLower = columnLower[iColumn];
00135 double nowUpper = columnUpper[iColumn];
00136 if (value > 0.0) {
00137
00138 if (lower>-large) {
00139 if (!infiniteUpper) {
00140 assert(nowUpper < large2);
00141 newBound = nowUpper +
00142 (lower - maximumUp) / value;
00143
00144 if (fabs(maximumUp)>1.0e8)
00145 newBound -= 1.0e-12*fabs(maximumUp);
00146 } else if (infiniteUpper==1&&nowUpper>large) {
00147 newBound = (lower -maximumUp) / value;
00148
00149 if (fabs(maximumUp)>1.0e8)
00150 newBound -= 1.0e-12*fabs(maximumUp);
00151 } else {
00152 newBound = -COIN_DBL_MAX;
00153 }
00154 if (newBound > nowLower + 1.0e-12&&newBound>-large) {
00155
00156 columnLower[iColumn] = newBound;
00157 markRow[iRow]=1;
00158 numberChanged++;
00159
00160 if (nowUpper - newBound <
00161 -100.0*tolerance) {
00162 numberInfeasible++;
00163 }
00164
00165 double now;
00166 if (nowLower<-large) {
00167 now=0.0;
00168 infiniteLower--;
00169 } else {
00170 now = nowLower;
00171 }
00172 maximumDown += (newBound-now) * value;
00173 nowLower = newBound;
00174 }
00175 }
00176 if (upper <large) {
00177 if (!infiniteLower) {
00178 assert(nowLower >- large2);
00179 newBound = nowLower +
00180 (upper - maximumDown) / value;
00181
00182 if (fabs(maximumDown)>1.0e8)
00183 newBound += 1.0e-12*fabs(maximumDown);
00184 } else if (infiniteLower==1&&nowLower<-large) {
00185 newBound = (upper - maximumDown) / value;
00186
00187 if (fabs(maximumDown)>1.0e8)
00188 newBound += 1.0e-12*fabs(maximumDown);
00189 } else {
00190 newBound = COIN_DBL_MAX;
00191 }
00192 if (newBound < nowUpper - 1.0e-12&&newBound<large) {
00193
00194 columnUpper[iColumn] = newBound;
00195 markRow[iRow]=1;
00196 numberChanged++;
00197
00198 if (newBound - nowLower <
00199 -100.0*tolerance) {
00200 numberInfeasible++;
00201 }
00202
00203 double now;
00204 if (nowUpper>large) {
00205 now=0.0;
00206 infiniteUpper--;
00207 } else {
00208 now = nowUpper;
00209 }
00210 maximumUp += (newBound-now) * value;
00211 nowUpper = newBound;
00212 }
00213 }
00214 } else {
00215
00216 if (lower>-large) {
00217 if (!infiniteUpper) {
00218 assert(nowLower < large2);
00219 newBound = nowLower +
00220 (lower - maximumUp) / value;
00221
00222 if (fabs(maximumUp)>1.0e8)
00223 newBound += 1.0e-12*fabs(maximumUp);
00224 } else if (infiniteUpper==1&&nowLower<-large) {
00225 newBound = (lower -maximumUp) / value;
00226
00227 if (fabs(maximumUp)>1.0e8)
00228 newBound += 1.0e-12*fabs(maximumUp);
00229 } else {
00230 newBound = COIN_DBL_MAX;
00231 }
00232 if (newBound < nowUpper - 1.0e-12&&newBound<large) {
00233
00234 columnUpper[iColumn] = newBound;
00235 markRow[iRow]=1;
00236 numberChanged++;
00237
00238 if (newBound - nowLower <
00239 -100.0*tolerance) {
00240 numberInfeasible++;
00241 }
00242
00243 double now;
00244 if (nowUpper>large) {
00245 now=0.0;
00246 infiniteLower--;
00247 } else {
00248 now = nowUpper;
00249 }
00250 maximumDown += (newBound-now) * value;
00251 nowUpper = newBound;
00252 }
00253 }
00254 if (upper <large) {
00255 if (!infiniteLower) {
00256 assert(nowUpper < large2);
00257 newBound = nowUpper +
00258 (upper - maximumDown) / value;
00259
00260 if (fabs(maximumDown)>1.0e8)
00261 newBound -= 1.0e-12*fabs(maximumDown);
00262 } else if (infiniteLower==1&&nowUpper>large) {
00263 newBound = (upper - maximumDown) / value;
00264
00265 if (fabs(maximumDown)>1.0e8)
00266 newBound -= 1.0e-12*fabs(maximumDown);
00267 } else {
00268 newBound = -COIN_DBL_MAX;
00269 }
00270 if (newBound > nowLower + 1.0e-12&&newBound>-large) {
00271
00272 columnLower[iColumn] = newBound;
00273 markRow[iRow]=1;
00274 numberChanged++;
00275
00276 if (nowUpper - newBound <
00277 -100.0*tolerance) {
00278 numberInfeasible++;
00279 }
00280
00281 double now;
00282 if (nowLower<-large) {
00283 now=0.0;
00284 infiniteUpper--;
00285 } else {
00286 now = nowLower;
00287 }
00288 maximumUp += (newBound-now) * value;
00289 nowLower = newBound;
00290 }
00291 }
00292 }
00293 }
00294 }
00295 }
00296 }
00297 totalTightened += numberChanged;
00298 if (iPass==1)
00299 numberCheck=numberChanged>>4;
00300 if (numberInfeasible) break;
00301 }
00302 if (!numberInfeasible) {
00303 for (iRow = 0; iRow < numberRows; iRow++) {
00304
00305 if ((rowLower[iRow]>-large||rowUpper[iRow]<large)&&markRow[iRow]<0&&rowLength[iRow]>0) {
00306
00307
00308 int infiniteUpper = 0;
00309 int infiniteLower = 0;
00310 double maximumUp = 0.0;
00311 double maximumDown = 0.0;
00312 CoinBigIndex rStart = rowStart[iRow];
00313 CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
00314 CoinBigIndex j;
00315
00316
00317 for (j = rStart; j < rEnd; ++j) {
00318 double value=element[j];
00319 iColumn = column[j];
00320 if (value > 0.0) {
00321 if (columnUpper[iColumn] >= large) {
00322 ++infiniteUpper;
00323 } else {
00324 maximumUp += columnUpper[iColumn] * value;
00325 }
00326 if (columnLower[iColumn] <= -large) {
00327 ++infiniteLower;
00328 } else {
00329 maximumDown += columnLower[iColumn] * value;
00330 }
00331 } else if (value<0.0) {
00332 if (columnUpper[iColumn] >= large) {
00333 ++infiniteLower;
00334 } else {
00335 maximumDown += columnUpper[iColumn] * value;
00336 }
00337 if (columnLower[iColumn] <= -large) {
00338 ++infiniteUpper;
00339 } else {
00340 maximumUp += columnLower[iColumn] * value;
00341 }
00342 }
00343 }
00344
00345 maximumUp += 1.0e-8*fabs(maximumUp);
00346 maximumDown -= 1.0e-8*fabs(maximumDown);
00347 double maxUp = maximumUp+infiniteUpper*1.0e31;
00348 double maxDown = maximumDown-infiniteLower*1.0e31;
00349 if (maxUp <= rowUpper[iRow] + tolerance &&
00350 maxDown >= rowLower[iRow] - tolerance) {
00351
00352
00353 useless_rows[nuseless_rows++] = iRow;
00354 prob->addRow(iRow);
00355
00356 }
00357 }
00358 }
00359 }
00360 if (nuseless_rows)
00361 next = useless_constraint_action::presolve(prob,
00362 useless_rows, nuseless_rows,
00363 next);
00364
00365 delete[]useless_rows;
00366
00367 delete [] columnLower;
00368 delete [] columnUpper;
00369 delete [] markRow;
00370 return next;
00371 }
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405 const CoinPresolveAction *implied_free_action::presolve(CoinPresolveMatrix *prob,
00406 const CoinPresolveAction *next,
00407 int & fill_level)
00408 {
00409 double *colels = prob->colels_;
00410 int *hrow = prob->hrow_;
00411 const CoinBigIndex *mcstrt = prob->mcstrt_;
00412 int *hincol = prob->hincol_;
00413 const int ncols = prob->ncols_;
00414
00415 const double *clo = prob->clo_;
00416 const double *cup = prob->cup_;
00417
00418 const double *rowels = prob->rowels_;
00419 const int *hcol = prob->hcol_;
00420 const CoinBigIndex *mrstrt = prob->mrstrt_;
00421 int *hinrow = prob->hinrow_;
00422 int nrows = prob->nrows_;
00423
00424 double *rlo = prob->rlo_;
00425 double *rup = prob->rup_;
00426
00427 double *cost = prob->cost_;
00428
00429 presolvehlink *rlink = prob->rlink_;
00430 presolvehlink *clink = prob->clink_;
00431
00432 const char *integerType = prob->integerType_;
00433
00434 const double tol = prob->feasibilityTolerance_;
00435 #if 1
00436
00437 int numberInfeasible;
00438 next = testRedundant(prob,next,numberInfeasible);
00439 if (numberInfeasible) {
00440
00441 prob->status_|= 1;
00442 return (next);
00443 }
00444 #endif
00445
00446
00447 action *actions = new action [ncols];
00448 int nactions = 0;
00449
00450 int *implied_free = new int[ncols];
00451 int i;
00452 for (i=0;i<ncols;i++)
00453 implied_free[i]=-1;
00454
00455
00456 int * infiniteDown = new int [nrows];
00457 int * infiniteUp = new int [nrows];
00458 double * maxDown = new double[nrows];
00459 double * maxUp = new double[nrows];
00460
00461
00462
00463 for (i=0;i<nrows;i++) {
00464 if (hinrow[i]>1)
00465 infiniteUp[i]=-1;
00466 else
00467 infiniteUp[i]=-2;
00468 }
00469 double large=1.0e20;
00470
00471 int numberLook = prob->numberColsToDo_;
00472 int iLook;
00473 int * look = prob->colsToDo_;
00474 int * look2 = NULL;
00475
00476 if (fill_level<0) {
00477 look2 = new int[ncols];
00478 look=look2;
00479 for (iLook=0;iLook<ncols;iLook++)
00480 look[iLook]=iLook;
00481 numberLook=ncols;
00482 }
00483
00484 for (iLook=0;iLook<numberLook;iLook++) {
00485 int j=look[iLook];
00486 if ((hincol[j] <= 3) && !integerType[j]) {
00487 if (hincol[j]>1) {
00488
00489 CoinBigIndex kcs = mcstrt[j];
00490 CoinBigIndex kce = kcs + hincol[j];
00491 bool possible = false;
00492 bool singleton = false;
00493 CoinBigIndex k;
00494 double largestElement=0.0;
00495 for (k=kcs; k<kce; ++k) {
00496 int row = hrow[k];
00497 double coeffj = colels[k];
00498
00499
00500 if (hinrow[row] > 1 ) {
00501 if ( fabs(rlo[row] - rup[row]) < tol &&
00502 fabs(coeffj) > ZTOLDP) {
00503 possible=true;
00504 }
00505 largestElement = max(largestElement,fabs(coeffj));
00506 } else {
00507 singleton=true;
00508 }
00509 }
00510 if (possible&&!singleton) {
00511 double low=-COIN_DBL_MAX;
00512 double high=COIN_DBL_MAX;
00513
00514 for (k=kcs; k<kce; ++k) {
00515 int row = hrow[k];
00516 double coeffj = colels[k];
00517 if (fabs(coeffj) > ZTOLDP) {
00518 if (infiniteUp[row]==-1) {
00519
00520 CoinBigIndex krs = mrstrt[row];
00521 CoinBigIndex kre = krs + hinrow[row];
00522 int infiniteUpper = 0;
00523 int infiniteLower = 0;
00524 double maximumUp = 0.0;
00525 double maximumDown = 0.0;
00526 CoinBigIndex kk;
00527
00528 for (kk = krs; kk < kre; ++kk) {
00529 double value=rowels[kk];
00530 int iColumn = hcol[kk];
00531 if (value > 0.0) {
00532 if (cup[iColumn] >= large) {
00533 ++infiniteUpper;
00534 } else {
00535 maximumUp += cup[iColumn] * value;
00536 }
00537 if (clo[iColumn] <= -large) {
00538 ++infiniteLower;
00539 } else {
00540 maximumDown += clo[iColumn] * value;
00541 }
00542 } else if (value<0.0) {
00543 if (cup[iColumn] >= large) {
00544 ++infiniteLower;
00545 } else {
00546 maximumDown += cup[iColumn] * value;
00547 }
00548 if (clo[iColumn] <= -large) {
00549 ++infiniteUpper;
00550 } else {
00551 maximumUp += clo[iColumn] * value;
00552 }
00553 }
00554 }
00555 double maxUpx = maximumUp+infiniteUpper*1.0e31;
00556 double maxDownx = maximumDown-infiniteLower*1.0e31;
00557 if (maxUpx <= rup[row] + tol &&
00558 maxDownx >= rlo[row] - tol) {
00559
00560
00561 infiniteUp[row]=-3;
00562
00563 } else if (maxUpx < rlo[row] -tol ) {
00564
00565 prob->status_|= 1;
00566 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00567 prob->messages())
00568 <<row
00569 <<rlo[row]
00570 <<rup[row]
00571 <<CoinMessageEol;
00572 infiniteUp[row]=-3;
00573 break;
00574 } else if ( maxDownx > rup[row]+tol) {
00575
00576 prob->status_|= 1;
00577 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00578 prob->messages())
00579 <<row
00580 <<rlo[row]
00581 <<rup[row]
00582 <<CoinMessageEol;
00583 infiniteUp[row]=-3;
00584 break;
00585 } else {
00586 infiniteUp[row]=infiniteUpper;
00587 infiniteDown[row]=infiniteLower;
00588 maxUp[row]=maximumUp;
00589 maxDown[row]=maximumDown;
00590 }
00591 }
00592 if (infiniteUp[row]>=0) {
00593 double lower = rlo[row];
00594 double upper = rup[row];
00595 double value=coeffj;
00596 double nowLower = clo[j];
00597 double nowUpper = cup[j];
00598 double newBound;
00599 int infiniteUpper=infiniteUp[row];
00600 int infiniteLower=infiniteDown[row];
00601 double maximumUp = maxUp[row];
00602 double maximumDown = maxDown[row];
00603 if (value > 0.0) {
00604
00605 if (lower>-large) {
00606 if (!infiniteUpper) {
00607 assert(nowUpper < large);
00608 newBound = nowUpper +
00609 (lower - maximumUp) / value;
00610
00611 if (fabs(maximumUp)>1.0e8)
00612 newBound -= 1.0e-12*fabs(maximumUp);
00613 } else if (infiniteUpper==1&&nowUpper>large) {
00614 newBound = (lower -maximumUp) / value;
00615
00616 if (fabs(maximumUp)>1.0e8)
00617 newBound -= 1.0e-12*fabs(maximumUp);
00618 } else {
00619 newBound = -COIN_DBL_MAX;
00620 }
00621 if (newBound > nowLower + 1.0e-12) {
00622
00623
00624 double now;
00625 if (nowLower<-large) {
00626 now=0.0;
00627 infiniteLower--;
00628 } else {
00629 now = nowLower;
00630 }
00631 maximumDown += (newBound-now) * value;
00632 nowLower = newBound;
00633 }
00634 low=max(low,newBound);
00635 }
00636 if (upper <large) {
00637 if (!infiniteLower) {
00638 assert(nowLower >- large);
00639 newBound = nowLower +
00640 (upper - maximumDown) / value;
00641
00642 if (fabs(maximumDown)>1.0e8)
00643 newBound += 1.0e-12*fabs(maximumDown);
00644 } else if (infiniteLower==1&&nowLower<-large) {
00645 newBound = (upper - maximumDown) / value;
00646
00647 if (fabs(maximumDown)>1.0e8)
00648 newBound += 1.0e-12*fabs(maximumDown);
00649 } else {
00650 newBound = COIN_DBL_MAX;
00651 }
00652 if (newBound < nowUpper - 1.0e-12) {
00653
00654
00655 double now;
00656 if (nowUpper>large) {
00657 now=0.0;
00658 infiniteUpper--;
00659 } else {
00660 now = nowUpper;
00661 }
00662 maximumUp += (newBound-now) * value;
00663 nowUpper = newBound;
00664 }
00665 high=min(high,newBound);
00666 }
00667 } else {
00668
00669 if (lower>-large) {
00670 if (!infiniteUpper) {
00671 assert(nowLower >- large);
00672 newBound = nowLower +
00673 (lower - maximumUp) / value;
00674
00675 if (fabs(maximumUp)>1.0e8)
00676 newBound += 1.0e-12*fabs(maximumUp);
00677 } else if (infiniteUpper==1&&nowLower<-large) {
00678 newBound = (lower -maximumUp) / value;
00679
00680 if (fabs(maximumUp)>1.0e8)
00681 newBound += 1.0e-12*fabs(maximumUp);
00682 } else {
00683 newBound = COIN_DBL_MAX;
00684 }
00685 if (newBound < nowUpper - 1.0e-12) {
00686
00687
00688 double now;
00689 if (nowUpper>large) {
00690 now=0.0;
00691 infiniteLower--;
00692 } else {
00693 now = nowUpper;
00694 }
00695 maximumDown += (newBound-now) * value;
00696 nowUpper = newBound;
00697 }
00698 high=min(high,newBound);
00699 }
00700 if (upper <large) {
00701 if (!infiniteLower) {
00702 assert(nowUpper < large);
00703 newBound = nowUpper +
00704 (upper - maximumDown) / value;
00705
00706 if (fabs(maximumDown)>1.0e8)
00707 newBound -= 1.0e-12*fabs(maximumDown);
00708 } else if (infiniteLower==1&&nowUpper>large) {
00709 newBound = (upper - maximumDown) / value;
00710
00711 if (fabs(maximumDown)>1.0e8)
00712 newBound -= 1.0e-12*fabs(maximumDown);
00713 } else {
00714 newBound = -COIN_DBL_MAX;
00715 }
00716 if (newBound > nowLower + 1.0e-12) {
00717
00718
00719 double now;
00720 if (nowLower<-large) {
00721 now=0.0;
00722 infiniteUpper--;
00723 } else {
00724 now = nowLower;
00725 }
00726 maximumUp += (newBound-now) * value;
00727 nowLower = newBound;
00728 }
00729 low = max(low,newBound);
00730 }
00731 }
00732 } else if (infiniteUp[row]==-3) {
00733
00734 high=COIN_DBL_MAX;
00735 low=-COIN_DBL_MAX;
00736 break;
00737 }
00738 }
00739 }
00740 if (clo[j] <= low && high <= cup[j]) {
00741
00742
00743
00744 largestElement *= 0.1;
00745 int krow=-1;
00746 int ninrow=ncols+1;
00747 for (k=kcs; k<kce; ++k) {
00748 int row = hrow[k];
00749 double coeffj = colels[k];
00750 if ( fabs(rlo[row] - rup[row]) < tol &&
00751 fabs(coeffj) > largestElement) {
00752 if (hinrow[row]<ninrow) {
00753 ninrow=hinrow[row];
00754 krow=row;
00755 }
00756 }
00757 }
00758 if (krow>=0) {
00759 implied_free[j] = krow;
00760
00761 infiniteUp[krow]=-3;
00762
00763
00764 }
00765 }
00766 }
00767 } else if (hincol[j]) {
00768
00769 CoinBigIndex k = mcstrt[j];
00770 int row = hrow[k];
00771 double coeffj = colels[k];
00772 if ((!cost[j]||rlo[row]==rup[row])&&hinrow[row]>1&&
00773 fabs(coeffj) > ZTOLDP&&infiniteUp[row]!=-3) {
00774
00775 CoinBigIndex krs = mrstrt[row];
00776 CoinBigIndex kre = krs + hinrow[row];
00777
00778 double maxup, maxdown, ilow, iup;
00779 implied_bounds(rowels, clo, cup, hcol,
00780 krs, kre,
00781 &maxup, &maxdown,
00782 j, rlo[row], rup[row], &ilow, &iup);
00783
00784
00785 if (maxup < PRESOLVE_INF && maxup + tol < rlo[row]) {
00786
00787 prob->status_|= 1;
00788 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00789 prob->messages())
00790 <<row
00791 <<rlo[row]
00792 <<rup[row]
00793 <<CoinMessageEol;
00794 break;
00795 } else if (-PRESOLVE_INF < maxdown && rup[row] < maxdown - tol) {
00796
00797 prob->status_|= 1;
00798 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00799 prob->messages())
00800 <<row
00801 <<rlo[row]
00802 <<rup[row]
00803 <<CoinMessageEol;
00804 break;
00805 } else if (clo[j] <= ilow && iup <= cup[j]) {
00806
00807
00808 implied_free[j] = row;
00809 infiniteUp[row]=-3;
00810
00811
00812 }
00813 }
00814 }
00815 }
00816 }
00817
00818
00819 delete [] infiniteDown;
00820 delete [] infiniteUp;
00821 delete [] maxDown;
00822 delete [] maxUp;
00823
00824 int isolated_row = -1;
00825
00826
00827
00828
00829
00830 for (iLook=0;iLook<numberLook;iLook++) {
00831 int j=look[iLook];
00832 if (hincol[j] == 1 && implied_free[j] >=0) {
00833 CoinBigIndex kcs = mcstrt[j];
00834 int row = hrow[kcs];
00835 double coeffj = colels[kcs];
00836
00837 CoinBigIndex krs = mrstrt[row];
00838 CoinBigIndex kre = krs + hinrow[row];
00839
00840
00841
00842 {
00843 int n = 0;
00844 for (CoinBigIndex k=krs; k<kre; ++k)
00845 n += hincol[hcol[k]];
00846 if (n==hinrow[row]) {
00847 isolated_row = row;
00848 break;
00849 }
00850 }
00851
00852 const bool nonzero_cost = (cost[j] != 0.0&&fabs(rup[row]-rlo[row])<=tol);
00853
00854 double *save_costs = nonzero_cost ? new double[hinrow[row]] : NULL;
00855
00856 {
00857 action *s = &actions[nactions++];
00858
00859 s->row = row;
00860 s->col = j;
00861
00862 s->clo = clo[j];
00863 s->cup = cup[j];
00864 s->rlo = rlo[row];
00865 s->rup = rup[row];
00866
00867 s->ninrow = hinrow[row];
00868 s->rowels = presolve_duparray(&rowels[krs], hinrow[row]);
00869 s->rowcols = presolve_duparray(&hcol[krs], hinrow[row]);
00870 s->costs = save_costs;
00871 }
00872
00873 if (nonzero_cost) {
00874 double rhs = rlo[row];
00875 double costj = cost[j];
00876
00877 #if DEBUG_PRESOLVE
00878 printf("FREE COSTS: %g ", costj);
00879 #endif
00880 for (CoinBigIndex k=krs; k<kre; k++) {
00881 int jcol = hcol[k];
00882 save_costs[k-krs] = cost[jcol];
00883
00884 if (jcol != j) {
00885 double coeff = rowels[k];
00886
00887 #if DEBUG_PRESOLVE
00888 printf("%g %g ", cost[jcol], coeff/coeffj);
00889 #endif
00890
00891
00892
00893
00894
00895 cost[jcol] += costj * (-coeff / coeffj);
00896 }
00897 }
00898 #if DEBUG_PRESOLVE
00899 printf("\n");
00900
00901
00902 printf("BIAS??????? %g %g %g %g\n",
00903 costj * rhs / coeffj,
00904 costj, rhs, coeffj);
00905 #endif
00906 prob->change_bias(costj * rhs / coeffj);
00907
00908 cost[j] = 0.0;
00909 }
00910
00911
00912 for (CoinBigIndex k=krs; k<kre; k++) {
00913 int jcol=hcol[k];
00914 prob->addCol(jcol);
00915 presolve_delete_from_row(jcol, row, mcstrt, hincol, hrow, colels);
00916 }
00917 PRESOLVE_REMOVE_LINK(rlink, row);
00918 hinrow[row] = 0;
00919
00920
00921 rlo[row] = 0.0;
00922 rup[row] = 0.0;
00923
00924 PRESOLVE_REMOVE_LINK(clink, j);
00925 hincol[j] = 0;
00926
00927 implied_free[j] = -1;
00928 }
00929 }
00930
00931 delete [] look2;
00932 if (nactions) {
00933 #if PRESOLVE_SUMMARY
00934 printf("NIMPLIED FREE: %d\n", nactions);
00935 #endif
00936 action *actions1 = new action[nactions];
00937 CoinMemcpyN(actions, nactions, actions1);
00938 next = new implied_free_action(nactions, actions1, next);
00939 }
00940 delete [] actions;
00941
00942 if (isolated_row != -1) {
00943 const CoinPresolveAction *nextX = isolated_constraint_action::presolve(prob,
00944 isolated_row, next);
00945 if (nextX)
00946 next = nextX;
00947 }
00948
00949 if (fill_level) {
00950 next = subst_constraint_action::presolve(prob, implied_free, next,fill_level);
00951 }
00952 delete[]implied_free;
00953
00954 return (next);
00955 }
00956
00957
00958
00959 const char *implied_free_action::name() const
00960 {
00961 return ("implied_free_action");
00962 }
00963
00964 void implied_free_action::postsolve(CoinPostsolveMatrix *prob) const
00965 {
00966 const action *const actions = actions_;
00967 const int nactions = nactions_;
00968
00969 double *elementByColumn = prob->colels_;
00970 int *hrow = prob->hrow_;
00971 CoinBigIndex *columnStart = prob->mcstrt_;
00972 int *numberInColumn = prob->hincol_;
00973 int *link = prob->link_;
00974
00975 double *clo = prob->clo_;
00976 double *cup = prob->cup_;
00977
00978 double *rlo = prob->rlo_;
00979 double *rup = prob->rup_;
00980
00981 double *sol = prob->sol_;
00982
00983 double *rcosts = prob->rcosts_;
00984 double *dcost = prob->cost_;
00985
00986 double *acts = prob->acts_;
00987 double *rowduals = prob->rowduals_;
00988
00989
00990
00991 const double maxmin = prob->maxmin_;
00992
00993 char *cdone = prob->cdone_;
00994 char *rdone = prob->rdone_;
00995 CoinBigIndex free_list = prob->free_list_;
00996
00997 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
00998
00999 int irow = f->row;
01000 int icol = f->col;
01001
01002 int ninrow = f->ninrow;
01003 const double *rowels = f->rowels;
01004 const int *rowcols = f->rowcols;
01005 const double *save_costs = f->costs;
01006
01007
01008
01009 {
01010 for (int k = 0; k<ninrow; k++) {
01011 int jcol = rowcols[k];
01012 double coeff = rowels[k];
01013
01014 if (save_costs) {
01015 rcosts[jcol] += maxmin*(save_costs[k]-dcost[jcol]);
01016 dcost[jcol] = save_costs[k];
01017 }
01018 {
01019 CoinBigIndex kk = free_list;
01020 free_list = link[free_list];
01021
01022 check_free_list(free_list);
01023
01024 link[kk] = columnStart[jcol];
01025 columnStart[jcol] = kk;
01026 elementByColumn[kk] = coeff;
01027 hrow[kk] = irow;
01028 }
01029
01030 if (jcol == icol) {
01031
01032 numberInColumn[jcol] = 1;
01033 clo[icol] = f->clo;
01034 cup[icol] = f->cup;
01035
01036 cdone[icol] = IMPLIED_FREE;
01037 } else {
01038 numberInColumn[jcol]++;
01039 }
01040 }
01041 rdone[irow] = IMPLIED_FREE;
01042
01043 rlo[irow] = f->rlo;
01044 rup[irow] = f->rup;
01045 }
01046 deleteAction( save_costs,double*);
01047
01048
01049
01050 {
01051 double act = 0.0;
01052 double coeff = 0.0;
01053
01054 for (int k = 0; k<ninrow; k++)
01055 if (rowcols[k] == icol)
01056 coeff = rowels[k];
01057 else {
01058 int jcol = rowcols[k];
01059 PRESOLVE_STMT(CoinBigIndex kk = presolve_find_row2(irow, columnStart[jcol], numberInColumn[jcol], hrow, link));
01060 act += rowels[k] * sol[jcol];
01061 }
01062
01063 PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
01064 double thisCost = maxmin*dcost[icol];
01065 double loActivity,upActivity;
01066 if (coeff>0) {
01067 loActivity = (rlo[irow]-act)/coeff;
01068 upActivity = (rup[irow]-act)/coeff;
01069 } else {
01070 loActivity = (rup[irow]-act)/coeff;
01071 upActivity = (rlo[irow]-act)/coeff;
01072 }
01073 loActivity = max(loActivity,clo[icol]);
01074 upActivity = min(upActivity,cup[icol]);
01075 int where;
01076 const double tolCheck = 0.1*prob->ztolzb_;
01077 if (loActivity<clo[icol]+tolCheck/fabs(coeff)&&thisCost>=0.0)
01078 where=-1;
01079 else if (upActivity>cup[icol]-tolCheck/fabs(coeff)&&thisCost<0.0)
01080 where=1;
01081 else
01082 where =0;
01083
01084 double possibleDual = thisCost/coeff;
01085 if (where) {
01086 double worst= prob->ztoldj_;
01087 for (int k = 0; k<ninrow; k++) {
01088 int jcol = rowcols[k];
01089 if (jcol!=icol) {
01090 CoinPrePostsolveMatrix::Status status = prob->getColumnStatus(jcol);
01091
01092 if (status==CoinPrePostsolveMatrix::basic) {
01093 if (fabs(rcosts[jcol])>worst)
01094 worst=fabs(rcosts[jcol]);
01095 } else if (sol[jcol]<clo[jcol]+ZTOLDP) {
01096 if (-rcosts[jcol]>worst)
01097 worst=-rcosts[jcol];
01098 } else if (sol[jcol]>cup[jcol]-ZTOLDP) {
01099 if (rcosts[jcol]>worst)
01100 worst=rcosts[jcol];
01101 }
01102 }
01103 }
01104 if (worst>prob->ztoldj_) {
01105
01106 double worst2 = prob->ztoldj_;
01107 for (int k = 0; k<ninrow; k++) {
01108 int jcol = rowcols[k];
01109 if (jcol!=icol) {
01110 double coeff = rowels[k];
01111 double newDj = rcosts[jcol]-possibleDual*coeff;
01112 CoinPrePostsolveMatrix::Status status = prob->getColumnStatus(jcol);
01113
01114 if (status==CoinPrePostsolveMatrix::basic) {
01115 if (fabs(newDj)>worst2)
01116 worst2=fabs(newDj);
01117 } else if (sol[jcol]<clo[jcol]+ZTOLDP) {
01118 if (-newDj>worst2)
01119 worst2=-newDj;
01120 } else if (sol[jcol]>cup[jcol]-ZTOLDP) {
01121 if (newDj>worst2)
01122 worst2=newDj;
01123 }
01124 }
01125 }
01126 if (worst2<worst)
01127 where=0;
01128 }
01129 }
01130 if (!where) {
01131
01132 rowduals[irow] = possibleDual;
01133 if ((rlo[irow] < rup[irow] && rowduals[irow] < 0.0)
01134 || rlo[irow]< -1.0e20) {
01135 if (rlo[irow]<-1.0e20&&rowduals[irow]>ZTOLDP)
01136 printf("IMP %g %g %g\n",rlo[irow],rup[irow],rowduals[irow]);
01137 sol[icol] = (rup[irow] - act) / coeff;
01138 assert (sol[icol]>=clo[icol]-1.0e-5&&sol[icol]<=cup[icol]+1.0e-5);
01139 acts[irow] = rup[irow];
01140 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atUpperBound);
01141 } else {
01142 sol[icol] = (rlo[irow] - act) / coeff;
01143 assert (sol[icol]>=clo[icol]-1.0e-5&&sol[icol]<=cup[icol]+1.0e-5);
01144 acts[irow] = rlo[irow];
01145 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
01146 }
01147 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::basic);
01148 for (int k = 0; k<ninrow; k++) {
01149 int jcol = rowcols[k];
01150 double coeff = rowels[k];
01151 rcosts[jcol] -= possibleDual*coeff;
01152 }
01153 rcosts[icol] = 0.0;
01154 } else {
01155 rowduals[irow] = 0.0;
01156 rcosts[icol] = thisCost;
01157 prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
01158 if (where<0) {
01159
01160 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
01161 sol[icol]=clo[icol];
01162 } else {
01163
01164 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
01165 sol[icol]=cup[icol];
01166 }
01167 acts[irow] = act + sol[icol]*coeff;
01168 assert (acts[irow]>=rlo[irow]-1.0e-5&&acts[irow]<=rup[irow]+1.0e-5);
01169 }
01170 #if DEBUG_PRESOLVE
01171 {
01172 double *colels = prob->colels_;
01173 int *hrow = prob->hrow_;
01174 const CoinBigIndex *mcstrt = prob->mcstrt_;
01175 int *hincol = prob->hincol_;
01176 for (int j = 0; j<ninrow; j++) {
01177 int jcol = rowcols[j];
01178 CoinBigIndex k = mcstrt[jcol];
01179 int nx = hincol[jcol];
01180 double dj = dcost[jcol];
01181 for (int i=0; i<nx; ++i) {
01182 int row = hrow[k];
01183 double coeff = colels[k];
01184 k = link[k];
01185 dj -= rowduals[row] * coeff;
01186
01187
01188 }
01189 if (fabs(dj-rcosts[jcol])>1.0e-3)
01190 printf("changed\n");
01191 }
01192 }
01193 #endif
01194 }
01195 }
01196 prob->free_list_ = free_list;
01197 }
01198 implied_free_action::~implied_free_action()
01199 {
01200 int i;
01201 for (i=0;i<nactions_;i++) {
01202
01203
01204 deleteAction(actions_[i].rowcols,int *);
01205 deleteAction(actions_[i].rowels,int *);
01206
01207 }
01208
01209 deleteAction(actions_,action *);
01210 }