00001
00002
00003 #include <stdio.h>
00004 #include <math.h>
00005
00006 #include "CoinPresolveMatrix.hpp"
00007 #include "CoinPresolveFixed.hpp"
00008 #include "CoinPresolveDual.hpp"
00009 #include "CoinMessage.hpp"
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 const CoinPresolveAction *remove_dual_action::presolve(CoinPresolveMatrix *prob,
00027 const CoinPresolveAction *next)
00028 {
00029 double *colels = prob->colels_;
00030 int *hrow = prob->hrow_;
00031 CoinBigIndex *mcstrt = prob->mcstrt_;
00032 int *hincol = prob->hincol_;
00033 int ncols = prob->ncols_;
00034
00035 double *clo = prob->clo_;
00036 double *cup = prob->cup_;
00037
00038
00039
00040
00041
00042 int nrows = prob->nrows_;
00043
00044 double *rlo = prob->rlo_;
00045 double *rup = prob->rup_;
00046
00047 double *dcost = prob->cost_;
00048
00049 const double maxmin = prob->maxmin_;
00050
00051 const double ekkinf = 1e28;
00052 const double ekkinf2 = 1e20;
00053 const double ztoldj = prob->ztoldj_;
00054
00055 double *rdmin = new double[nrows];
00056 double *rdmax = new double[nrows];
00057
00058
00059
00060
00061
00062 int i;
00063 for ( i = 0; i < nrows; i++) {
00064 double sup = -rlo[i];
00065 double slo = -rup[i];
00066 bool no_lb = (slo <= -ekkinf);
00067 bool no_ub = (sup >= ekkinf);
00068
00069
00070
00071
00072 rdmin[i] = (no_lb && !no_ub) ? 0.0 : -PRESOLVE_INF;
00073
00074 rdmax[i] = (no_ub && !no_lb) ? 0.0 : PRESOLVE_INF;
00075 }
00076
00077
00078
00079 int j;
00080 for ( j = 0; j<ncols; j++) {
00081 bool no_ub = (cup[j] >= ekkinf);
00082 bool no_lb = (clo[j] <= -ekkinf);
00083
00084 if (hincol[j] == 1 &&
00085
00086
00087 (no_ub ^ no_lb)) {
00088 int row = hrow[mcstrt[j]];
00089 double coeff = colels[mcstrt[j]];
00090
00091 PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
00092
00093
00094
00095
00096 double dprice = maxmin * dcost[j] / coeff;
00097
00098
00099
00100
00101
00102 if ((coeff > 0.0) == no_ub) {
00103
00104
00105
00106
00107 if (rdmax[row] > dprice)
00108 rdmax[row] = dprice;
00109 } else {
00110 if (rdmin[row] < dprice)
00111 rdmin[row] = dprice;
00112 }
00113 }
00114 }
00115
00116 int *fixup_cols = new int[ncols];
00117
00118 int *fixdown_cols = new int[ncols];
00119
00120 #if PRESOLVE_TIGHTEN_DUALS
00121 double *djmin = new double[ncols];
00122 double *djmax = new double[ncols];
00123 #endif
00124 int nfixup_cols = 0;
00125 int nfixdown_cols = 0;
00126
00127 while (1) {
00128 int tightened = 0;
00129
00130 for (int j = 0; j<ncols; j++) {
00131 if (hincol[j] > 0) {
00132 CoinBigIndex kcs = mcstrt[j];
00133 CoinBigIndex kce = kcs + hincol[j];
00134
00135 int nflagu = 0;
00136 int nflagl = 0;
00137
00138 int nordu = 0;
00139 int nordl = 0;
00140 double ddjlo = maxmin * dcost[j];
00141 double ddjhi = ddjlo;
00142
00143 for (CoinBigIndex k = kcs; k < kce; k++) {
00144 int i = hrow[k];
00145 double coeff = colels[k];
00146
00147 if (coeff > 0.0) {
00148 if (rdmin[i] >= -ekkinf2) {
00149 ddjhi -= coeff * rdmin[i];
00150 nordu ++;
00151 } else {
00152 nflagu ++;
00153 }
00154 if (rdmax[i] <= ekkinf2) {
00155 ddjlo -= coeff * rdmax[i];
00156 nordl ++;
00157 } else {
00158 nflagl ++;
00159 }
00160 } else {
00161 if (rdmax[i] <= ekkinf2) {
00162 ddjhi -= coeff * rdmax[i];
00163 nordu ++;
00164 } else {
00165 nflagu ++;
00166 }
00167 if (rdmin[i] >= -ekkinf2) {
00168 ddjlo -= coeff * rdmin[i];
00169 nordl ++;
00170 } else {
00171 nflagl ++;
00172 }
00173 }
00174 }
00175
00176 if (cup[j]>ekkinf) {
00177
00178 if (nflagu==1&&ddjhi<-ztoldj) {
00179
00180 for (CoinBigIndex k = kcs; k < kce; k++) {
00181 int i = hrow[k];
00182 double coeff = colels[k];
00183
00184 if (coeff > 0.0&&rdmin[i] < -ekkinf2) {
00185
00186 if (ddjhi<rdmax[i]*coeff-ztoldj) {
00187 double newValue = ddjhi/coeff;
00188
00189 if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
00190 nflagl--;
00191 ddjlo -= coeff * newValue;
00192 } else if (rdmax[i] <= ekkinf2) {
00193 ddjlo -= coeff * (newValue-rdmax[i]);
00194 }
00195 rdmax[i] = newValue;
00196 tightened++;
00197 #if DEBUG_PRESOLVE
00198 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
00199 #endif
00200 }
00201 } else if (coeff < 0.0 && rdmax[i] > ekkinf2) {
00202
00203 if (ddjhi<rdmin[i]*coeff-ztoldj) {
00204 double newValue = ddjhi/coeff;
00205
00206 if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
00207 nflagl--;
00208 ddjlo -= coeff * newValue;
00209 } else if (rdmax[i] >= -ekkinf2) {
00210 ddjlo -= coeff * (newValue-rdmin[i]);
00211 }
00212 rdmin[i] = newValue;
00213 tightened++;
00214 #if DEBUG_PRESOLVE
00215 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
00216 #endif
00217 ddjlo = 0.0;
00218 }
00219 }
00220 }
00221 } else if (nflagl==0&&nordl==1&&ddjlo<-ztoldj) {
00222
00223 for (CoinBigIndex k = kcs; k < kce; k++) {
00224 int i = hrow[k];
00225 double coeff = colels[k];
00226
00227 if (coeff > 0.0) {
00228 rdmax[i] += ddjlo/coeff;
00229 ddjlo =0.0;
00230 tightened++;
00231 #if DEBUG_PRESOLVE
00232 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
00233 #endif
00234 } else if (coeff < 0.0 ) {
00235 rdmin[i] += ddjlo/coeff;
00236 ddjlo =0.0;
00237 tightened++;
00238 #if DEBUG_PRESOLVE
00239 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
00240 #endif
00241 }
00242 }
00243 }
00244 }
00245 #if 0
00246 if (clo[j]<-ekkinf) {
00247
00248 if (ddjlo>ztoldj&&nflagl==1) {
00249
00250 for (CoinBigIndex k = kcs; k < kce; k++) {
00251 int i = hrow[k];
00252 double coeff = colels[k];
00253
00254 if (coeff < 0.0&&rdmin[i] < -ekkinf2) {
00255
00256 if (ddjlo>rdmax[i]*coeff+ztoldj) {
00257 double newValue = ddjlo/coeff;
00258
00259 if (rdmax[i] > ekkinf2 && newValue <= ekkinf2) {
00260 nflagu--;
00261 ddjhi -= coeff * newValue;
00262 } else if (rdmax[i] <= ekkinf2) {
00263 ddjhi -= coeff * (newValue-rdmax[i]);
00264 }
00265 rdmax[i] = newValue;
00266 tightened++;
00267 #if DEBUG_PRESOLVE
00268 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
00269 #endif
00270 }
00271 } else if (coeff > 0.0 && rdmax[i] > ekkinf2) {
00272
00273 if (ddjlo>rdmin[i]*coeff+ztoldj) {
00274 double newValue = ddjlo/coeff;
00275
00276 if (rdmin[i] < -ekkinf2 && newValue >= -ekkinf2) {
00277 nflagu--;
00278 ddjhi -= coeff * newValue;
00279 } else if (rdmax[i] >= -ekkinf2) {
00280 ddjhi -= coeff * (newValue-rdmin[i]);
00281 }
00282 rdmin[i] = newValue;
00283 tightened++;
00284 #if DEBUG_PRESOLVE
00285 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
00286 #endif
00287 }
00288 }
00289 }
00290 } else if (nflagu==0&&nordu==1&&ddjhi>ztoldj) {
00291
00292 for (CoinBigIndex k = kcs; k < kce; k++) {
00293 int i = hrow[k];
00294 double coeff = colels[k];
00295
00296 if (coeff < 0.0) {
00297 rdmax[i] += ddjhi/coeff;
00298 ddjhi =0.0;
00299 tightened++;
00300 #if DEBUG_PRESOLVE
00301 printf("Col %d, row %d max pi now %g\n",j,i,rdmax[i]);
00302 #endif
00303 } else if (coeff > 0.0 ) {
00304 rdmin[i] += ddjhi/coeff;
00305 ddjhi =0.0;
00306 tightened++;
00307 #if DEBUG_PRESOLVE
00308 printf("Col %d, row %d min pi now %g\n",j,i,rdmin[i]);
00309 #endif
00310 }
00311 }
00312 }
00313 }
00314 #endif
00315
00316 #if PRESOLVE_TIGHTEN_DUALS
00317 djmin[j] = (nflagl ? -PRESOLVE_INF : ddjlo);
00318 djmax[j] = (nflagu ? PRESOLVE_INF : ddjhi);
00319 #endif
00320
00321 if (ddjlo > ztoldj && nflagl == 0&&!prob->colProhibited2(j)) {
00322
00323 if (clo[j] <= -ekkinf) {
00324 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
00325 prob->messages())
00326 <<j
00327 <<CoinMessageEol;
00328 prob->status_ |= 2;
00329 break;
00330 } else {
00331 fixdown_cols[nfixdown_cols++] = j;
00332 }
00333 } else if (ddjhi < -ztoldj && nflagu == 0&&!prob->colProhibited2(j)) {
00334
00335 if (cup[j] >= ekkinf) {
00336 prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
00337 prob->messages())
00338 <<j
00339 <<CoinMessageEol;
00340 prob->status_ |= 2;
00341 break;
00342 } else {
00343 fixup_cols[nfixup_cols++] = j;
00344 }
00345 }
00346 }
00347 }
00348
00349
00350
00351
00352 #if PRESOLVE_TIGHTEN_DUALS
00353 const double *rowels = prob->rowels_;
00354 const int *hcol = prob->hcol_;
00355 const CoinBigIndex *mrstrt = prob->mrstrt_;
00356 int *hinrow = prob->hinrow_;
00357
00358 for (int i = 0; i < nrows; i++) {
00359 bool no_ub = (rup[i] >= ekkinf);
00360 bool no_lb = (rlo[i] <= -ekkinf);
00361
00362 if ((no_ub ^ no_lb) == true) {
00363 CoinBigIndex krs = mrstrt[i];
00364 CoinBigIndex kre = krs + hinrow[i];
00365 double rmax = rdmax[i];
00366 double rmin = rdmin[i];
00367
00368
00369 for (CoinBigIndex k=krs; k<kre; k++) {
00370 double coeff = rowels[k];
00371 int icol = hcol[k];
00372 double djmax0 = djmax[icol];
00373 double djmin0 = djmin[icol];
00374
00375 if (no_ub) {
00376
00377 if (coeff > ZTOLDP && djmax0 <PRESOLVE_INF && cup[icol]>=ekkinf) {
00378 double bnd = djmax0 / coeff;
00379 if (rmax > bnd) {
00380 #if DEBUG_PRESOLVE
00381 printf("MAX TIGHT[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
00382 #endif
00383 rdmax[i] = rmax = bnd;
00384 tightened ++;;
00385 }
00386 } else if (coeff < -ZTOLDP && djmax0 <PRESOLVE_INF && cup[icol] >= ekkinf) {
00387 double bnd = djmax0 / coeff ;
00388 if (rmin < bnd) {
00389 #if DEBUG_PRESOLVE
00390 printf("MIN TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
00391 #endif
00392 rdmin[i] = rmin = bnd;
00393 tightened ++;;
00394 }
00395 }
00396 } else {
00397
00398 if (coeff > ZTOLDP && djmin0 > -PRESOLVE_INF && clo[icol]<=-ekkinf) {
00399 double bnd = djmin0 / coeff ;
00400 if (rmin < bnd) {
00401 #if DEBUG_PRESOLVE
00402 printf("MIN1 TIGHT[%d,%d]: %g --> %g\n", i, hrow[k], rdmin[i], bnd);
00403 #endif
00404 rdmin[i] = rmin = bnd;
00405 tightened ++;;
00406 }
00407 } else if (coeff < -ZTOLDP && djmin0 > -PRESOLVE_INF && clo[icol] <= -ekkinf) {
00408 double bnd = djmin0 / coeff ;
00409 if (rmax > bnd) {
00410 #if DEBUG_PRESOLVE
00411 printf("MAX TIGHT1[%d,%d]: %g --> %g\n", i,hrow[k], rdmax[i], bnd);
00412 #endif
00413 rdmax[i] = rmax = bnd;
00414 tightened ++;;
00415 }
00416 }
00417 }
00418 }
00419 }
00420 }
00421 #endif
00422
00423 if (!tightened||nfixdown_cols||nfixup_cols)
00424 break;
00425 #if PRESOLVE_TIGHTEN_DUALS
00426 else
00427 printf("DUAL TIGHTENED! %d\n", tightened);
00428 #endif
00429 }
00430
00431 if (nfixup_cols) {
00432 #if DEBUG_PRESOLVE
00433 printf("NDUAL: %d\n", nfixup_cols);
00434 #endif
00435 next = make_fixed_action::presolve(prob, fixup_cols, nfixup_cols, false, next);
00436 }
00437
00438 if (nfixdown_cols) {
00439 #if DEBUG_PRESOLVE
00440 printf("NDUAL: %d\n", nfixdown_cols);
00441 #endif
00442 next = make_fixed_action::presolve(prob, fixdown_cols, nfixdown_cols, true, next);
00443 }
00444
00445
00446
00447 int * canFix = (int *) rdmin;
00448 for ( i = 0; i < nrows; i++) {
00449 bool no_lb = (rlo[i] <= -ekkinf);
00450 bool no_ub = (rup[i] >= ekkinf);
00451 canFix[i]=0;
00452 if (no_ub && !no_lb ) {
00453 if ( rdmin[i]>0.0)
00454 canFix[i]=-1;
00455 else
00456 canFix[i]=-2;
00457 } else if (no_lb && !no_ub ) {
00458 if (rdmax[i]<0.0)
00459 canFix[i]=1;
00460 else
00461 canFix[i]=2;
00462 }
00463 }
00464 for (j = 0; j<ncols; j++) {
00465 if (hincol[j]<=1)
00466 continue;
00467 CoinBigIndex kcs = mcstrt[j];
00468 CoinBigIndex kce = kcs + hincol[j];
00469 int bindingUp=-1;
00470 int bindingDown=-1;
00471 if (cup[j]<ekkinf)
00472 bindingUp=-2;
00473 if (clo[j]>-ekkinf)
00474 bindingDown=-2;
00475 for (CoinBigIndex k = kcs; k < kce; k++) {
00476 int i = hrow[k];
00477 if (abs(canFix[i])!=2) {
00478 bindingUp=-2;
00479 bindingDown=-2;
00480 break;
00481 }
00482 double coeff = colels[k];
00483 if (coeff>0.0) {
00484 if (canFix[i]==2) {
00485
00486 if (bindingUp==-1)
00487 bindingUp = i;
00488 else
00489 bindingUp = -2;
00490 } else {
00491
00492 if (bindingDown==-1)
00493 bindingDown = i;
00494 else
00495 bindingDown = -2;
00496 }
00497 } else {
00498 if (canFix[i]==2) {
00499
00500 if (bindingDown==-1)
00501 bindingDown = i;
00502 else
00503 bindingDown = -2;
00504 } else {
00505
00506 if (bindingUp==-1)
00507 bindingUp = i;
00508 else
00509 bindingUp = -2;
00510 }
00511 }
00512 }
00513 double cost = maxmin * dcost[j];
00514 if (bindingUp>-2&&cost<=0.0) {
00515
00516 if (bindingUp>=0) {
00517 canFix[bindingUp] /= 2;
00518
00519 } else {
00520
00521 }
00522 } else if (bindingDown>-2 &&cost>=0.0) {
00523
00524 if (bindingDown>=0) {
00525 canFix[bindingDown] /= 2;
00526
00527 } else {
00528
00529 }
00530 }
00531 }
00532 for ( i = 0; i < nrows; i++) {
00533 if (canFix[i]==1) {
00534 rlo[i]=rup[i];
00535 prob->addRow(i);
00536 } else if (canFix[i]==-1) {
00537 rup[i]=rlo[i];
00538 prob->addRow(i);
00539 }
00540 }
00541
00542 delete[]rdmin;
00543 delete[]rdmax;
00544
00545 delete[]fixup_cols;
00546 delete[]fixdown_cols;
00547
00548 #if PRESOLVE_TIGHTEN_DUALS
00549 delete[]djmin;
00550 delete[]djmax;
00551 #endif
00552
00553 return (next);
00554 }