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