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