// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <stdio.h>
#include <math.h>
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
#include "CoinPresolvePsdebug.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveZeros.hpp"
#include "CoinPresolveSubst.hpp"
#include "CoinMessage.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinSort.hpp"
#include "CoinError.hpp"
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
#include "CoinPresolvePsdebug.hpp"
#endif
namespace { // begin unnamed file-local namespace
inline void prepend_elem(int jcol, double coeff, int irow,
CoinBigIndex *mcstrt,
double *colels,
int *hrow,
int *link, CoinBigIndex *free_listp)
{
CoinBigIndex kk = *free_listp;
assert(kk >= 0) ;
*free_listp = link[*free_listp];
link[kk] = mcstrt[jcol];
mcstrt[jcol] = kk;
colels[kk] = coeff;
hrow[kk] = irow;
}
// add coeff_factor * rowy to rowx
static bool add_row(CoinBigIndex *mrstrt,
double *rlo, double * acts, double *rup,
double *rowels,
int *hcol,
int *hinrow,
presolvehlink *rlink, int nrows,
double coeff_factor,
int irowx, int irowy,
int *x_to_y)
{
CoinBigIndex krs = mrstrt[irowy];
CoinBigIndex kre = krs + hinrow[irowy];
CoinBigIndex krsx = mrstrt[irowx];
CoinBigIndex krex = krsx + hinrow[irowx];
// const int maxk = mrstrt[nrows]; // (22)
// if irowx is very long, the searching gets very slow,
// so we always sort.
// whatever sorts rows should handle almost-sorted data efficiently
// (quicksort may not)
CoinSort_2(hcol+krsx,hcol+krsx+hinrow[irowx],rowels+krsx);
CoinSort_2(hcol+krs,hcol+krs+hinrow[irowy],rowels+krs);
//ekk_sort2(hcol+krsx, rowels+krsx, hinrow[irowx]);
//ekk_sort2(hcol+krs, rowels+krs, hinrow[irowy]);
//printf("%s x=%d y=%d cf=%g nx=%d ny=%d\n",
// "ADD_ROW:",
// irowx, irowy, coeff_factor, hinrow[irowx], hinrow[irowy]);
# if PRESOLVE_DEBUG
printf("%s x=%d y=%d cf=%g nx=%d ycols=(",
"ADD_ROW:",
irowx, irowy, coeff_factor, hinrow[irowx]);
# endif
// adjust row bounds of rowx;
// analogous to adjusting bounds info of colx in doubleton,
// or perhaps adjustment to rlo/rup in elim_doubleton
//
// I believe that since we choose a column that is implied free,
// no other column bounds need to be updated.
// This is what would happen in doubleton if y's bounds were implied free;
// in that case,
// lo1 would never improve clo[icolx] and
// up1 would never improve cup[icolx].
{
double rhsy = rlo[irowy];
// (1)
if (-PRESOLVE_INF < rlo[irowx]) {
# if PRESOLVE_DEBUG
if (rhsy * coeff_factor)
printf("ELIM_ROW RLO: %g -> %g\n",
rlo[irowx],
rlo[irowx] + rhsy * coeff_factor);
# endif
rlo[irowx] += rhsy * coeff_factor;
}
// (2)
if (rup[irowx] < PRESOLVE_INF) {
# if PRESOLVE_DEBUG
if (rhsy * coeff_factor)
printf("ELIM_ROW RUP: %g -> %g\n",
rup[irowx],
rup[irowx] + rhsy * coeff_factor);
# endif
rup[irowx] += rhsy * coeff_factor;
}
if (acts)
{ acts[irowx] += rhsy * coeff_factor ; }
}
CoinBigIndex kcolx = krsx;
CoinBigIndex krex0 = krex;
int x_to_y_i = 0;
for (CoinBigIndex krowy=krs; krowy<kre; krowy++) {
int jcol = hcol[krowy];
// even though these values are updated, they remain consistent
PRESOLVEASSERT(krex == krsx + hinrow[irowx]);
// see if row appears in colx
// do NOT look beyond the original elements of rowx
//CoinBigIndex kcolx = presolve_find_col1(jcol, krsx, krex, hcol);
while (kcolx < krex0 && hcol[kcolx] < jcol)
kcolx++;
# if PRESOLVE_DEBUG
printf("%d%s ", jcol, (kcolx < krex0 && hcol[kcolx] == jcol) ? "+" : "");
# endif
if (kcolx < krex0 && hcol[kcolx] == jcol) {
// before: both x and y are in the jcol
// after: only x is in the jcol
// so: number of elems in col x unchanged, and num elems in jcol is one less
// update row rep - just modify coefficent
// column y is deleted as a whole at the end of the loop
# if PRESOLVE_DEBUG
printf("CHANGING %g + %g -> %g\n",
rowels[kcolx],
rowels[krowy],
rowels[kcolx] + rowels[krowy] * coeff_factor);
# endif
rowels[kcolx] += rowels[krowy] * coeff_factor;
// this is where this element in rowy ended up
x_to_y[x_to_y_i++] = kcolx - krsx;
kcolx++;
} else {
// before: only y is in the jcol
// after: only x is in the jcol
// so: number of elems in col x is one greater, but num elems in jcol remains same
{
bool outOfSpace=presolve_expand_row(mrstrt,rowels,hcol,hinrow,rlink,nrows,irowx) ;
if (outOfSpace)
return true;
// this may force a compaction
// this will be called excessively if the rows are packed too tightly
// have to adjust various induction variables
krowy = mrstrt[irowy] + (krowy - krs);
krs = mrstrt[irowy]; // do this for ease of debugging
kre = mrstrt[irowy] + hinrow[irowy];
kcolx = mrstrt[irowx] + (kcolx - krsx); // don't really need to do this
krex0 = mrstrt[irowx] + (krex0 - krsx);
krsx = mrstrt[irowx];
krex = mrstrt[irowx] + hinrow[irowx];
}
// this is where this element in rowy ended up
x_to_y[x_to_y_i++] = krex - krsx;
// there is now an unused entry in the memory after the column - use it
// mrstrt[nrows] == penultimate index of arrays hcol/rowels
hcol[krex] = jcol;
rowels[krex] = rowels[krowy] * coeff_factor;
hinrow[irowx]++, krex++; // expand the col
// do NOT increment kcolx
}
}
# if PRESOLVE_DEBUG
printf(")\n");
# endif
return false;
}
// It is common in osl to copy from one representation to another
// (say from a col rep to a row rep).
// One such routine is ekkclcp.
// This is similar, except that it does not assume that the
// representation is packed, and it adds some slack space
// in the target rep.
// It assumes both hincol/hinrow are correct.
// Note that such routines automatically sort the target rep by index,
// because they sweep the rows in ascending order.
void copyrep(const int * mrstrt, const int *hcol, const double *rowels,
const int *hinrow, int nrows,
int *mcstrt, int *hrow, double *colels,
int *hincol, int ncols)
{
int pos = 0;
for (int j = 0; j < ncols; ++j) {
mcstrt[j] = pos;
pos += hincol[j];
pos += CoinMin(hincol[j], 10); // slack
hincol[j] = 0;
}
for (int i = 0; i < nrows; ++i) {
CoinBigIndex krs = mrstrt[i];
CoinBigIndex kre = krs + hinrow[i];
for (CoinBigIndex kr = krs; kr < kre; ++kr) {
int icol = hcol[kr];
int iput = hincol[icol];
hincol[icol] = iput + 1;
iput += mcstrt[icol];
hrow[iput] = i;
colels[iput] = rowels[kr];
}
}
}
} // end unnamed file-local namespace
const char *subst_constraint_action::name() const
{
return ("subst_constraint_action");
}
// add -x/y times row y to row x, thus cancelling out one column of rowx;
// afterwards, that col will be singleton for rowy, so we drop the row.
//
// This no longer maintains the col rep as it goes along.
// Instead, it reconstructs it from scratch afterward.
//
// This implements the functionality of ekkrdc3.
/*
This routine is called only from implied_free_action. There are several
oddities and redundancies in the relationship. The two routines need a good
grooming.
try_fill_level limits the allowable number of coefficients in a column
under consideration for substitution. There's some sort of hack going on
that has the following effect: if try_fill_level comes in as 2, and that
seems overly limiting (number of substitutions < 30), try increasing it to
3. To trigger a wider examination of columns, this is actually passed back
as -3. The next entry of implied_free_action (and then this routine) will
override ColsToDo and examine all columns.
Hence the initial loop triggered when try_fill_level < 0. Other positive
values of fill_level will have no effect. A value of -3 will be converted
(and passed back out) as +3. Arbitrary negative values of try_fill_level
will also trigger the expansion of search and be converted to positive
values.
I would have thought that the columns considered by implied_free_action
should also be limited by fill_level, but that's not currently the case.
It's hard-wired to consider columns with 1 to 3 coefficients.
There must be a better way. -- lh, 040818 --
*/
const CoinPresolveAction *subst_constraint_action::presolve(CoinPresolveMatrix *prob,
int *implied_free,
const CoinPresolveAction *next,
int &try_fill_level)
{
double *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
const int ncols = prob->ncols_;
double *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
const int nrows = prob->nrows_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
double *acts = prob->acts_;
double *dcost = prob->cost_;
presolvehlink *clink = prob->clink_;
presolvehlink *rlink = prob->rlink_;
const double tol = prob->feasibilityTolerance_;
action *actions = new action [ncols];
# ifdef ZEROFAULT
CoinZeroN(reinterpret_cast<char *>(actions),ncols*sizeof(action)) ;
# endif
int nactions = 0;
int *zerocols = new int[ncols];
int nzerocols = 0;
int *x_to_y = new int[ncols];
#if 0
// follmer.mps presents a challenge, since it has some very
// long rows. I started experimenting with how to deal with it,
// but haven't yet finished.
// The idea was to space out the rows to add some padding between them.
// Ideally, we shouldn't have to do this just here, but could try to
// do it a little everywhere.
// sort the row rep by reconstructing from col rep
copyrep(mcstrt, hrow, colels, hincol, ncols,
mrstrt, hcol, rowels, hinrow, nrows);
presolve_make_memlists(mrstrt, hinrow, rlink, nrows);
// NEED SOME ASSERTION ABOUT NELEMS
copyrep(mrstrt, hcol, rowels, hinrow, nrows,
mcstrt, hrow, colels, hincol, ncols);
presolve_make_memlists(mcstrt, hincol, clink, ncols);
#endif
// in the original presolve, I don't think the two representations were
// kept in sync. It may be useful not to do that here, either,
// but rather just keep the columns with nfill_level rows accurate
// and resync at the end of the function.
// DEBUGGING
#if PRESOLVE_DEBUGx
int maxsubst = atoi(getenv("MAXSUBST"));
#else
const int maxsubst = 1000000;
#endif
int nsubst = 0;
// This loop does very nearly the same thing as
// the first loop in implied_free_action::presolve.
// We have to do it again in case constraints change while we process
// them (???).
/*
No --- given the hack with -3 coming in to implied_free_action and overriding
ColsToDo, we could have columns in implied_free that aren't in ColsToDo.
-- lh, 040818 --
*/
int numberLook = prob->numberColsToDo_;
int iLook;
int * look = prob->colsToDo_;
int fill_level = try_fill_level;
int * look2 = NULL;
// if gone from 2 to 3 look at all
if (fill_level<0) {
fill_level=-fill_level;
try_fill_level=fill_level;
look2 = new int[ncols];
look=look2;
if (!prob->anyProhibited()) {
for (iLook=0;iLook<ncols;iLook++)
look[iLook]=iLook;
numberLook=ncols;
} else {
// some prohibited
numberLook=0;
for (iLook=0;iLook<ncols;iLook++)
if (!prob->colProhibited(iLook))
look[numberLook++]=iLook;
}
}
for (iLook=0;iLook<numberLook;iLook++) {
int jcoly=look[iLook];
bool looksGood=false;
if (hincol[jcoly] > 1 && hincol[jcoly] <= fill_level &&
implied_free[jcoly] >=0) {
looksGood=true;
CoinBigIndex kcs = mcstrt[jcoly];
CoinBigIndex kce = kcs + hincol[jcoly];
int bestrowy_size = 0;
int bestrowy_row=-1;
int bestrowy_k=-1;
double bestrowy_coeff=0.0;
CoinBigIndex k;
for (k=kcs; k<kce; ++k) {
int row = hrow[k];
double coeffj = colels[k];
// we don't clean up zeros in the middle of the routine.
// if there is one, skip this candidate.
if (fabs(coeffj) <= ZTOLDP) {
bestrowy_size = 0;
break;
}
if (row==implied_free[jcoly]) {
// if its row is an equality constraint...
if (hinrow[row] > 1 && // don't bother with singleton rows
fabs(rlo[row] - rup[row]) < tol &&
!prob->rowUsed(row)) {
// both column bounds implied by the constraint bounds
// we want coeffy to be smaller than x, BACKWARDS from in doubleton
bestrowy_size = hinrow[row];
bestrowy_row = row;
bestrowy_coeff = coeffj;
bestrowy_k = k;
}
}
}
if (bestrowy_size == 0)
continue;
bool all_ok = true;
for (k=kcs; k<kce; ++k) {
double coeff_factor = fabs(colels[k] / bestrowy_coeff);
if (fabs(coeff_factor) > 10.0)
all_ok = false;
}
#if 0 // block A
// check fill-in
if (all_ok && hincol[jcoly] == 3) {
// compute fill-in
int row1 = -1;
int row2=-1;
CoinBigIndex kk;
for (kk=kcs; kk<kce; ++kk)
if (kk != bestrowy_k) {
if (row1 == -1)
row1 = hrow[kk];
else
row2 = hrow[kk];
}
CoinBigIndex krs = mrstrt[bestrowy_row];
CoinBigIndex kre = krs + hinrow[bestrowy_row];
CoinBigIndex krs1 = mrstrt[row1];
CoinBigIndex kre1 = krs + hinrow[row1];
CoinBigIndex krs2 = mrstrt[row2];
CoinBigIndex kre2 = krs + hinrow[row2];
CoinSort_2(hcol+krs,hcol+krs+hinrow[bestrowy_row],rowels+krs);
CoinSort_2(hcol+krs1,hcol+krs1+hinrow[row1],rowels+krs1);
CoinSort_2(hcol+krs2,hcol+krs2+hinrow[row2],rowels+krs2);
//ekk_sort2(hcol+krs, rowels+krs, hinrow[bestrowy_row]);
//ekk_sort2(hcol+krs1, rowels+krs1, hinrow[row1]);
//ekk_sort2(hcol+krs2, rowels+krs2, hinrow[row2]);
int nfill = -hinrow[bestrowy_row];
CoinBigIndex kcol1 = krs1;
for (kk=krs; kk<kre; ++kk) {
int jcol = hcol[kk];
while (kcol1 < kre1 && hcol[kcol1] < jcol)
kcol1++;
if (! (kcol1 < kre1 && hcol[kcol1] == jcol))
nfill++;
}
CoinBigIndex kcol2 = krs2;
for (kk=krs; kk<kre; ++kk) {
int jcol = hcol[kk];
while (kcol2 < kre2 && hcol[kcol2] < jcol)
kcol2++;
if (! (kcol2 < kre2 && hcol[kcol2] == jcol))
nfill++;
}
#if PRESOLVE_DEBUG
printf("FILL: %d\n", nfill);
#endif
#if 0
static int maxfill = atoi(getenv("MAXFILL"));
if (nfill > maxfill)
all_ok = false;
#endif
// not too much
if (nfill <= 0)
ngood++;
#if 0
static int nts = 0;
if (++nts > atoi(getenv("NTS")))
all_ok = false;
else
nt++;
#endif
}
#endif // end block A
// probably never happens
if (all_ok && nzerocols + hinrow[bestrowy_row] >= ncols)
all_ok = false;
if (nsubst >= maxsubst) {
all_ok = false;
}
if (all_ok) {
nsubst++;
#if 0
// debug
if (numberLook<ncols&&iLook==numberLook-1) {
printf("found last one?? %d\n", jcoly);
}
#endif
CoinBigIndex kcs = mcstrt[jcoly];
int rowy = bestrowy_row;
double coeffy = bestrowy_coeff;
PRESOLVEASSERT(fabs(colels[kcs]) > ZTOLDP);
PRESOLVEASSERT(fabs(colels[kcs+1]) > ZTOLDP);
PRESOLVEASSERT(hinrow[rowy] > 1);
const bool nonzero_cost = (fabs(dcost[jcoly]) > tol);
double *costsx = nonzero_cost ? new double[hinrow[rowy]] : 0;
int ntotels = 0;
for (k=kcs; k<kce; ++k) {
int irow = hrow[k];
ntotels += hinrow[irow];
// mark row as contaminated
prob->setRowUsed(irow);
}
{
action *ap = &actions[nactions++];
int nincol = hincol[jcoly];
ap->col = jcoly;
ap->rowy = rowy;
ap->nincol = nincol;
ap->rows = new int[nincol];
ap->rlos = new double[nincol];
ap->rups = new double[nincol];
// coefficients in deleted col
ap->coeffxs = new double[nincol];
ap->ninrowxs = new int[nincol];
ap->rowcolsxs = new int[ntotels];
ap->rowelsxs = new double[ntotels];
ap->costsx = costsx;
// copy all the rows for restoring later - wasteful
{
int nel = 0;
for (CoinBigIndex k=kcs; k<kce; ++k) {
int irow = hrow[k];
CoinBigIndex krs = mrstrt[irow];
// CoinBigIndex kre = krs + hinrow[irow];
prob->addRow(irow);
ap->rows[k-kcs] = irow;
ap->ninrowxs[k-kcs] = hinrow[irow];
ap->rlos[k-kcs] = rlo[irow];
ap->rups[k-kcs] = rup[irow];
ap->coeffxs[k-kcs] = colels[k];
memcpy( &ap->rowcolsxs[nel], &hcol[krs],hinrow[irow]*sizeof(int));
memcpy( &ap->rowelsxs[nel], &rowels[krs],hinrow[irow]*sizeof(double));
nel += hinrow[irow];
}
}
}
// rowy is supposed to be an equality row
PRESOLVEASSERT(fabs(rup[rowy] - rlo[rowy]) < ZTOLDP);
// now adjust for the implied free row - COPIED
if (nonzero_cost) {
#if 0&&PRESOLVE_DEBUG
printf("NONZERO SUBST COST: %d %g\n", jcoly, dcost[jcoly]);
#endif
double *cost = dcost;
double *save_costs = costsx;
double coeffj = coeffy;
CoinBigIndex krs = mrstrt[rowy];
CoinBigIndex kre = krs + hinrow[rowy];
double rhs = rlo[rowy];
double costj = cost[jcoly];
for (CoinBigIndex k=krs; k<kre; k++) {
int jcol = hcol[k];
prob->addCol(jcol);
save_costs[k-krs] = cost[jcol];
if (jcol != jcoly) {
double coeff = rowels[k];
/*
* Similar to eliminating doubleton:
* cost1 x = cost1 (c - b y) / a = (c cost1)/a - (b cost1)/a
* cost[icoly] += cost[icolx] * (-coeff2 / coeff1);
*/
cost[jcol] += costj * (-coeff / coeffj);
}
}
// I'm not sure about this
prob->change_bias(costj * rhs / coeffj);
// ??
cost[jcoly] = 0.0;
}
#if 0&&PRESOLVE_DEBUG
if (hincol[jcoly] == 3) {
CoinBigIndex krs = mrstrt[rowy];
CoinBigIndex kre = krs + hinrow[rowy];
printf("HROW0 (%d): ", rowy);
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
double coeff = rowels[k];
printf("%d:%g (%d) ", jcol, coeff, hincol[jcol]);
}
printf("\n");
}
#endif
if (hincol[jcoly] != 2) {
CoinBigIndex krs = mrstrt[rowy];
// CoinBigIndex kre = krs + hinrow[rowy];
CoinSort_2(hcol+krs,hcol+krs+hinrow[rowy],rowels+krs);
//ekk_sort2(hcol+krs, rowels+krs, hinrow[rowy]);
}
// substitute away jcoly in the other rows
// Use ap as mcstrt etc may move if compacted
kce = hincol[jcoly];
action *ap = &actions[nactions-1];
for (k=0; k<kce; ++k) {
int rowx = ap->rows[k];
//assert(rowx==hrow[k+kcs]);
//assert(ap->coeffxs[k]==colels[k+kcs]);
if (rowx != rowy) {
double coeffx = ap->coeffxs[k];
double coeff_factor = -coeffx / coeffy; // backwards from doubleton
#if 0&&PRESOLVE_DEBUG
{
CoinBigIndex krs = mrstrt[rowx];
CoinBigIndex kre = krs + hinrow[rowx];
printf("HROW (%d %d %d): ", rowx, hinrow[rowx], jcoly);
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
double coeff = rowels[k];
printf("%d ", jcol);
}
printf("\n");
#if 0
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
prob->addCol(jcol);
double coeff = rowels[k];
printf("%g ", coeff);
}
printf("\n");
#endif
}
#endif
{
CoinBigIndex krsx = mrstrt[rowx];
CoinBigIndex krex = krsx + hinrow[rowx];
int i;
for (i=krsx;i<krex;i++)
prob->addCol(hcol[i]);
if (hincol[jcoly] != 2)
CoinSort_2(hcol+krsx,hcol+krsx+hinrow[rowx],rowels+krsx);
//ekk_sort2(hcol+krsx, rowels+krsx, hinrow[rowx]);
}
// add (coeff_factor * <rowy>) to rowx
// does not affect rowy
// may introduce (or cancel) elements in rowx
bool outOfSpace = add_row(mrstrt,
rlo, acts, rup,
rowels, hcol,
hinrow,
rlink, nrows,
coeff_factor,
rowx, rowy,
x_to_y);
if (outOfSpace)
throwCoinError("out of memory",
"CoinImpliedFree::presolve");
// update col rep of rowx from row rep:
// for every col in rowy, copy the elem for that col in rowx
// from the row rep to the col rep
{
CoinBigIndex krs = mrstrt[rowy];
// CoinBigIndex kre = krs + hinrow[rowy];
int niny = hinrow[rowy];
CoinBigIndex krsx = mrstrt[rowx];
// CoinBigIndex krex = krsx + hinrow[rowx];
for (CoinBigIndex ki=0; ki<niny; ++ki) {
CoinBigIndex k = krs + ki;
int jcol = hcol[k];
prob->addCol(jcol);
CoinBigIndex kcs = mcstrt[jcol];
CoinBigIndex kce = kcs + hincol[jcol];
//double coeff = rowels[presolve_find_col(jcol, krsx, krex, hcol)];
if (hcol[krsx + x_to_y[ki]] != jcol)
abort();
double coeff = rowels[krsx + x_to_y[ki]];
// see if rowx appeared in jcol in the col rep
CoinBigIndex k2 = presolve_find_row1(rowx, kcs, kce, hrow);
//PRESOLVEASSERT(fabs(coeff) > ZTOLDP);
if (k2 < kce) {
// yes - just update the entry
colels[k2] = coeff;
} else {
// no - make room, then append
bool outOfSpace=presolve_expand_row(mcstrt,colels,hrow,hincol,
clink,ncols,jcol) ;
if (outOfSpace)
throwCoinError("out of memory",
"CoinImpliedFree::presolve");
krsx = mrstrt[rowx];
krs = mrstrt[rowy];
kcs = mcstrt[jcol];
kce = kcs + hincol[jcol];
hrow[kce] = rowx;
colels[kce] = coeff;
hincol[jcol]++;
}
}
}
// now colels[k] == 0.0
#if 1
// now remove jcoly from rowx in the row rep
// better if this were first
presolve_delete_from_row(rowx, jcoly, mrstrt, hinrow, hcol, rowels);
#endif
#if 0&&PRESOLVE_DEBUG
{
CoinBigIndex krs = mrstrt[rowx];
CoinBigIndex kre = krs + hinrow[rowx];
printf("HROW (%d %d %d): ", rowx, hinrow[rowx], jcoly);
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
double coeff = rowels[k];
printf("%d ", jcol);
}
printf("\n");
#if 0
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
double coeff = rowels[k];
printf("%g ", coeff);
}
printf("\n");
#endif
}
#endif
// don't have to update col rep, since entire col deleted later
}
}
#if 0&&PRESOLVE_DEBUG
printf("\n");
#endif
// the addition of rows may have created zero coefficients
memcpy( &zerocols[nzerocols], &hcol[mrstrt[rowy]],hinrow[rowy]*sizeof(int));
nzerocols += hinrow[rowy];
// delete rowy in col rep
{
CoinBigIndex krs = mrstrt[rowy];
CoinBigIndex kre = krs + hinrow[rowy];
for (CoinBigIndex k=krs; k<kre; ++k) {
int jcol = hcol[k];
// delete rowy from the jcol
presolve_delete_from_col(rowy,jcol,mcstrt,hincol,hrow,colels) ;
if (hincol[jcol] == 0)
{ PRESOLVE_REMOVE_LINK(clink,jcol) ; }
}
}
// delete rowy in row rep
hinrow[rowy] = 0;
// This last is entirely dual to doubleton, but for the cost adjustment
// eliminate col entirely from the col rep
PRESOLVE_REMOVE_LINK(clink, jcoly);
hincol[jcoly] = 0;
// eliminate rowy entirely from the row rep
PRESOLVE_REMOVE_LINK(rlink, rowy);
//cost[irowy] = 0.0;
rlo[rowy] = 0.0;
rup[rowy] = 0.0;
#if 0 && PRESOLVE_DEBUG
printf("ROWY COLS: ");
for (CoinBigIndex k=0; k<save_ninrowy; ++k)
if (rowycols[k] != col) {
printf("%d ", rowycols[k]);
(void)presolve_find_col(rowycols[k], mrstrt[rowx], mrstrt[rowx]+hinrow[rowx],
hcol);
}
printf("\n");
#endif
# if PRESOLVE_CONSISTENCY
presolve_links_ok(prob) ;
presolve_consistent(prob) ;
# endif
}
}
}
// Clear row used flags
for (int iRow=0;iRow<nrows;iRow++)
prob->unsetRowUsed(iRow);
// general idea - only do doubletons until there are almost none left
if (nactions < 30&&fill_level<prob->maxSubstLevel_)
try_fill_level = -fill_level-1;
if (nactions) {
# if PRESOLVE_SUMMARY
printf("NSUBSTS: %d\n", nactions);
//printf("NT: %d NGOOD: %d FILL_LEVEL: %d\n", nt, ngood, fill_level);
# endif
next = new subst_constraint_action(nactions, CoinCopyOfArray(actions,nactions), next);
next = drop_zero_coefficients_action::presolve(prob, zerocols, nzerocols, next);
}
delete [] look2;
deleteAction(actions,action*);
delete[]x_to_y;
delete[]zerocols;
return (next);
}
void subst_constraint_action::postsolve(CoinPostsolveMatrix *prob) const
{
const action *const actions = actions_;
const int nactions = nactions_;
double *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *link = prob->link_;
// int ncols = prob->ncols_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
double *dcost = prob->cost_;
double *sol = prob->sol_;
double *rcosts = prob->rcosts_;
double *acts = prob->acts_;
double *rowduals = prob->rowduals_;
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
char *cdone = prob->cdone_;
char *rdone = prob->rdone_;
# endif
CoinBigIndex &free_list = prob->free_list_;
// const double ztoldj = prob->ztoldj_;
const double maxmin = prob->maxmin_;
int k;
for (const action *f = &actions[nactions-1]; actions<=f; f--) {
int icol = f->col;
int nincoly = f->nincol;
double *rlos = f->rlos;
double *rups = f->rups;
int *rows = f->rows;
double *coeffxs = f->coeffxs;
int jrowy = f->rowy;
int *ninrowxs = f->ninrowxs;
const int *rowcolsxs = f->rowcolsxs;
const double *rowelsxs = f->rowelsxs;
/* the row was in the reduced problem */
for (int i=0; i<nincoly; ++i) {
if (rows[i] != jrowy)
PRESOLVEASSERT(rdone[rows[i]]);
}
PRESOLVEASSERT(cdone[icol]==DROP_COL);
PRESOLVEASSERT(rdone[jrowy]==DROP_ROW);
// DEBUG CHECK
#if 0 && PRESOLVE_DEBUG
{
double actx = 0.0;
for (int j=0; j<ncols; ++j)
if (hincol[j] > 0 && cdone[j]) {
CoinBigIndex krow = presolve_find_row1(jrowx, mcstrt[j], mcstrt[j] + hincol[j], hrow);
if (krow < mcstrt[j] + hincol[j])
actx += colels[krow] * sol[j];
}
if (! (fabs(acts[jrowx] - actx) < 100*ztolzb))
printf("BAD ACTSX: acts[%d]==%g != %g\n",
jrowx, acts[jrowx], actx);
if (! (rlo[jrowx] - 100*ztolzb <= actx && actx <= rup[jrowx] + 100*ztolzb))
printf("ACTSX NOT IN RANGE: %d %g %g %g\n",
jrowx, rlo[jrowx], actx, rup[jrowx]);
}
#endif
int ninrowy=-1;
const int *rowcolsy=NULL;
const double *rowelsy=NULL;
double coeffy=0.0;
double rloy=1.0e50;
{
int nel = 0;
for (int i=0; i<nincoly; ++i) {
int row = rows[i];
rlo[row] = rlos[i];
rup[row] = rups[i];
if (row == jrowy) {
ninrowy = ninrowxs[i];
rowcolsy = &rowcolsxs[nel];
rowelsy = &rowelsxs[nel];
coeffy = coeffxs[i];
rloy = rlo[row];
}
nel += ninrowxs[i];
}
}
double rhsy = rloy;
// restore costs
{
const double *costs = f->costsx;
if (costs)
for (int i = 0; i<ninrowy; ++i) {
dcost[rowcolsy[i]] = costs[i];
}
}
// solve for the equality to find the solution for the eliminated col
// this is why we want coeffx < coeffy (55)
{
double sol0 = rloy;
sol[icol] = 0.0; // to avoid condition in loop
for (k = 0; k<ninrowy; ++k) {
int jcolx = rowcolsy[k];
double coeffx = rowelsy[k];
sol0 -= coeffx * sol[jcolx];
}
sol[icol] = sol0 / coeffy;
# if PRESOLVE_DEBUG
const double ztolzb = prob->ztolzb_;
double *clo = prob->clo_;
double *cup = prob->cup_;
if (! (sol[icol] > clo[icol] - ztolzb &&
cup[icol] + ztolzb > sol[icol]))
printf("NEW SOL OUT-OF-TOL: %g %g %g\n", clo[icol],
sol[icol], cup[icol]);
# endif
}
// since this row is fixed
acts[jrowy] = rloy;
// acts[irow] always ok, since slack is fixed
prob->setRowStatus(jrowy,CoinPrePostsolveMatrix::atLowerBound);
// remove old rowx from col rep
// we don't explicitly store what the current rowx is;
// however, after the presolve, rowx contains a col for every
// col in either the original rowx or the original rowy.
// If there were cancellations, those were handled in subsequent
// presolves.
{
// erase those cols in the other rows that occur in rowy
// (with the exception of icol, which was deleted);
// the other rows *must* contain these cols
for (k = 0; k<ninrowy; ++k) {
int col = rowcolsy[k];
// remove jrowx from col in the col rep
// icol itself was deleted, so won't be there
if (col != icol)
for (int i = 0; i<nincoly; ++i) {
if (rows[i] != jrowy)
presolve_delete_from_col2(rows[i],col,mcstrt,hincol,hrow,colels,
link,&free_list) ;
}
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
// initialize this for loops below
hincol[icol] = 0;
// now restore the original rows (other than rowy).
// those cols that were also in rowy were just removed;
// otherwise, they *must* already be there.
// This loop and the next automatically create the rep for the new col.
{
const int *rowcolsx = rowcolsxs;
const double *rowelsx = rowelsxs;
for (int i = 0; i<nincoly; ++i) {
int ninrowx = ninrowxs[i];
int jrowx = rows[i];
if (jrowx != jrowy)
for (k = 0; k<ninrowx; ++k) {
int col = rowcolsx[k];
CoinBigIndex kcolx = presolve_find_row3(jrowx, mcstrt[col], hincol[col], hrow, link);
if (kcolx != -1) {
PRESOLVEASSERT(presolve_find_col1(col, 0, ninrowy, rowcolsy) == ninrowy);
// overwrite the existing entry
colels[kcolx] = rowelsx[k];
} else {
PRESOLVEASSERT(presolve_find_col1(col, 0, ninrowy, rowcolsy) < ninrowy);
{
CoinBigIndex kk = free_list;
assert(kk >= 0 && kk < prob->bulk0_) ;
free_list = link[free_list];
link[kk] = mcstrt[col];
mcstrt[col] = kk;
colels[kk] = rowelsx[k];
hrow[kk] = jrowx;
}
++hincol[col];
}
}
rowcolsx += ninrowx;
rowelsx += ninrowx;
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
}
// finally, add original rowy elements
for (k = 0; k<ninrowy; ++k) {
int col = rowcolsy[k];
{
prepend_elem(col, rowelsy[k], jrowy, mcstrt, colels, hrow, link, &free_list);
++hincol[col];
}
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
}
// my guess is that the CLAIM in doubleton generalizes to
// equations with more than one x-style variable.
// Since I can't see how to distinguish among them,
// I assume that any of them will do.
{
// CoinBigIndex k;
double dj = maxmin*dcost[icol];
double bounds_factor = rhsy/coeffy;
for (int i=0; i<nincoly; ++i)
if (rows[i] != jrowy) {
int row = rows[i];
double coeff = coeffxs[i];
// PROBABLY DOESN'T MAKE SENSE
acts[row] += coeff * bounds_factor;
dj -= rowduals[row] * coeff;
}
// DEBUG CHECK
double acty = 0.0;
for (k = 0; k<ninrowy; ++k) {
int col = rowcolsy[k];
acty += rowelsy[k] * sol[col];
}
PRESOLVEASSERT(fabs(acty - acts[jrowy]) < 100*ZTOLDP);
// RECOMPUTING
{
const int *rowcolsx = rowcolsxs;
const double *rowelsx = rowelsxs;
for (int i=0; i<nincoly; ++i) {
int ninrowx = ninrowxs[i];
if (rows[i] != jrowy) {
int jrowx = rows[i];
double actx = 0.0;
for (k = 0; k<ninrowx; ++k) {
int col = rowcolsx[k];
actx += rowelsx[k] * sol[col];
}
PRESOLVEASSERT(rlo[jrowx] - prob->ztolzb_ <= actx
&& actx <= rup[jrowx] + prob->ztolzb_);
acts[jrowx] = actx;
}
rowcolsx += ninrowx;
rowelsx += ninrowx;
}
}
// this is the coefficient we need to force col y's reduced cost to 0.0;
// for example, this is obviously true if y is a singleton column
rowduals[jrowy] = dj / coeffy;
rcosts[icol] = 0.0;
// furthermore, this should leave rcosts[colx] for all colx
// in jrowx unchanged (????).
}
// Unlike doubleton, there should never be a problem with keeping
// the reduced costs the way they were, because the other
// variable's bounds are never changed, since col was implied free.
//rowstat[jrowy] = 0;
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::basic);
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
cdone[icol] = SUBST_ROW;
rdone[jrowy] = SUBST_ROW;
# endif
}
return ;
}
subst_constraint_action::~subst_constraint_action()
{
const action *actions = actions_;
for (int i=0; i<nactions_; ++i) {
delete[]actions[i].rows;
delete[]actions[i].rlos;
delete[]actions[i].rups;
delete[]actions[i].coeffxs;
delete[]actions[i].ninrowxs;
delete[]actions[i].rowcolsxs;
delete[]actions[i].rowelsxs;
//delete [](double*)actions[i].costsx;
deleteAction(actions[i].costsx,double*);
}
// Must add cast to placate MS compiler
//delete [] (subst_constraint_action::action*)actions_;
deleteAction(actions_,subst_constraint_action::action*);
}
syntax highlighted by Code2HTML, v. 0.9.1