00001
00002
00003 #include <stdio.h>
00004 #include <math.h>
00005
00006 #include "CoinPresolveMatrix.hpp"
00007 #include "CoinPresolveFixed.hpp"
00008 #include "CoinPresolveDupcol.hpp"
00009 #include "CoinSort.hpp"
00010 #include "CoinPresolveUseless.hpp"
00011 #include "CoinMessage.hpp"
00012
00013 #define DSEED2 2147483647.0
00014
00015 const char *dupcol_action::name() const
00016 {
00017 return ("dupcol_action");
00018 }
00019
00020 const char *duprow_action::name() const
00021 {
00022 return ("duprow_action");
00023 }
00024
00025 void init_random_vec(double *work, int n)
00026 {
00027 double deseed = 12345678.0;
00028
00029 for (int i = 0; i < n; ++i) {
00030 deseed *= 16807.;
00031 int jseed = (int) (deseed / DSEED2);
00032 deseed -= (double) jseed * DSEED2;
00033 double random = deseed / DSEED2;
00034
00035 work[i]=random;
00036 }
00037 }
00038
00039 int compute_sums(const int *len, const CoinBigIndex *starts,
00040 int *index, double *elems,
00041 const double *work,
00042 double *sums, int *sorted, int n)
00043 {
00044 int nlook=0;
00045 for (int i = 0; i < n; ++i)
00046 if (len[i] > 0 ) {
00047 CoinBigIndex kcs = starts[i];
00048 CoinBigIndex kce = kcs + len[i];
00049
00050 double value=0.0;
00051
00052
00053 CoinSort_2(index+kcs,index+kcs+len[i],elems+kcs);
00054
00055
00056 for (CoinBigIndex k=kcs;k<kce;k++) {
00057 int irow=index[k];
00058 value += work[irow]*elems[k];
00059 }
00060 sums[nlook]=value;
00061 sorted[nlook]=i;
00062 ++nlook;
00063 }
00064 return (nlook);
00065 }
00066
00067 dupcol_action::dupcol_action(int nactions,
00068 const action *actions,
00069 const CoinPresolveAction *next) :
00070 CoinPresolveAction(next),
00071 nactions_(nactions), actions_(actions)
00072 {
00073 }
00074
00075
00076
00077
00078 const CoinPresolveAction *dupcol_action::presolve(CoinPresolveMatrix *prob,
00079 const CoinPresolveAction *next)
00080 {
00081 double *colels = prob->colels_;
00082 int *hrow = prob->hrow_;
00083 CoinBigIndex *mcstrt = prob->mcstrt_;
00084 int *hincol = prob->hincol_;
00085 int ncols = prob->ncols_;
00086
00087 double *clo = prob->clo_;
00088 double *cup = prob->cup_;
00089 double * sol = prob->sol_;
00090
00091 double *rowels = prob->rowels_;
00092 int *hcol = prob->hcol_;
00093 const CoinBigIndex *mrstrt = prob->mrstrt_;
00094 int *hinrow = prob->hinrow_;
00095 int nrows = prob->nrows_;
00096
00097
00098
00099
00100 double *dcost = prob->cost_;
00101
00102 const char *integerType = prob->integerType_;
00103
00104 double maxmin = prob->maxmin_;
00105
00106 action *actions = new action [ncols];
00107 int nactions = 0;
00108
00109 int *fixed_down = new int[ncols];
00110 int nfixed_down = 0;
00111 int *fixed_up = new int[ncols];
00112 int nfixed_up = 0;
00113
00114 double * workcol = new double[ncols];
00115 double * workrow = new double[nrows];
00116 int * sort = new int[ncols];
00117
00118
00119 init_random_vec(workrow, nrows);
00120
00121
00122
00123
00124
00125
00126 int nlook = compute_sums(hincol, mcstrt, hrow, colels, workrow, workcol, sort, ncols);
00127
00128 #if 1
00129
00130 if (maxmin < 0.0)
00131 return (next);
00132 #endif
00133
00134 CoinSort_2(workcol,workcol+nlook,sort);
00135
00136
00137 #if 0
00138
00139
00140
00141
00142 {
00143 double dval = workcol[0];
00144 int first = 0;
00145 for (int jj = 1; jj < nlook; jj++) {
00146 while (workcol[jj]==dval)
00147 jj++;
00148
00149 if (first + 1 < jj) {
00150 double buf[jj - first];
00151 for (int i=first; i<jj; ++i)
00152 buf[i-first] = dcost[sort[i]]*maxmin;
00153
00154 CoinSort_2(buf,buf+jj-first,sort+first);
00155
00156 }
00157 }
00158 }
00159 #endif
00160
00161
00162
00163
00164 for (int jj = 1; jj < nlook; jj++) {
00165 if (workcol[jj]==workcol[jj-1]) {
00166 int ithis=sort[jj];
00167 int ilast=sort[jj-1];
00168 CoinBigIndex kcs = mcstrt[ithis];
00169 CoinBigIndex kce = kcs + hincol[ithis];
00170
00171 if (hincol[ithis] == hincol[ilast]&&!prob->colProhibited2(ithis)) {
00172 int ishift = mcstrt[ilast] - kcs;
00173 CoinBigIndex k;
00174 for (k=kcs;k<kce;k++) {
00175 if (hrow[k] != hrow[k+ishift] ||
00176 colels[k] != colels[k+ishift]) {
00177 break;
00178 }
00179 }
00180 if (k == kce) {
00181
00182
00183
00184 double clo1=clo[ilast];
00185 double cup1=cup[ilast];
00186 double clo2=clo[ithis];
00187 double cup2=cup[ithis];
00188 double dcost1=dcost[ilast]*maxmin;
00189 double dcost2=dcost[ithis]*maxmin;
00190 double newSolution = sol[ilast]+sol[ithis];
00191
00192 if (clo1==cup1||clo2==cup2)
00193 abort();
00194
00195 if (
00196 integerType[ilast] || integerType[ithis]) {
00197 continue;
00198 } else if (dcost1==dcost2) {
00199
00200
00201
00202 double l_j = clo[ithis];
00203 double u_j = cup[ithis];
00204 double l_k = clo[ilast];
00205 double u_k = cup[ilast];
00206
00207 if (! (l_j + u_k <= l_k + u_j)) {
00208 swap(ilast, ithis);
00209 swap(clo1, clo2);
00210 swap(cup1, cup2);
00211
00212 }
00213
00214
00215
00216
00217 clo1 += clo2;
00218 cup1 += cup2;
00219 if (clo1<-1.0e20) {
00220 clo1=-1.0e31;
00221 }
00222 if (cup1>1.0e20) {
00223 cup1=1.0e31;
00224 }
00225
00226
00227
00228 {
00229 action *s = &actions[nactions];
00230 nactions++;
00231
00232 s->thislo = clo[ithis];
00233 s->thisup = cup[ithis];
00234 s->lastlo = clo[ilast];
00235 s->lastup = cup[ilast];
00236 s->ithis = ithis;
00237 s->ilast = ilast;
00238
00239 s->nincol = hincol[ithis];
00240 s->colels = copyOfArray(&colels[mcstrt[ithis]], hincol[ithis]);
00241 s->colrows= copyOfArray(&hrow[mcstrt[ithis]], hincol[ithis]);
00242 }
00243
00244
00245 clo[ilast]=clo1;
00246 cup[ilast]=cup1;
00247 sol[ilast] = newSolution;
00248 if (prob->getColumnStatus(ilast)==CoinPrePostsolveMatrix::basic||
00249 prob->getColumnStatus(ithis)==CoinPrePostsolveMatrix::basic)
00250 prob->setColumnStatus(ilast,CoinPrePostsolveMatrix::basic);
00251
00252
00253 dcost[ithis] = 0.0;
00254 sol[ithis]=clo2;
00255 {
00256 CoinBigIndex kcs = mcstrt[ithis];
00257 CoinBigIndex kce = kcs + hincol[ithis];
00258 for (CoinBigIndex k=kcs; k<kce; ++k)
00259
00260 presolve_delete_from_row(hrow[k], ithis, mrstrt, hinrow, hcol, rowels);
00261 }
00262 hincol[ithis] = 0;
00263
00264
00265 sort[jj] = ilast;
00266
00267
00268
00269 } else {
00270
00271
00272
00273
00274 if (clo1<-1.0e20 && clo2<-1.0e20) {
00275
00276
00277 #ifdef PRINT_DEBUG
00278 printf("** Odd columns %d and %d\n", ilast,ithis);
00279 #endif
00280 continue;
00281 } else if (clo1<-1.0e20) {
00282
00283 swap(ilast, ithis);
00284
00285 swap(clo1, clo2);
00286 swap(cup1, cup2);
00287 swap(dcost1, dcost2);
00288 }
00289
00290
00291 if (cup1>1.0e20) {
00292
00293 if (clo2<-1.0e20) {
00294 if (dcost1<dcost2) {
00295 #ifdef PRINT_DEBUG
00296 printf("** Problem seems unbounded %d and %d\n", ilast,ithis);
00297 #endif
00298 break;
00299 } else {
00300 continue;
00301 }
00302 } else if (cup2>1.0e20) {
00303
00304
00305
00306 if (dcost1>dcost2) {
00307 swap(ilast, ithis);
00308 swap(clo1, clo2);
00309 swap(cup1, cup2);
00310 swap(dcost1, dcost2);
00311 }
00312 sol[ilast] = newSolution;
00313 if (prob->getColumnStatus(ilast)==CoinPrePostsolveMatrix::basic||
00314 prob->getColumnStatus(ithis)==CoinPrePostsolveMatrix::basic)
00315 prob->setColumnStatus(ilast,CoinPrePostsolveMatrix::basic);
00316
00317 sol[ithis]=clo2;
00318 fixed_down[nfixed_down++] = ithis ;
00319 sort[jj] = ilast;
00320 } else {
00321
00322 if (dcost2>dcost1) {
00323
00324
00325 fixed_down[nfixed_down++] = ithis;
00326 sort[jj] = ilast;
00327 sol[ilast] = newSolution;
00328 if (prob->getColumnStatus(ilast)==CoinPrePostsolveMatrix::basic||
00329 prob->getColumnStatus(ithis)==CoinPrePostsolveMatrix::basic)
00330 prob->setColumnStatus(ilast,CoinPrePostsolveMatrix::basic);
00331
00332 sol[ithis]=clo2;
00333 } else {
00334
00335 continue;
00336 }
00337 }
00338 } else {
00339
00340 if (clo2<-1.0e20) {
00341
00342 if (dcost2>dcost1) {
00343
00344
00345 fixed_up[nfixed_up++] = ilast;
00346 sort[jj] = ithis;
00347 sol[ithis] = newSolution;
00348 if (prob->getColumnStatus(ithis)==CoinPrePostsolveMatrix::basic||
00349 prob->getColumnStatus(ilast)==CoinPrePostsolveMatrix::basic)
00350 prob->setColumnStatus(ithis,CoinPrePostsolveMatrix::basic);
00351
00352 sol[ilast]=clo1;
00353 } else {
00354
00355 continue;
00356 }
00357 } else if (cup2>1.0e20) {
00358 if (dcost2<dcost1) {
00359
00360
00361
00362 fixed_down[nfixed_down++] = ilast;
00363 sort[jj] = ithis;
00364 sol[ithis] = newSolution;
00365 if (prob->getColumnStatus(ithis)==CoinPrePostsolveMatrix::basic||
00366 prob->getColumnStatus(ilast)==CoinPrePostsolveMatrix::basic)
00367 prob->setColumnStatus(ithis,CoinPrePostsolveMatrix::basic);
00368
00369 sol[ilast]=clo1;
00370 } else {
00371
00372 continue;
00373 }
00374 } else {
00375
00376 continue;
00377 }
00378 }
00379 }
00380 }
00381 }
00382 }
00383 }
00384
00385 delete[]workrow;
00386 delete[]workcol;
00387 delete[]sort;
00388
00389
00390 if (nactions) {
00391 #if PRESOLVE_SUMMARY
00392 printf("DUPLICATE COLS: %d\n", nactions);
00393 #endif
00394 next = new dupcol_action(nactions, copyOfArray(actions,nactions), next);
00395 }
00396 deleteAction(actions,action*);
00397 if (nfixed_down)
00398 next = make_fixed_action::presolve(prob, fixed_down, nfixed_down,
00399 true,
00400 next);
00401
00402 if (nfixed_up)
00403 next = make_fixed_action::presolve(prob, fixed_up, nfixed_up,
00404 false,
00405 next);
00406
00407 delete[]fixed_down;
00408 delete[]fixed_up;
00409
00410 return (next);
00411 }
00412
00413 void create_col(int col, int n, int *rows, double *els,
00414 CoinBigIndex *mcstrt, double *colels, int *hrow, int *link,
00415 CoinBigIndex *free_listp)
00416 {
00417 CoinBigIndex free_list = *free_listp;
00418 int xstart = NO_LINK;
00419 for (int i=0; i<n; ++i) {
00420 CoinBigIndex k = free_list;
00421 free_list = link[free_list];
00422
00423 check_free_list(free_list);
00424
00425 hrow[k] = rows[i];
00426 colels[k] = els[i];
00427 link[k] = xstart;
00428 xstart = k;
00429 }
00430 mcstrt[col] = xstart;
00431 *free_listp = free_list;
00432 }
00433
00434
00435 void dupcol_action::postsolve(CoinPostsolveMatrix *prob) const
00436 {
00437 const action *const actions = actions_;
00438 const int nactions = nactions_;
00439
00440 double *clo = prob->clo_;
00441 double *cup = prob->cup_;
00442
00443 double *sol = prob->sol_;
00444 double *dcost = prob->cost_;
00445
00446 double *colels = prob->colels_;
00447 int *hrow = prob->hrow_;
00448 CoinBigIndex *mcstrt = prob->mcstrt_;
00449 int *hincol = prob->hincol_;
00450 int *link = prob->link_;
00451
00452 double *rcosts = prob->rcosts_;
00453
00454 CoinBigIndex free_list = prob->free_list_;
00455
00456 for (const action *f = &actions[nactions-1]; actions<=f; f--) {
00457 int icol = f->ithis;
00458 int icol2 = f->ilast;
00459
00460 dcost[icol] = dcost[icol2];
00461 clo[icol] = f->thislo;
00462 cup[icol] = f->thisup;
00463 clo[icol2] = f->lastlo;
00464 cup[icol2] = f->lastup;
00465
00466 create_col(icol, f->nincol, f->colrows, f->colels,
00467 mcstrt, colels, hrow, link, &free_list);
00468 hincol[icol] = hincol[icol2];
00469
00470 int nfinite = ((fabs(f->thislo) < PRESOLVE_INF) +
00471 (fabs(f->thisup) < PRESOLVE_INF) +
00472 (fabs(f->lastlo) < PRESOLVE_INF) +
00473 (fabs(f->lastup) < PRESOLVE_INF));
00474
00475 if (nfinite > 1) {
00476
00477 double u_j = f->thisup;
00478 double l_k = f->lastlo;
00479 double u_k = f->lastup;
00480 double x_k_sol = sol[icol2];
00481
00482 prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
00483 if (x_k_sol <= l_k + u_j) {
00484 sol[icol2] = l_k;
00485 sol[icol] = x_k_sol - l_k;
00486 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
00487 } else {
00488 sol[icol2] = u_k;
00489 sol[icol] = x_k_sol - u_k;
00490 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
00491 }
00492 } else if (nfinite == 1) {
00493 double x_k_sol = sol[icol2];
00494
00495 if (fabs(f->thislo) < PRESOLVE_INF) {
00496 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
00497 sol[icol] = f->thislo;
00498 sol[icol2] = x_k_sol - sol[icol];
00499 } else if (fabs(f->thisup) < PRESOLVE_INF) {
00500 prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
00501 sol[icol] = f->thisup;
00502 sol[icol2] = x_k_sol - sol[icol];
00503 } else if (fabs(f->lastlo) < PRESOLVE_INF) {
00504 prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
00505 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
00506 sol[icol2] = f->lastlo;
00507 sol[icol] = x_k_sol - sol[icol2];
00508 } else {
00509
00510 prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
00511 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atUpperBound);
00512 sol[icol2] = f->lastup;
00513 sol[icol] = x_k_sol - sol[icol2];
00514 }
00515 } else {
00516
00517 sol[icol] = 0.0;
00518 prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::isFree);
00519 }
00520
00521
00522
00523 rcosts[icol] = rcosts[icol2];
00524
00525 #ifdef DEBUG_PRESOLVE
00526 const double ztolzb = prob->ztolzb_;
00527 if (! (clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
00528 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol], sol[icol], cup[icol]);
00529 if (! (clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
00530 printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
00531 #endif
00532 }
00533 prob->free_list_ = free_list;
00534 }
00535
00536
00537
00538
00539
00540
00541 const CoinPresolveAction *duprow_action::presolve(CoinPresolveMatrix *prob,
00542 const CoinPresolveAction *next)
00543 {
00544
00545
00546 int ncols = prob->ncols_;
00547
00548 double *rowels = prob->rowels_;
00549 int *hcol = prob->hcol_;
00550 const CoinBigIndex *mrstrt = prob->mrstrt_;
00551 int *hinrow = prob->hinrow_;
00552 int nrows = prob->nrows_;
00553
00554 double *rlo = prob->rlo_;
00555 double *rup = prob->rup_;
00556
00557 int nuseless_rows = 0;
00558
00559 double * workcol = new double[ncols];
00560 double * workrow = new double[nrows];
00561 int * sort = new int[nrows];
00562
00563 init_random_vec(workcol, ncols);
00564
00565 int nlook = compute_sums(hinrow, mrstrt, hcol, rowels, workcol, workrow, sort, nrows);
00566 CoinSort_2(workrow,workrow+nlook,sort);
00567
00568 double dval = workrow[0];
00569 for (int jj = 1; jj < nlook; jj++) {
00570 if (workrow[jj]==dval) {
00571 int ithis=sort[jj];
00572 int ilast=sort[jj-1];
00573 CoinBigIndex krs = mrstrt[ithis];
00574 CoinBigIndex kre = krs + hinrow[ithis];
00575 if (hinrow[ithis] == hinrow[ilast]) {
00576 int ishift = mrstrt[ilast] - krs;
00577 CoinBigIndex k;
00578 for (k=krs;k<kre;k++) {
00579 if (hcol[k] != hcol[k+ishift] ||
00580 rowels[k] != rowels[k+ishift]) {
00581 break;
00582 }
00583 }
00584 if (k == kre) {
00585
00586 double rlo1=rlo[ilast];
00587 double rup1=rup[ilast];
00588 double rlo2=rlo[ithis];
00589 double rup2=rup[ithis];
00590
00591 int idelete=-1;
00592 if (rlo1<=rlo2) {
00593 if (rup2<=rup1) {
00594
00595 idelete=ilast;
00596 } else if (fabs(rlo1-rlo2)<1.0e-12) {
00597
00598 idelete=ithis;
00599
00600 sort[jj-1]=ithis;
00601 sort[jj]=ilast;
00602 } else {
00603
00604 #ifdef PRINT_DEBUG
00605 printf("overlapping duplicate row %g %g, %g %g\n",
00606 rlo1,rup1,rlo2,rup2);
00607 #endif
00608 if (rup1<rlo2) {
00609
00610 prob->status_|= 1;
00611
00612 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00613 prob->messages())
00614 <<ithis
00615 <<rlo[ithis]
00616 <<rup[ithis]
00617 <<CoinMessageEol;
00618 break;
00619 }
00620 }
00621 } else {
00622
00623 if (rup1<=rup2) {
00624
00625 idelete=ithis;
00626
00627 sort[jj-1]=ithis;
00628 sort[jj]=ilast;
00629 } else {
00630
00631 #ifdef PRINT_DEBUG
00632 printf("overlapping duplicate row %g %g, %g %g\n",
00633 rlo1,rup1,rlo2,rup2);
00634 #endif
00635 if (rup2<rlo1) {
00636
00637 prob->status_|= 1;
00638
00639 prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
00640 prob->messages())
00641 <<ithis
00642 <<rlo[ithis]
00643 <<rup[ithis]
00644 <<CoinMessageEol;
00645 break;
00646 }
00647 }
00648 }
00649 if (idelete>=0)
00650 sort[nuseless_rows++]=idelete;
00651 }
00652 }
00653 }
00654 dval=workrow[jj];
00655 }
00656
00657 delete[]workrow;
00658 delete[]workcol;
00659
00660
00661 if (nuseless_rows) {
00662 #if PRESOLVE_SUMMARY
00663 printf("DUPLICATE ROWS: %d\n", nuseless_rows);
00664 #endif
00665 next = useless_constraint_action::presolve(prob,
00666 sort, nuseless_rows,
00667 next);
00668 }
00669 delete[]sort;
00670
00671 return (next);
00672 }
00673
00674 void duprow_action::postsolve(CoinPostsolveMatrix *prob) const
00675 {
00676 printf("STILL NO POSTSOLVE FOR DUPROW!\n");
00677 abort();
00678 }