00001
00002
00003
00004 #include <stdio.h>
00005 #include <math.h>
00006
00007 #include "CoinHelperFunctions.hpp"
00008 #include "CoinPresolveMatrix.hpp"
00009
00010 #include "CoinPresolveEmpty.hpp"
00011 #include "CoinPresolveZeros.hpp"
00012 #include "CoinPresolveFixed.hpp"
00013 #include "CoinPresolveDoubleton.hpp"
00014
00015 #include "CoinPresolvePsdebug.hpp"
00016 #include "CoinMessage.hpp"
00017
00018 static double * presolve_duparray(const double * element, const int * index,
00019 int length, int offset,int row)
00020 {
00021 int n;
00022 length--;
00023 if (sizeof(double)==2*sizeof(int))
00024 n = (3*length+1)>>1;
00025 else
00026 n = 2*length;
00027 double * dArray = new double [n];
00028 int * iArray = (int *) (dArray+length);
00029 length++;
00030 int i,j=0;
00031 index += offset;
00032 element += offset;
00033 for (i=0;i<length;i++) {
00034 int iRow = index[i];
00035 if (iRow!=row) {
00036 dArray[j]=element[i];
00037 iArray[j++]=index[i];
00038 }
00039 }
00040 return dArray;
00041 }
00042
00043
00044 void compact_rep(double *elems, int *indices, CoinBigIndex *starts, const int *lengths, int n,
00045 const presolvehlink *link)
00046 {
00047 #if PRESOLVE_SUMMARY
00048 printf("****COMPACTING****\n");
00049 #endif
00050
00051
00052 int i = n;
00053 while (link[i].pre != NO_LINK)
00054 i = link[i].pre;
00055
00056 int j = 0;
00057 for (; i != n; i = link[i].suc) {
00058 CoinBigIndex s = starts[i];
00059 CoinBigIndex e = starts[i] + lengths[i];
00060
00061
00062 starts[i] = j;
00063 for (CoinBigIndex k = s; k < e; k++) {
00064 elems[j] = elems[k];
00065 indices[j] = indices[k];
00066 j++;
00067 }
00068 }
00069 }
00070
00071
00072 static bool expand_col(CoinBigIndex *mcstrt,
00073 double *colels,
00074 int *hrow,
00075 int *hincol,
00076 presolvehlink *clink, int ncols,
00077
00078 int icolx)
00079 {
00080 CoinBigIndex kcsx = mcstrt[icolx];
00081 CoinBigIndex kcex = kcsx + hincol[icolx];
00082
00083 const int maxk = mcstrt[ncols];
00084
00085
00086 int nextcol = clink[icolx].suc;
00087
00088
00089 if (kcex + 1 < mcstrt[nextcol] || nextcol == ncols) {
00090 if (! (kcex + 1 < mcstrt[nextcol])) {
00091
00092 compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
00093
00094
00095 kcsx = mcstrt[icolx];
00096 kcex = kcsx + hincol[icolx];
00097
00098 if (! (kcex + 1 < mcstrt[nextcol])) {
00099 return (true);
00100 }
00101 }
00102 } else {
00103
00104
00105
00106
00107 int lastcol = clink[ncols].pre;
00108
00109
00110
00111 int newkcsx = mcstrt[lastcol] + hincol[lastcol];
00112
00113 if (newkcsx + hincol[icolx] + 1 >= maxk) {
00114 compact_rep(colels, hrow, mcstrt, hincol, ncols, clink);
00115
00116
00117 kcsx = mcstrt[icolx];
00118 kcex = kcsx + hincol[icolx];
00119
00120 newkcsx = mcstrt[lastcol] + hincol[lastcol];
00121
00122 if (newkcsx + hincol[icolx] + 1 >= maxk) {
00123 return (true);
00124 }
00125
00126 kcsx = mcstrt[icolx];
00127 kcex = mcstrt[icolx] + hincol[icolx];
00128 }
00129
00130
00131 memcpy((void*)&hrow[newkcsx], (void*)&hrow[kcsx], hincol[icolx] * sizeof(int));
00132 memcpy((void*)&colels[newkcsx], (void*)&colels[kcsx], hincol[icolx] * sizeof(double));
00133
00134
00135 PRESOLVE_REMOVE_LINK(clink, icolx);
00136 PRESOLVE_INSERT_LINK(clink, icolx, lastcol);
00137
00138
00139 mcstrt[icolx] = newkcsx;
00140 kcsx = newkcsx;
00141 kcex = newkcsx + hincol[icolx];
00142 }
00143
00144 return (false);
00145 }
00146
00147 #if 0
00148 static bool reject_doubleton(int *mcstrt,
00149 double *colels,
00150 int *hrow,
00151 int *hincol,
00152 double coeff_factor,
00153 double max_coeff_ratio,
00154 int row0, int icolx, int icoly)
00155 {
00156 CoinBigIndex kcs = mcstrt[icoly];
00157 CoinBigIndex kce = kcs + hincol[icoly];
00158 CoinBigIndex kcsx = mcstrt[icolx];
00159 CoinBigIndex kcex = kcsx + hincol[icolx];
00160
00161 for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
00162 int row = hrow[kcoly];
00163
00164 if (row != row0) {
00165
00166 CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
00167
00168 if (kcolx<kcex) {
00169
00170
00171
00172 double orig = fabs(colels[kcoly] * coeff_factor);
00173 double addin = fabs(colels[kcolx]);
00174
00175 if (max_coeff_ratio * min(orig,addin) < max(orig,addin)) {
00176 #if DEBUG_PRESOLVE
00177 printf("REJECTED %d %g %g\n", row0, orig, addin);
00178 #endif
00179 return (true);
00180 }
00181
00182
00183 double newval = fabs((colels[kcoly] * coeff_factor) + colels[kcolx]);
00184
00185 if (max_coeff_ratio * newval < orig) {
00186 #if DEBUG_PRESOLVE
00187 printf("REJECTED1 %d %g %g\n", row0,
00188 (colels[kcoly] * coeff_factor),
00189 colels[kcolx]);
00190 #endif
00191 return (true);
00192 }
00193
00194 }
00195 }
00196 }
00197 return (false);
00198 }
00199 #endif
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 static bool elim_doubleton(const char *msg,
00239 CoinBigIndex *mcstrt,
00240 double *rlo, double * acts, double *rup,
00241 double *colels,
00242 int *hrow, int *hcol,
00243 int *hinrow, int *hincol,
00244 presolvehlink *clink, int ncols,
00245 CoinBigIndex *mrstrt, double *rowels,
00246
00247 double coeff_factor,
00248 double bounds_factor,
00249 int row0, int icolx, int icoly)
00250 {
00251 CoinBigIndex kcs = mcstrt[icoly];
00252 CoinBigIndex kce = kcs + hincol[icoly];
00253 CoinBigIndex kcsx = mcstrt[icolx];
00254 CoinBigIndex kcex = kcsx + hincol[icolx];
00255
00256
00257
00258 #if DEBUG_PRESOLVE
00259 printf("%s %d x=%d y=%d cf=%g bf=%g nx=%d yrows=(", msg,
00260 row0, icolx, icoly, coeff_factor, bounds_factor, hincol[icolx]);
00261 #endif
00262 for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
00263 int row = hrow[kcoly];
00264
00265
00266 PRESOLVEASSERT(kcex == kcsx + hincol[icolx]);
00267
00268
00269 if (row != row0) {
00270
00271 CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
00272
00273 if (bounds_factor != 0.0) {
00274
00275 if (-PRESOLVE_INF < rlo[row])
00276 rlo[row] -= colels[kcoly] * bounds_factor;
00277
00278
00279 if (rup[row] < PRESOLVE_INF)
00280 rup[row] -= colels[kcoly] * bounds_factor;
00281
00282
00283 acts[row] -= colels[kcoly] * bounds_factor;
00284 }
00285
00286 #if DEBUG_PRESOLVE
00287 printf("%d%s ", row, (kcolx<kcex) ? "+" : "");
00288 #endif
00289
00290 if (kcolx<kcex) {
00291
00292
00293
00294
00295
00296
00297 colels[kcolx] += colels[kcoly] * coeff_factor;
00298
00299
00300
00301 CoinBigIndex k2 = presolve_find_row(icolx, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
00302 rowels[k2] = colels[kcolx];
00303
00304
00305 presolve_delete_from_row(row, icoly, mrstrt, hinrow, hcol, rowels);
00306 } else {
00307
00308
00309
00310
00311
00312 {
00313 CoinBigIndex k2 = presolve_find_row(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
00314 hcol[k2] = icolx;
00315 rowels[k2] = colels[kcoly] * coeff_factor;
00316 }
00317
00318 {
00319 bool no_mem = expand_col(mcstrt, colels, hrow, hincol, clink, ncols,
00320 icolx);
00321 if (no_mem)
00322 return (true);
00323
00324
00325 kcoly = mcstrt[icoly] + (kcoly - kcs);
00326 kcs = mcstrt[icoly];
00327 kce = mcstrt[icoly] + hincol[icoly];
00328
00329 kcolx = mcstrt[icolx] + (kcolx - kcs);
00330 kcsx = mcstrt[icolx];
00331 kcex = mcstrt[icolx] + hincol[icolx];
00332 }
00333
00334
00335
00336 hrow[kcex] = row;
00337 colels[kcex] = colels[kcoly] * coeff_factor;
00338 hincol[icolx]++, kcex++;
00339 }
00340 }
00341 }
00342
00343 #if DEBUG_PRESOLVE
00344 printf(")\n");
00345 #endif
00346
00347
00348 hincol[icoly] = 0;
00349
00350 return (false);
00351 }
00352
00353
00354 void update_other_rep_quick(int col,
00355 const int *mcstrt, const int *hrow, const double *colels,
00356 const int *hincol,
00357
00358 const int *mrstrt, int *hcol, double *rowels, int *hinrow)
00359 {
00360 CoinBigIndex kcs = mcstrt[col];
00361 CoinBigIndex kce = kcs + hincol[col];
00362
00363 for (CoinBigIndex k=kcs; k<kce; ++k) {
00364 int row = hrow[k];
00365 double coeff = colels[k];
00366
00367 CoinBigIndex krs = mrstrt[row];
00368 CoinBigIndex kre = krs + hinrow[row];
00369
00370
00371
00372 hcol[kre] = col;
00373 rowels[kre] = coeff;
00374 ++hinrow[row];
00375 }
00376 }
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 const CoinPresolveAction *doubleton_action::presolve(CoinPresolveMatrix *prob,
00389 const CoinPresolveAction *next)
00390 {
00391 double *colels = prob->colels_;
00392 int *hrow = prob->hrow_;
00393 CoinBigIndex *mcstrt = prob->mcstrt_;
00394 int *hincol = prob->hincol_;
00395 int ncols = prob->ncols_;
00396
00397 double *clo = prob->clo_;
00398 double *cup = prob->cup_;
00399
00400 double *rowels = prob->rowels_;
00401 int *hcol = prob->hcol_;
00402 CoinBigIndex *mrstrt = prob->mrstrt_;
00403 int *hinrow = prob->hinrow_;
00404 int nrows = prob->nrows_;
00405
00406 double *rlo = prob->rlo_;
00407 double *rup = prob->rup_;
00408
00409 presolvehlink *clink = prob->clink_;
00410 presolvehlink *rlink = prob->rlink_;
00411
00412 const char *integerType = prob->integerType_;
00413
00414 double *cost = prob->cost_;
00415
00416 int numberLook = prob->numberRowsToDo_;
00417 int iLook;
00418 int * look = prob->rowsToDo_;
00419 const double ztolzb = prob->ztolzb_;
00420
00421 action * actions = new action [nrows];
00422 int nactions = 0;
00423
00424 int *zeros = new int[ncols];
00425 int nzeros = 0;
00426
00427 int *fixed = new int[ncols];
00428 int nfixed = 0;
00429
00430
00431 unsigned char *rowstat = prob->rowstat_;
00432 double *acts = prob->acts_;
00433 double * sol = prob->sol_;
00434
00435
00436
00437 #if CHECK_CONSISTENCY
00438 presolve_links_ok(clink, mcstrt, hincol, ncols);
00439 #endif
00440
00441
00442 for (iLook=0;iLook<numberLook;iLook++) {
00443 int irow = look[iLook];
00444 if (hinrow[irow] == 2 &&
00445 fabs(rup[irow] - rlo[irow]) <= ZTOLDP) {
00446 double rhs = rlo[irow];
00447 CoinBigIndex krs = mrstrt[irow];
00448 CoinBigIndex kre = krs + hinrow[irow];
00449 int icolx, icoly;
00450 CoinBigIndex k;
00451
00452
00453 for (k=krs; k<kre; k++) {
00454 if (hincol[hcol[k]] > 0) {
00455 break;
00456 }
00457 }
00458 PRESOLVEASSERT(k<kre);
00459 if (fabs(rowels[k]) < ZTOLDP)
00460 continue;
00461 icolx = hcol[k];
00462
00463
00464 for (k++; k<kre; k++) {
00465 if (hincol[hcol[k]] > 0) {
00466 break;
00467 }
00468 }
00469 PRESOLVEASSERT(k<kre);
00470 if (fabs(rowels[k]) < ZTOLDP)
00471 continue;
00472 icoly = hcol[k];
00473
00474
00475 if (!(fabs(cup[icolx] - clo[icolx]) < ZTOLDP) &&
00476 !(fabs(cup[icoly] - clo[icoly]) < ZTOLDP)) {
00477 double coeffx, coeffy;
00478
00479 CoinBigIndex krowx = presolve_find_row(irow, mcstrt[icolx], mcstrt[icolx] + hincol[icolx], hrow);
00480 CoinBigIndex krowy = presolve_find_row(irow, mcstrt[icoly], mcstrt[icoly] + hincol[icoly], hrow);
00481
00482
00483
00484
00485 int integerStatus=0;
00486 if (integerType[icolx]) {
00487 if (integerType[icoly]) {
00488
00489 int good = 0;
00490 double rhs2 = rhs;
00491 double value;
00492 value=colels[krowx];
00493 if (value<0.0) {
00494 value = - value;
00495 rhs2 += 1;
00496 }
00497 if (cup[icolx]==1.0&&clo[icolx]==0.0&&fabs(value-1.0)<1.0e-7)
00498 good =1;
00499 value=colels[krowy];
00500 if (value<0.0) {
00501 value = - value;
00502 rhs2 += 1;
00503 }
00504 if (cup[icoly]==1.0&&clo[icoly]==0.0&&fabs(value-1.0)<1.0e-7)
00505 good |= 2;
00506 if (good==3&&fabs(rhs2-1.0)<1.0e-7)
00507 integerStatus = 3;
00508 else
00509 integerStatus=-1;
00510 } else {
00511 integerStatus = 1;
00512 }
00513 } else if (integerType[icoly]) {
00514 integerStatus = 2;
00515 }
00516 if (integerStatus<0)
00517 continue;
00518 if (integerStatus == 2) {
00519 swap(icoly,icolx);
00520 swap(krowy,krowx);
00521 }
00522
00523
00524
00525
00526
00527
00528
00529 if (!integerStatus) {
00530 if (fabs(colels[krowy]) < fabs(colels[krowx])) {
00531 swap(icoly,icolx);
00532 swap(krowy,krowx);
00533 }
00534 }
00535
00536 #if 0
00537
00538 if (integerType[icolx] &&
00539 clo[icoly] != -PRESOLVE_INF &&
00540 cup[icoly] != PRESOLVE_INF) {
00541 continue;
00542 }
00543 #endif
00544
00545 {
00546 CoinBigIndex kcs = mcstrt[icoly];
00547 CoinBigIndex kce = kcs + hincol[icoly];
00548 for (k=kcs; k<kce; k++) {
00549 if (hinrow[hrow[k]] == 1) {
00550 break;
00551 }
00552 }
00553
00554 if (k<kce)
00555 continue;
00556 }
00557
00558 coeffx = colels[krowx];
00559 coeffy = colels[krowy];
00560
00561
00562
00563 if (hincol[icolx] == 1 && hincol[icoly] == 1)
00564 continue;
00565
00566
00567
00568
00569
00570 #if 0
00571 if (fabs(coeffx) * max_coeff_factor <= fabs(coeffy))
00572 continue;
00573 #endif
00574
00575 #if 0
00576 if (only_zero_rhs && rhs != 0)
00577 continue;
00578
00579 if (reject_doubleton(mcstrt, colels, hrow, hincol,
00580 -coeffx / coeffy,
00581 max_coeff_ratio,
00582 irow, icolx, icoly))
00583 continue;
00584 #endif
00585
00586
00587 {
00588 action *s = &actions[nactions];
00589 nactions++;
00590
00591 s->row = irow;
00592 s->icolx = icolx;
00593
00594 s->clox = clo[icolx];
00595 s->cupx = cup[icolx];
00596 s->costx = cost[icolx];
00597
00598 s->icoly = icoly;
00599 s->costy = cost[icoly];
00600
00601 s->rlo = rlo[irow];
00602
00603 s->coeffx = coeffx;
00604
00605 s->coeffy = coeffy;
00606
00607 s->ncolx = hincol[icolx];
00608
00609 s->ncoly = hincol[icoly];
00610 if (s->ncoly<s->ncolx) {
00611
00612 s->colel = presolve_duparray(colels, hrow, hincol[icoly],
00613 mcstrt[icoly],irow);
00614 s->ncolx=0;
00615 } else {
00616 s->colel = presolve_duparray(colels, hrow, hincol[icolx],
00617 mcstrt[icolx],irow);
00618 s->ncoly=0;
00619 }
00620 }
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634 {
00635 double lo1 = -PRESOLVE_INF;
00636 double up1 = PRESOLVE_INF;
00637
00638
00639
00640 if (-PRESOLVE_INF < clo[icoly]) {
00641 if (coeffx * coeffy < 0)
00642 lo1 = (coeffy * clo[icoly] - rhs) / -coeffx;
00643 else
00644 up1 = (coeffy * clo[icoly] - rhs) / -coeffx;
00645 }
00646
00647 if (cup[icoly] < PRESOLVE_INF) {
00648 if (coeffx * coeffy < 0)
00649 up1 = (coeffy * cup[icoly] - rhs) / -coeffx;
00650 else
00651 lo1 = (coeffy * cup[icoly] - rhs) / -coeffx;
00652 }
00653
00654
00655
00656 cost[icolx] += cost[icoly] * (-coeffx / coeffy);
00657
00658 prob->change_bias(cost[icoly] * rhs / coeffy);
00659
00660 if (0 ) {
00661 abort();
00662
00663 #if 0
00664 lo1 = trunc(lo1);
00665 up1 = trunc(up1);
00666
00667
00668
00669
00670
00671 if (lo1 > clo[icolx]) {
00672 (clo[icolx] <= 0.0)
00673 clo[icolx] = ? ilo
00674
00675 clo[icolx] = ilo;
00676 cup[icolx] = iup;
00677 }
00678 #endif
00679 } else {
00680 double lo2 = max(clo[icolx], lo1);
00681 double up2 = min(cup[icolx], up1);
00682 if (lo2 > up2) {
00683 if (lo2 <= up2 + prob->feasibilityTolerance_) {
00684
00685 double nearest = floor(lo2+0.5);
00686 if (fabs(nearest-lo2)<2.0*prob->feasibilityTolerance_) {
00687 lo2 = nearest;
00688 up2 = nearest;
00689 } else {
00690 lo2 = up2;
00691 }
00692 } else {
00693 prob->status_ |= 1;
00694 prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
00695 prob->messages())
00696 <<icolx
00697 <<lo2
00698 <<up2
00699 <<CoinMessageEol;
00700 break;
00701 }
00702 }
00703 clo[icolx] = lo2;
00704 cup[icolx] = up2;
00705
00706 if (rowstat) {
00707
00708 int basisChoice=0;
00709 int numberBasic=0;
00710 if (prob->columnIsBasic(icolx))
00711 numberBasic++;
00712 if (prob->columnIsBasic(icoly))
00713 numberBasic++;
00714 if (prob->rowIsBasic(irow))
00715 numberBasic++;
00716 if (sol[icolx]<=lo2+ztolzb) {
00717 sol[icolx] =lo2;
00718 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::atLowerBound);
00719 } else if (sol[icolx]>=up2-ztolzb) {
00720 sol[icolx] =up2;
00721 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::atUpperBound);
00722 } else {
00723 basisChoice=1;
00724 }
00725 if (numberBasic>1)
00726 prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::basic);
00727 }
00728 if (lo2 == up2)
00729 fixed[nfixed++] = icolx;
00730 }
00731 }
00732
00733
00734 {
00735 prob->addCol(icolx);
00736 int i,kcs,kce;
00737 kcs = mcstrt[icoly];
00738 kce = kcs + hincol[icoly];
00739 for (i=kcs;i<kce;i++) {
00740 int row = hrow[i];
00741 prob->addRow(row);
00742 }
00743 kcs = mcstrt[icolx];
00744 kce = kcs + hincol[icolx];
00745 for (i=kcs;i<kce;i++) {
00746 int row = hrow[i];
00747 prob->addRow(row);
00748 }
00749 }
00750
00751
00752 bool no_mem = elim_doubleton("ELIMD",
00753 mcstrt, rlo, acts, rup, colels,
00754 hrow, hcol, hinrow, hincol,
00755 clink, ncols,
00756 mrstrt, rowels,
00757 -coeffx / coeffy,
00758 rhs / coeffy,
00759 irow, icolx, icoly);
00760 if (no_mem)
00761 throwCoinError("out of memory",
00762 "doubleton_action::presolve");
00763
00764
00765
00766 presolve_delete_from_row(icolx, irow, mcstrt, hincol, hrow, colels);
00767
00768
00769 hinrow[irow] = 0;
00770
00771
00772 PRESOLVE_REMOVE_LINK(rlink, irow);
00773
00774
00775 PRESOLVE_REMOVE_LINK(clink, icoly);
00776 cost[icoly] = 0.0;
00777
00778 rlo[irow] = 0.0;
00779 rup[irow] = 0.0;
00780
00781 zeros[nzeros++] = icolx;
00782
00783
00784
00785
00786
00787 }
00788
00789 #if 0
00790 presolve_links_ok(clink, mcstrt, ncols);
00791 presolve_links_ok(rlink, mrstrt, nrows);
00792 prob->consistent();
00793 #endif
00794 }
00795 }
00796
00797 if (nactions) {
00798 #if PRESOLVE_SUMMARY
00799 printf("NDOUBLETONS: %d\n", nactions);
00800 #endif
00801 action *actions1 = new action[nactions];
00802 CoinMemcpyN(actions, nactions, actions1);
00803
00804 next = new doubleton_action(nactions, actions1, next);
00805
00806 if (nzeros)
00807 next = drop_zero_coefficients_action::presolve(prob, zeros, nzeros, next);
00808 if (nfixed)
00809 next = remove_fixed_action::presolve(prob, fixed, nfixed, next);
00810 }
00811
00812 delete[]zeros;
00813 delete[]fixed;
00814 deleteAction(actions,action*);
00815
00816 return (next);
00817 }
00818
00819
00820 void doubleton_action::postsolve(CoinPostsolveMatrix *prob) const
00821 {
00822 const action *const actions = actions_;
00823 const int nactions = nactions_;
00824
00825 double *colels = prob->colels_;
00826 int *hrow = prob->hrow_;
00827 CoinBigIndex *mcstrt = prob->mcstrt_;
00828 int *hincol = prob->hincol_;
00829 int *link = prob->link_;
00830
00831 double *clo = prob->clo_;
00832 double *cup = prob->cup_;
00833
00834 double *rlo = prob->rlo_;
00835 double *rup = prob->rup_;
00836
00837 double *dcost = prob->cost_;
00838
00839 double *sol = prob->sol_;
00840 double *rcosts = prob->rcosts_;
00841
00842 double *acts = prob->acts_;
00843 double *rowduals = prob->rowduals_;
00844
00845 unsigned char *colstat = prob->colstat_;
00846 unsigned char *rowstat = prob->rowstat_;
00847
00848 const double maxmin = prob->maxmin_;
00849
00850 char *cdone = prob->cdone_;
00851 char *rdone = prob->rdone_;
00852
00853 CoinBigIndex free_list = prob->free_list_;
00854
00855 const double ztolzb = prob->ztolzb_;
00856 const double ztoldj = prob->ztoldj_;
00857
00858
00859 int nrows = prob->nrows_;
00860 int * index1 = new int[nrows];
00861 double * element1 = new double[nrows];
00862 memset(element1,0,nrows*sizeof(double));
00863
00864 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
00865 int irow = f->row;
00866 double lo0 = f->clox;
00867 double up0 = f->cupx;
00868
00869
00870 double coeffx = f->coeffx;
00871 double coeffy = f->coeffy;
00872 int jcolx = f->icolx;
00873 int jcoly = f->icoly;
00874
00875
00876 double rhs = f->rlo;
00877
00878
00879 PRESOLVEASSERT(cdone[jcolx] && rdone[irow]==DROP_ROW);
00880 PRESOLVEASSERT(cdone[jcoly]==DROP_COL);
00881
00882
00883 rlo[irow] = f->rlo;
00884 rup[irow] = f->rlo;
00885
00886 clo[jcolx] = lo0;
00887 cup[jcolx] = up0;
00888
00889 dcost[jcolx] = f->costx;
00890 dcost[jcoly] = f->costy;
00891
00892 #if DEBUG_PRESOLVE
00893
00894 if ((rhs < 0) == ((coeffx * sol[jcolx]) < 0) &&
00895 fabs(rhs - coeffx * sol[jcolx]) * 100 < rhs &&
00896 fabs(rhs - coeffx * sol[jcolx]) * 100 < (coeffx * sol[jcolx]))
00897 printf("DANGEROUS RHS??? %g %g %g\n",
00898 rhs, coeffx * sol[jcolx],
00899 (rhs - coeffx * sol[jcolx]));
00900 #endif
00901
00902 sol[jcoly] = (rhs - coeffx * sol[jcolx]) / coeffy;
00903
00904
00905 acts[irow] = rhs;
00906
00907
00908 if (rowstat)
00909 prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
00910
00911
00912
00913
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934 double djy = maxmin * dcost[jcoly];
00935 double djx = maxmin * dcost[jcolx];
00936 double bounds_factor = rhs/coeffy;
00937 if (f->ncoly) {
00938
00939 int ncoly=f->ncoly-1;
00940 double multiplier = coeffx/coeffy;
00941
00942 int * indy = (int *) (f->colel+ncoly);
00943 int ystart = NO_LINK;
00944 int nX=0;
00945 int i,iRow;
00946 for (i=0; i<ncoly; ++i) {
00947 int iRow = indy[i];
00948 double yValue = f->colel[i];
00949 CoinBigIndex k = free_list;
00950 free_list = link[free_list];
00951
00952 check_free_list(free_list);
00953
00954
00955
00956 if (-PRESOLVE_INF < rlo[iRow])
00957 rlo[iRow] += yValue * bounds_factor;
00958
00959
00960 if (rup[iRow] < PRESOLVE_INF)
00961 rup[iRow] += yValue * bounds_factor;
00962
00963 acts[iRow] += yValue * bounds_factor;
00964
00965 djy -= rowduals[iRow] * yValue;
00966
00967 hrow[k] = iRow;
00968 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
00969 colels[k] = yValue;
00970 link[k] = ystart;
00971 ystart = k;
00972 yValue *= multiplier;
00973 element1[iRow]=yValue;
00974 index1[nX++]=iRow;
00975 }
00976
00977 {
00978 double yValue = coeffy;
00979 CoinBigIndex k = free_list;
00980 free_list = link[free_list];
00981
00982 check_free_list(free_list);
00983
00984 hrow[k] = irow;
00985 colels[k] = yValue;
00986 link[k] = ystart;
00987 ystart = k;
00988 yValue *= multiplier;
00989 element1[irow]=yValue;
00990 index1[nX++]=irow;
00991 }
00992 mcstrt[jcoly] = ystart;
00993 hincol[jcoly] = f->ncoly;
00994
00995 CoinBigIndex k=mcstrt[jcolx];
00996 CoinBigIndex last = NO_LINK;
00997 int numberInColumn = hincol[jcolx];
00998 int numberToDo=numberInColumn;
00999 for (i=0; i<numberToDo; ++i) {
01000 iRow = hrow[k];
01001 assert (iRow>=0&&iRow<nrows);
01002 double value = colels[k]+element1[iRow];
01003 element1[iRow]=0.0;
01004 if (fabs(value)>=1.0e-15) {
01005 colels[k]=value;
01006 last=k;
01007 k = link[k];
01008 if (iRow != irow)
01009 djx -= rowduals[iRow] * value;
01010 } else {
01011 numberInColumn--;
01012
01013 int nextk = link[k];
01014 assert(free_list>=0);
01015 link[k]=free_list;
01016 free_list=k;
01017 assert (k>=0);
01018 k=nextk;
01019 if (last!=NO_LINK)
01020 link[last]=k;
01021 else
01022 mcstrt[jcolx]=k;
01023 }
01024 }
01025 for (i=0;i<nX;i++) {
01026 int iRow = index1[i];
01027 double xValue = element1[iRow];
01028 element1[iRow]=0.0;
01029 if (fabs(xValue)>=1.0e-15) {
01030 if (iRow != irow)
01031 djx -= rowduals[iRow] * xValue;
01032 numberInColumn++;
01033 CoinBigIndex k = free_list;
01034 free_list = link[free_list];
01035
01036 check_free_list(free_list);
01037
01038 hrow[k] = iRow;
01039 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
01040 colels[k] = xValue;
01041 link[last] = k;
01042 last = k;
01043 }
01044 }
01045 link[last]=NO_LINK;
01046 assert(numberInColumn);
01047 hincol[jcolx] = numberInColumn;
01048 } else {
01049
01050 int ncolx=f->ncolx-1;
01051 double multiplier = -coeffy/coeffx;
01052 int * indx = (int *) (f->colel+ncolx);
01053
01054
01055 CoinBigIndex k=mcstrt[jcolx];
01056 int nX=0;
01057 int i,iRow;
01058 for (i=0; i<hincol[jcolx]-1; ++i) {
01059 if (colels[k]) {
01060 iRow = hrow[k];
01061 index1[nX++]=iRow;
01062 element1[iRow]=multiplier*colels[k];
01063 }
01064 k = link[k];
01065 }
01066 iRow = hrow[k];
01067 index1[nX++]=iRow;
01068 element1[iRow]=multiplier*colels[k];
01069 multiplier = - multiplier;
01070 link[k] = free_list;
01071 free_list = mcstrt[jcolx];
01072 int xstart = NO_LINK;
01073 for (i=0; i<ncolx; ++i) {
01074 int iRow = indx[i];
01075 double xValue = f->colel[i];
01076
01077 CoinBigIndex k = free_list;
01078 free_list = link[free_list];
01079
01080 check_free_list(free_list);
01081
01082 hrow[k] = iRow;
01083 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
01084 colels[k] = xValue;
01085 link[k] = xstart;
01086 xstart = k;
01087 xValue *= multiplier;
01088 if (!element1[iRow]) {
01089 element1[iRow]=xValue;
01090 index1[nX++]=iRow;
01091 } else {
01092 element1[iRow]+=xValue;
01093 }
01094 }
01095
01096 {
01097 double xValue = coeffx;
01098 CoinBigIndex k = free_list;
01099 free_list = link[free_list];
01100
01101 check_free_list(free_list);
01102
01103 hrow[k] = irow;
01104 colels[k] = xValue;
01105 link[k] = xstart;
01106 xstart = k;
01107 xValue *= multiplier;
01108 if (!element1[irow]) {
01109 element1[irow]=xValue;
01110 index1[nX++]=irow;
01111 } else {
01112 element1[irow]+=xValue;
01113 }
01114 }
01115 mcstrt[jcolx] = xstart;
01116 hincol[jcolx] = f->ncolx;
01117 int ystart = NO_LINK;
01118 int n=0;
01119 for (i=0;i<nX;i++) {
01120 int iRow = index1[i];
01121 double yValue = element1[iRow];
01122 element1[iRow]=0.0;
01123 if (fabs(yValue)>=1.0e-12) {
01124 n++;
01125 CoinBigIndex k = free_list;
01126 free_list = link[free_list];
01127
01128 check_free_list(free_list);
01129
01130 hrow[k] = iRow;
01131 PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
01132 colels[k] = yValue;
01133 link[k] = ystart;
01134 ystart = k;
01135 }
01136 }
01137 mcstrt[jcoly] = ystart;
01138 assert(n);
01139 hincol[jcoly] = n;
01140
01141 k = mcstrt[jcoly];
01142 int ny = hincol[jcoly];
01143
01144 for (i=0; i<ny; ++i) {
01145 int row = hrow[k];
01146 double coeff = colels[k];
01147 k = link[k];
01148
01149 if (row != irow) {
01150
01151
01152
01153 if (-PRESOLVE_INF < rlo[row])
01154 rlo[row] += coeff * bounds_factor;
01155
01156
01157 if (rup[row] < PRESOLVE_INF)
01158 rup[row] += coeff * bounds_factor;
01159
01160 acts[row] += coeff * bounds_factor;
01161
01162 djy -= rowduals[row] * coeff;
01163 }
01164 }
01165 k = mcstrt[jcolx];
01166 int nx = hincol[jcolx];
01167
01168 for ( i=0; i<nx; ++i) {
01169 int row = hrow[k];
01170 double coeff = colels[k];
01171 k = link[k];
01172
01173 if (row != irow) {
01174 djx -= rowduals[row] * coeff;
01175 }
01176 }
01177 }
01178 assert (fabs(coeffx-f->coeffx)<1.0e-6&&fabs(coeffy-f->coeffy)<1.0e-6);
01179
01180
01181
01182
01183
01184
01185 if (colstat) {
01186 if (prob->columnIsBasic(jcolx) ||
01187 (fabs(lo0 - sol[jcolx]) < ztolzb && rcosts[jcolx] >= -ztoldj) ||
01188 (fabs(up0 - sol[jcolx]) < ztolzb && rcosts[jcolx] <= ztoldj)) {
01189
01190
01191 prob->setColumnStatus(jcoly,CoinPrePostsolveMatrix::basic);
01192
01193
01194 rowduals[irow] = djy / coeffy;
01195 rcosts[jcolx] = djx - rowduals[irow] * coeffx;
01196 #if 0
01197 if (prob->columnIsBasic(jcolx))
01198 assert (fabs(rcosts[jcolx])<1.0e-5);
01199 #endif
01200 rcosts[jcoly] = 0.0;
01201 } else {
01202 prob->setColumnStatus(jcolx,CoinPrePostsolveMatrix::basic);
01203 prob->setColumnStatusUsingValue(jcoly);
01204
01205
01206 rowduals[irow] = djx / coeffx;
01207 rcosts[jcoly] = djy - rowduals[irow] * coeffy;
01208 rcosts[jcolx] = 0.0;
01209 }
01210 } else {
01211
01212
01213
01214 rowduals[irow] = djy / coeffy;
01215 rcosts[jcoly] = 0.0;
01216 }
01217
01218
01219 #if DEBUG_PRESOLVE
01220 {
01221 CoinBigIndex k = mcstrt[jcolx];
01222 int nx = hincol[jcolx];
01223 double dj = maxmin * dcost[jcolx];
01224
01225 for (int i=0; i<nx; ++i) {
01226 int row = hrow[k];
01227 double coeff = colels[k];
01228 k = link[k];
01229
01230 dj -= rowduals[row] * coeff;
01231 }
01232 if (! (fabs(rcosts[jcolx] - dj) < 100*ZTOLDP))
01233 printf("BAD DOUBLE X DJ: %d %d %g %g\n",
01234 irow, jcolx, rcosts[jcolx], dj);
01235 rcosts[jcolx]=dj;
01236 }
01237 {
01238 CoinBigIndex k = mcstrt[jcoly];
01239 int ny = hincol[jcoly];
01240 double dj = maxmin * dcost[jcoly];
01241
01242 for (int i=0; i<ny; ++i) {
01243 int row = hrow[k];
01244 double coeff = colels[k];
01245 k = link[k];
01246
01247 dj -= rowduals[row] * coeff;
01248
01249
01250 }
01251 if (! (fabs(rcosts[jcoly] - dj) < 100*ZTOLDP))
01252 printf("BAD DOUBLE Y DJ: %d %d %g %g\n",
01253 irow, jcoly, rcosts[jcoly], dj);
01254 rcosts[jcoly]=dj;
01255
01256 }
01257 #endif
01258
01259 cdone[jcoly] = DOUBLETON;
01260 rdone[irow] = DOUBLETON;
01261 }
01262 delete [] index1;
01263 delete [] element1;
01264 prob->free_list_ = free_list;
01265 }
01266
01267
01268 doubleton_action::~doubleton_action()
01269 {
01270 for (int i=nactions_-1; i>=0; i--) {
01271 delete[]actions_[i].colel;
01272 }
01273 deleteAction(actions_,action*);
01274 }
01275
01276
01277
01278 static double *doubleton_mult;
01279 static int *doubleton_id;
01280 void check_doubletons(const CoinPresolveAction * paction)
01281 {
01282 const CoinPresolveAction * paction0 = paction;
01283
01284 if (paction) {
01285 check_doubletons(paction->next);
01286
01287 if (strcmp(paction0->name(), "doubleton_action") == 0) {
01288 const doubleton_action *daction = (const doubleton_action *)paction0;
01289 for (int i=daction->nactions_-1; i>=0; --i) {
01290 int icolx = daction->actions_[i].icolx;
01291 int icoly = daction->actions_[i].icoly;
01292 double coeffx = daction->actions_[i].coeffx;
01293 double coeffy = daction->actions_[i].coeffy;
01294
01295 doubleton_mult[icoly] = -coeffx/coeffy;
01296 doubleton_id[icoly] = icolx;
01297 }
01298 }
01299 }
01300 }
01301
01302 void check_doubletons1(const CoinPresolveAction * paction,
01303 int ncols)
01304 {
01305 #if DEBUG_PRESOLVE
01306 doubleton_mult = new double[ncols];
01307 doubleton_id = new int[ncols];
01308 for (int i=0; i<ncols; ++i)
01309 doubleton_id[i] = i;
01310 check_doubletons(paction);
01311 double minmult = 1.0;
01312 int minid = -1;
01313 for (int i=0; i<ncols; ++i) {
01314 double mult = 1.0;
01315 int j = i;
01316 if (doubleton_id[j] != j) {
01317 printf("MULTS (%d): ", j);
01318 while (doubleton_id[j] != j) {
01319 printf("%d %g, ", doubleton_id[j], doubleton_mult[j]);
01320 mult *= doubleton_mult[j];
01321 j = doubleton_id[j];
01322 }
01323 printf(" == %g\n", mult);
01324 if (minmult > fabs(mult)) {
01325 minmult = fabs(mult);
01326 minid = i;
01327 }
01328 }
01329 }
01330 if (minid != -1)
01331 printf("MIN MULT: %d %g\n", minid, minmult);
01332 #endif
01333 }