// Copyright (C) 2002, International Business Machines // Corporation and others. All Rights Reserved. #include #include #include "CoinPresolveMatrix.hpp" #include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW #include "CoinPresolveFixed.hpp" #include "CoinPresolveSubst.hpp" #include "CoinHelperFunctions.hpp" #include "CoinPresolveUseless.hpp" #include "CoinPresolveForcing.hpp" #include "CoinMessage.hpp" /*static*/ void implied_bounds(const double *els, const double *clo, const double *cup, const int *hcol, CoinBigIndex krs, CoinBigIndex kre, double *maxupp, double *maxdownp, int jcol, double rlo, double rup, double *iclb, double *icub) { if (rlo<=-PRESOLVE_INF&&rup>=PRESOLVE_INF) { *iclb = -PRESOLVE_INF; *icub = PRESOLVE_INF; return; } bool posinf = false; bool neginf = false; double maxup = 0.0; double maxdown = 0.0; int jcolk = -1; // compute sum of all bounds except for jcol CoinBigIndex kk; for (kk=krs; kk ZTOLDP); double ilb = (rlo - maxup) / coeff; bool finite_ilb = (-PRESOLVE_INF < rlo && !posinf && PRESOLVEFINITE(maxup)); double iub = (rup - maxdown) / coeff; bool finite_iub = ( rup < PRESOLVE_INF && !neginf && PRESOLVEFINITE(maxdown)); if (coeff > 0.0) { *iclb = (finite_ilb ? ilb : -PRESOLVE_INF); *icub = (finite_iub ? iub : PRESOLVE_INF); } else { *iclb = (finite_iub ? iub : -PRESOLVE_INF); *icub = (finite_ilb ? ilb : PRESOLVE_INF); } } if (coeff > 0.0) { if (PRESOLVE_INF <= ub) { posinf = true; if (neginf) break; // pointless } else maxup += ub * coeff; if (lb <= -PRESOLVE_INF) { neginf = true; if (posinf) break; // pointless } else maxdown += lb * coeff; } else { if (PRESOLVE_INF <= ub) { neginf = true; if (posinf) break; // pointless } else maxdown += ub * coeff; if (lb <= -PRESOLVE_INF) { posinf = true; if (neginf) break; // pointless } else maxup += lb * coeff; } } // If we broke from the loop, then the bounds are infinite. // However, since we put the column whose implied bounds we want // to know at the end, and it doesn't matter if its own bounds // are infinite, don't worry about breaking at the last iteration. if (kktuning_) { startTime = CoinCpuTime(); startEmptyRows = prob->countEmptyRows(); startEmptyColumns = prob->countEmptyCols(); } double *clo = prob->clo_; double *cup = prob->cup_; const double *rowels = prob->rowels_; const int *hcol = prob->hcol_; const CoinBigIndex *mrstrt = prob->mrstrt_; const int *hinrow = prob->hinrow_; const int nrows = prob->nrows_; const double *rlo = prob->rlo_; const double *rup = prob->rup_; // const char *integerType = prob->integerType_; const double tol = ZTOLDP; const double inftol = prob->feasibilityTolerance_; const int ncols = prob->ncols_; int *fixed_cols = new int[ncols]; int nfixed_cols = 0; action *actions = new action [nrows]; int nactions = 0; int *useless_rows = new int[nrows]; int nuseless_rows = 0; int numberLook = prob->numberRowsToDo_; int iLook; int * look = prob->rowsToDo_; bool fixInfeasibility = (prob->presolveOptions_&16384)!=0; /* Open a loop to scan the constraints of interest. There must be variables left in the row. */ for (iLook=0;iLook 0) { CoinBigIndex krs = mrstrt[irow]; CoinBigIndex kre = krs + hinrow[irow]; /* Calculate upper and lower bounds on the row activity based on upper and lower bounds on the variables. If these are finite and incompatible with the given row bounds, we have infeasibility. */ double maxup, maxdown; implied_row_bounds(rowels, clo, cup, hcol, krs, kre, &maxup, &maxdown); if (maxup < PRESOLVE_INF && maxup + inftol < rlo[irow]&&!fixInfeasibility) { /* there is an upper bound and it can't be reached */ prob->status_|= 1; prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, prob->messages()) <status_|= 1; prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, prob->messages()) <= PRESOLVE_INF || (maxup < PRESOLVE_INF && rup[irow] >= maxup))) { /* Original comment: I'm not sure that these transforms don't intefere with each other. We can get it next time. Well, I'll argue that bounds are never really loosened (at worst, they're transferred onto some other variable, or inferred to be unnecessary. Once useless, always useless. Leaving this hook in place allows for a sort of infinite loop where this routine keeps queuing the same constraints over and over. -- lh, 040901 -- if (some_col_was_fixed(hcol, krs, kre, clo, cup)) { prob->addRow(irow); continue; } */ // this constraint must always be satisfied - drop it useless_rows[nuseless_rows++] = irow; } else if ((maxup < PRESOLVE_INF && fabs(rlo[irow] - maxup) < tol) || (-PRESOLVE_INF < maxdown && fabs(rup[irow] - maxdown) < tol)) { // the lower bound can just be reached, or // the upper bound can just be reached; // called a "forcing constraint" in the paper (p. 226) const int lbound_tight = (maxup < PRESOLVE_INF && fabs(rlo[irow] - maxup) < tol); /* Original comment and rebuttal as above. if (some_col_was_fixed(hcol, krs, kre, clo, cup)) { // make sure on next time prob->addRow(irow); continue; } */ // out of space - this probably never happens (but this routine will // often put duplicates in the fixed column list) if (nfixed_cols + (kre-krs) >= ncols) break; double *bounds = new double[hinrow[irow]]; int *rowcols = new int[hinrow[irow]]; int lk = krs; // load fix-to-down in front int uk = kre; // load fix-to-up in back CoinBigIndex k; for ( k=krs; kaddCol(jcol); double coeff = rowels[k]; PRESOLVEASSERT(fabs(coeff) > ZTOLDP); // one of the two contributed to maxup - set the other to that if (lbound_tight == (coeff > 0.0)) { --uk; bounds[uk-krs] = clo[jcol]; rowcols[uk-krs] = jcol; clo[jcol] = cup[jcol]; } else { bounds[lk-krs] = cup[jcol]; rowcols[lk-krs] = jcol; ++lk; cup[jcol] = clo[jcol]; } fixed_cols[nfixed_cols++] = jcol; } PRESOLVEASSERT(uk == lk); action *f = &actions[nactions]; nactions++; f->row = irow; f->nlo = lk-krs; f->nup = kre-uk; f->rowcols = rowcols; f->bounds = bounds; } } } if (nactions) { #if PRESOLVE_SUMMARY printf("NFORCED: %d\n", nactions); #endif next = new forcing_constraint_action(nactions, CoinCopyOfArray(actions,nactions), next); } deleteAction(actions,action*); if (nuseless_rows) { next = useless_constraint_action::presolve(prob, useless_rows, nuseless_rows, next); } delete[]useless_rows; /* We need to remove duplicates here, or we get into trouble in remove_fixed_action::postsolve when we try to reinstate a column multiple times. */ if (nfixed_cols) { if (nfixed_cols > 1) { std::sort(fixed_cols,fixed_cols+nfixed_cols) ; int *end = std::unique(fixed_cols,fixed_cols+nfixed_cols) ; nfixed_cols = static_cast(end-fixed_cols) ; } next = remove_fixed_action::presolve(prob,fixed_cols,nfixed_cols,next) ; } delete[]fixed_cols ; if (prob->tuning_) { double thisTime=CoinCpuTime(); int droppedRows = prob->countEmptyRows() - startEmptyRows ; int droppedColumns = prob->countEmptyCols() - startEmptyColumns; printf("CoinPresolveForcing(32) - %d rows, %d columns dropped in time %g, total %g\n", droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_); } return (next); } void forcing_constraint_action::postsolve(CoinPostsolveMatrix *prob) const { const action *const actions = actions_; const int nactions = nactions_; const double *colels = prob->colels_; const int *hrow = prob->hrow_; const CoinBigIndex *mcstrt = prob->mcstrt_; const int *hincol = prob->hincol_; const int *link = prob->link_; // CoinBigIndex free_list = prob->free_list_; double *clo = prob->clo_; double *cup = prob->cup_; const double *sol = prob->sol_; double *rcosts = prob->rcosts_; //double *acts = prob->acts_; double *rowduals = prob->rowduals_; const double ztoldj = prob->ztoldj_; const double ztolzb = prob->ztolzb_; for (const action *f = &actions[nactions-1]; actions<=f; f--) { const int irow = f->row; const int nlo = f->nlo; const int nup = f->nup; const int ninrow = nlo + nup; const int *rowcols = f->rowcols; const double *bounds= f->bounds; int k; /* Original comment: When we restore bounds here, we need to allow for the possibility that the restored bound is infinite. This implies a check for viable status. Hmmm ... I'm going to argue that in fact we have no choice: the status of the variable must reflect the value it was fixed at, else we lose feasibility. We don't care what the other bound does. -- lh, 040903 -- */ for (k=0; ksetColumnStatus(jcol,CoinPrePostsolveMatrix::atLowerBound) ; /* PRESOLVEASSERT(prob->getColumnStatus(jcol)!=CoinPrePostsolveMatrix::basic); if (cup[jcol] >= PRESOLVE_INF) { CoinPrePostsolveMatrix::Status statj = prob->getColumnStatus(jcol) ; if (statj == CoinPrePostsolveMatrix::atUpperBound) { if (clo[jcol] > -PRESOLVE_INF) { statj = CoinPrePostsolveMatrix::atLowerBound ; } else { statj = CoinPrePostsolveMatrix::isFree ; } prob->setColumnStatus(jcol,statj) ; } } */ } for (k=nlo; ksetColumnStatus(jcol,CoinPrePostsolveMatrix::atUpperBound) ; /* PRESOLVEASSERT(prob->getColumnStatus(jcol)!=CoinPrePostsolveMatrix::basic); if (clo[jcol] <= -PRESOLVE_INF) { CoinPrePostsolveMatrix::Status statj = prob->getColumnStatus(jcol) ; if (statj == CoinPrePostsolveMatrix::atLowerBound) { if (cup[jcol] < PRESOLVE_INF) { statj = CoinPrePostsolveMatrix::atUpperBound ; } else { statj = CoinPrePostsolveMatrix::isFree ; } prob->setColumnStatus(jcol,statj) ; } } */ } PRESOLVEASSERT(prob->getRowStatus(irow)==CoinPrePostsolveMatrix::basic); PRESOLVEASSERT(rowduals[irow] == 0.0); // this is a lazy implementation. // we tightened the col bounds, then let them be eliminated // by repeated uses of FIX_VARIABLE and a final DROP_ROW. // Therefore, by this point the row has been marked basic, // the rowdual of this row is 0.0, // and the reduced costs for the cols may or may not be ok // for the relaxed column bounds. // // find the one most out of whack and fix it. int whacked = -1; double whack = 0.0; for (k=0; k ztoldj && !(fabs(sol[jcol] - clo[jcol]) <= ztolzb)) || (rcosts[jcol] < -ztoldj && !(fabs(sol[jcol] - cup[jcol]) <= ztolzb))) && fabs(whack0) > fabs(whack)) { whacked = jcol; whack = whack0; } } if (whacked != -1) { prob->setColumnStatus(whacked,CoinPrePostsolveMatrix::basic); prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound); rowduals[irow] = whack; for (k=0; kfeasibilityTolerance_; for (int irow=0; irow 0.0) { if (PRESOLVE_INF <= ub) { if (ub_inf_index == -1) { ub_inf_index = k; } else { ub_inf_index = -2; if (lb_inf_index == -2) break; // pointless } } else maxup += ub * coeff; if (lb <= -PRESOLVE_INF) { if (lb_inf_index == -1) { lb_inf_index = k; } else { lb_inf_index = -2; if (ub_inf_index == -2) break; // pointless } } else maxdown += lb * coeff; } else { if (PRESOLVE_INF <= ub) { if (lb_inf_index == -1) { lb_inf_index = k; } else { lb_inf_index = -2; if (ub_inf_index == -2) break; // pointless } } else maxdown += ub * coeff; if (lb <= -PRESOLVE_INF) { if (ub_inf_index == -1) { ub_inf_index = k; } else { ub_inf_index = -2; if (lb_inf_index == -2) break; // pointless } } else maxup += lb * coeff; } } // ub_inf says whether the sum of the "other" ub terms is infinite // in the loop below. // In the case where we only saw one infinite term, the loop // will only cover that case, in which case the other terms // are *not* infinite. // With two or more terms, it is infinite. // If we only saw one infinite term, then if (ub_inf_index == -2) maxup = PRESOLVE_INF; if (lb_inf_index == -2) maxdown = -PRESOLVE_INF; const bool maxup_finite = PRESOLVEFINITE(maxup); const bool maxdown_finite = PRESOLVEFINITE(maxdown); if (ub_inf_index == -1 && maxup_finite && maxup + tol < rlo[irow]&&!fixInfeasibility) { /* infeasible */ prob->status_|= 1; prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, prob->messages()) <status_|= 1; prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS, prob->messages()) < ZTOLDP && !integerType[jcol]) { double maxup1 = (ub_inf_index == -1 || ub_inf_index == k ? maxup : PRESOLVE_INF); bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k ? maxup_finite : false); double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k ? maxdown : PRESOLVE_INF); bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k ? maxdown_finite : false); double ilb = (irlo - maxup1) / coeff; bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1); double iub = (irup - maxdown1) / coeff; bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1); double ilb1 = (coeff > 0.0 ? (finite_ilb ? ilb : -PRESOLVE_INF) : (finite_iub ? iub : -PRESOLVE_INF)); if (ilbound[jcol] < ilb1) { ilbound[jcol] = ilb1; if (jcol == 278001) printf("JCOL LB %g\n", ilb1); } } } for (k = krs; k ZTOLDP && !integerType[jcol]) { double maxup1 = (ub_inf_index == -1 || ub_inf_index == k ? maxup : PRESOLVE_INF); bool maxup_finite1 = (ub_inf_index == -1 || ub_inf_index == k ? maxup_finite : false); double maxdown1 = (lb_inf_index == -1 || lb_inf_index == k ? maxdown : PRESOLVE_INF); bool maxdown_finite1 = (ub_inf_index == -1 || ub_inf_index == k ? maxdown_finite : false); double ilb = (irlo - maxup1) / coeff; bool finite_ilb = (-PRESOLVE_INF < irlo && maxup_finite1); double iub = (irup - maxdown1) / coeff; bool finite_iub = ( irup < PRESOLVE_INF && maxdown_finite1); double iub1 = (coeff > 0.0 ? (finite_iub ? iub : PRESOLVE_INF) : (finite_ilb ? ilb : PRESOLVE_INF)); if (iub1 < iubound[jcol]) { iubound[jcol] = iub1; if (jcol == 278001) printf("JCOL UB %g\n", iub1); } } } } } #if 0 // (B) postsolve for implied_bound { double lo0 = pa->clo; double up0 = pa->cup; int irow = pa->irow; int jcol = pa->icol; int *rowcols = pa->rowcols; int ninrow = pa->ninrow; clo[jcol] = lo0; cup[jcol] = up0; if ((colstat[jcol] & PRESOLVE_XBASIC) == 0 && fabs(lo0 - sol[jcol]) > ztolzb && fabs(up0 - sol[jcol]) > ztolzb) { // this non-basic variable is now away from its bound // it is ok just to force it to be basic // informally: if this variable is at its implied bound, // then the other variables must be at their bounds, // which means the bounds will stop them even if the aren't basic. if (rowstat[irow] & PRESOLVE_XBASIC) rowstat[irow] = 0; else { int k; for (k=0; k= -ztoldj) || (fabs(cup[col] - sol[col]) <= ztolzb && rcosts[col] <= ztoldj))) break; } if (k %d\n", irow, col, jcol); #endif colstat[col] = 0; // since all vars were at their bounds, the slack must be 0 PRESOLVEASSERT(fabs(acts[irow]) < ZTOLDP); rowstat[irow] = PRESOLVE_XBASIC; } else { // should never happen? abort(); } // get rid of any remaining basic structurals, since their rcosts // are going to become non-zero in a second. abort(); ///////////////////zero_pivot(); } double rdual_adjust; { CoinBigIndex kk = presolve_find_row(irow, mcstrt[jcol], mcstrt[jcol] + hincol[jcol], hrow); // adjust rowdual to cancel out reduced cost // should probably search for col with largest factor rdual_adjust = (rcosts[jcol] / colels[kk]); rowduals[irow] += rdual_adjust; colstat[jcol] = PRESOLVE_XBASIC; } for (k=0; k