// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <stdio.h>
#include <math.h>
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveDupcol.hpp"
#include "CoinSort.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinPresolveUseless.hpp"
#include "CoinMessage.hpp"
//#define PRESOLVE_DEBUG 1
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
#include "CoinPresolvePsdebug.hpp"
#endif
#define DSEED2 2147483647.0
namespace { // begin unnamed file-local namespace
void init_random_vec(double *work, int n)
{
double deseed = 12345678.0;
for (int i = 0; i < n; ++i) {
deseed *= 16807.;
int jseed = (int) (deseed / DSEED2);
deseed -= (double) jseed * DSEED2;
double random = deseed / DSEED2;
work[i]=random;
}
}
/*
For each candidate major-dimension vector in majcands, calculate the sum
over the vector, with each minor dimension weighted by a random amount.
(E.g., calculate column sums with each row weighted by a random amount.)
The sums are returned in the corresponding entries of majsums.
*/
void compute_sums (int n, const int *majlens, const CoinBigIndex *majstrts,
int *minndxs, double *elems, const double *minmuls,
int *majcands, double *majsums, int nlook)
{ for (int cndx = 0 ; cndx < nlook ; ++cndx)
{ int i = majcands[cndx] ;
PRESOLVEASSERT(majlens[i] > 0) ;
CoinBigIndex kcs = majstrts[i] ;
CoinBigIndex kce = kcs + majlens[i] ;
double value = 0.0 ;
for (CoinBigIndex k = kcs ; k < kce ; k++)
{ int irow = minndxs[k] ;
value += minmuls[irow]*elems[k] ; }
majsums[cndx] = value ; }
return ; }
void create_col (int col, int n, double *els,
CoinBigIndex *mcstrt, double *colels, int *hrow, int *link,
CoinBigIndex *free_listp)
{
int *rows = reinterpret_cast<int *>(els+n) ;
CoinBigIndex free_list = *free_listp;
int xstart = NO_LINK;
for (int i=0; i<n; ++i) {
CoinBigIndex k = free_list;
assert(k >= 0) ;
free_list = link[free_list];
hrow[k] = rows[i];
colels[k] = els[i];
link[k] = xstart;
xstart = k;
}
mcstrt[col] = xstart;
*free_listp = free_list;
}
} // end unnamed file-local namespace
const char *dupcol_action::name () const
{
return ("dupcol_action");
}
/*
Original comment: This is just ekkredc5, adapted into the new framework.
The datasets scorpion.mps and allgrade.mps have duplicate columns.
In case you don't have your OSL manual handy, a somewhat more informative
explanation: We're looking for an easy-to-detect special case of linearly
dependent columns, where the coefficients of the duplicate columns are
exactly equal. The idea for locating such columns is to generate pseudo-
random weights for each row and then calculate the weighted sum of
coefficients of each column. Columns with equal sums are checked more
thoroughly.
Analysis of the situation says there are two major cases:
* If the columns have equal objective coefficients, we can combine
them.
* If the columns have unequal objective coefficients, we may be able to
fix one at bound. If the required bound doesn't exist, we have dual
infeasibility (hence one of primal infeasibility or unboundedness).
In the comments below are a few fragments of code from the original
routine. I don't think they make sense, but I've left them for the nonce in
case someone else recognises the purpose. -- lh, 040909 --
*/
const CoinPresolveAction
*dupcol_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
double startTime = 0.0;
int startEmptyRows=0;
int startEmptyColumns = 0;
if (prob->tuning_) {
startTime = CoinCpuTime();
startEmptyRows = prob->countEmptyRows();
startEmptyColumns = prob->countEmptyCols();
}
double maxmin = prob->maxmin_ ;
double *colels = prob->colels_ ;
int *hrow = prob->hrow_ ;
CoinBigIndex *mcstrt = prob->mcstrt_ ;
int *hincol = prob->hincol_ ;
int ncols = prob->ncols_ ;
int nrows = prob->nrows_ ;
double *clo = prob->clo_ ;
double *cup = prob->cup_ ;
double *sol = prob->sol_ ;
double *rlo = prob->rlo_ ;
double *rup = prob->rup_ ;
// If all coefficients positive do more simply
bool allPositive=true;
double * rhs = new double[nrows];
memcpy(rhs,rup,nrows*sizeof(double));
/*
Scan the columns for candidates, and write the indices into sort. We're not
interested in columns that are empty, prohibited, or integral.
Question: Should we exclude singletons, which are useful in other transforms?
Question: Why are we excluding integral columns?
*/
// allow integral columns if asked for
bool allowIntegers = ( prob->presolveOptions_&1)!=0;
int *sort = new int[ncols] ;
int nlook = 0 ;
for (int j = 0 ; j < ncols ; j++) {
if (hincol[j] == 0) continue ;
// check all positive and adjust rhs
if (allPositive) {
double lower = clo[j];
if (lower<cup[j]) {
for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
double value=colels[k];
if (value<0.0)
allPositive=false;
else
rhs[hrow[k]] -= lower*value;
}
} else {
for (int k=mcstrt[j];k<mcstrt[j]+hincol[j];k++) {
double value=colels[k];
rhs[hrow[k]] -= lower*value;
}
}
}
if (prob->colProhibited2(j)) continue ;
if (prob->isInteger(j)&&!allowIntegers) continue ;
sort[nlook++] = j ; }
if (nlook == 0)
{ delete[] sort ;
delete [] rhs;
return (next) ; }
/*
Prep: add the coefficients of each candidate column. To reduce false
positives, multiply each row by a `random' multiplier when forming the
sums. On return from compute_sums, sort and colsum are loaded with the
indices and column sums, respectively, of candidate columns. The pair of
arrays are then sorted by sum so that equal sums are adjacent.
*/
double *colsum = new double[ncols] ;
double *rowmul = new double[nrows] ;
init_random_vec(rowmul,nrows) ;
compute_sums(ncols,hincol,mcstrt,hrow,colels,rowmul,sort,colsum,nlook) ;
CoinSort_2(colsum,colsum+nlook,sort) ;
/*
General prep --- unpack the various vectors we'll need, and allocate arrays
to record the results.
*/
presolvehlink *clink = prob->clink_ ;
double *rowels = prob->rowels_ ;
int *hcol = prob->hcol_ ;
const CoinBigIndex *mrstrt = prob->mrstrt_ ;
int *hinrow = prob->hinrow_ ;
double *dcost = prob->cost_ ;
action *actions = new action [nlook] ;
int nactions = 0 ;
# ifdef ZEROFAULT
memset(actions,0,nlook*sizeof(action)) ;
# endif
int *fixed_down = new int[nlook] ;
int nfixed_down = 0 ;
int *fixed_up = new int[nlook] ;
int nfixed_up = 0 ;
#if 0
// Excluded in the original routine. I'm guessing it's excluded because
// it's just not cost effective to worry about this. -- lh, 040908 --
// It may be the case that several columns are duplicate.
// If not all have the same cost, then we have to make sure
// that we set the most expensive one to its minimum
// now sort in each class by cost
{
double dval = colsum[0] ;
int first = 0 ;
for (int jj = 1; jj < nlook; jj++) {
while (colsum[jj]==dval)
jj++ ;
if (first + 1 < jj) {
double buf[jj - first] ;
for (int i=first; i<jj; ++i)
buf[i-first] = dcost[sort[i]]*maxmin ;
CoinSort_2(buf,buf+jj-first,sort+first) ;
//ekk_sortonDouble(buf,&sort[first],jj-first) ;
}
}
}
#endif
/*
Original comment: It appears to be the case that this loop is finished,
there may still be duplicate cols left. I haven't done anything
about that yet.
Open the main loop to compare column pairs. We'll compare sort[jj] to
sort[tgt]. This allows us to accumulate multiple columns into one. But
we don't manage all-pairs comparison when we can't combine columns.
We can quickly dismiss pairs which have unequal sums or lengths.
*/
int isorted = -1 ;
int tgt = 0 ;
for (int jj = 1 ; jj < nlook ; jj++)
{ if (colsum[jj] != colsum[jj-1]) {
tgt = jj; // Must update before continuing
continue ;
}
int j2 = sort[jj] ;
int j1 = sort[tgt] ;
int len2 = hincol[j2] ;
int len1 = hincol[j1] ;
if (len2 != len1)
{ tgt = jj ;
continue ; }
/*
The final test: sort the columns by row index and compare each index and
coefficient.
*/
CoinBigIndex kcs = mcstrt[j2] ;
CoinBigIndex kce = kcs+hincol[j2] ;
int ishift = mcstrt[j1]-kcs ;
if (len1 > 1 && isorted < j1)
{ CoinSort_2(hrow+mcstrt[j1],hrow+mcstrt[j1]+len1,
colels+mcstrt[j1]) ;
isorted = j1 ; }
if (len2 > 1 && isorted < j2)
{ CoinSort_2(hrow+kcs,hrow+kcs+len2,colels+kcs) ;
isorted = j2 ; }
CoinBigIndex k ;
for (k = kcs ; k < kce ; k++)
{ if (hrow[k] != hrow[k+ishift] || colels[k] != colels[k+ishift])
{ break ; } }
if (k != kce)
{ tgt = jj ;
continue ; }
/*
These really are duplicate columns. Grab values for convenient reference.
Convert the objective coefficients for minimization.
*/
double clo1 = clo[j1] ;
double cup1 = cup[j1] ;
double clo2 = clo[j2] ;
double cup2 = cup[j2] ;
double c1 = dcost[j1]*maxmin ;
double c2 = dcost[j2]*maxmin ;
PRESOLVEASSERT(!(clo1 == cup1 || clo2 == cup2)) ;
// Get reasonable bounds on sum of two variables
double lowerBound=-COIN_DBL_MAX;
double upperBound=COIN_DBL_MAX;
// For now only if lower bounds are zero
if (!clo1&&!clo2) {
if (!allPositive) {
for (k=kcs;k<kce;k++) {
int iRow = hrow[k];
bool posinf = false;
bool neginf = false;
double maxup = 0.0;
double maxdown = 0.0;
// compute sum of all bounds except for j1,j2
CoinBigIndex kk;
CoinBigIndex kre = mrstrt[iRow]+hinrow[iRow];
double value1=0.0;
for (kk=mrstrt[iRow]; kk<kre; kk++) {
int col = hcol[kk];
if (col == j1||col==j2) {
value1=rowels[kk];
continue;
}
double coeff = rowels[kk];
double lb = clo[col];
double ub = cup[col];
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 (kk==kre) {
assert (value1);
if (value1>1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
upperBound = CoinMin(upperBound,(rup[iRow]-maxdown)/value1);
if (!posinf&&rlo[iRow]>-1.0e10)
lowerBound = CoinMax(lowerBound,(rlo[iRow]-maxup)/value1);
} else if (value1<-1.0e-5) {
if (!neginf&&rup[iRow]<1.0e10)
lowerBound = CoinMax(lowerBound,(rup[iRow]-maxdown)/value1);
if (!posinf&&rlo[iRow]>-1.0e10)
upperBound = CoinMin(upperBound,(rlo[iRow]-maxup)/value1);
}
}
}
} else {
// can do faster
lowerBound=0.0;
for (k=kcs;k<kce;k++) {
int iRow = hrow[k];
double value=colels[k];
upperBound = CoinMin(upperBound,rhs[iRow]/value);
}
}
} else {
// Not sure what to do so give up
continue;
}
/*
There are two main cases: The objective coefficients are equal or unequal.
For equal objective coefficients c1 == c2, we can combine the columns by
making the substitution x<j1> = x'<j1> - x<j2>. This will eliminate column
sort[jj] = j2 and leave the new variable x' in column sort[tgt] = j1. tgt
doesn't move.
*/
if (c1 == c2)
{
/*
As far as the presolved lp, there's no difference between columns. But we
need this relation to hold in order to guarantee that we can split the
value of the combined column during postsolve without damaging the basis.
(The relevant case is when the combined column is basic --- we need to be
able to retain column j1 in the basis and make column j2 nonbasic.)
*/
if (!(clo2+cup1 <= clo1+cup2))
{ CoinSwap(j1,j2) ;
CoinSwap(clo1,clo2) ;
CoinSwap(cup1,cup2) ;
tgt = jj ; }
/*
Create the postsolve action before we start to modify the columns.
*/
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d += %d\n",j1,j2,j1,j2)) ;
action *s = &actions[nactions++] ;
s->thislo = clo[j2] ;
s->thisup = cup[j2] ;
s->lastlo = clo[j1] ;
s->lastup = cup[j1] ;
s->ithis = j2 ;
s->ilast = j1 ;
s->nincol = hincol[j2] ;
s->colels = presolve_dupmajor(colels,hrow,hincol[j2],mcstrt[j2]) ;
/*
Combine the columns into column j1. Upper and lower bounds and solution
simply add, and the coefficients are unchanged.
I'm skeptical of pushing a bound to infinity like this, but leave it for now.
-- lh, 040908 --
*/
clo1 += clo2 ;
if (clo1 < -1.0e20)
{ clo1 = -PRESOLVE_INF ; }
clo[j1] = clo1 ;
cup1 += cup2 ;
if (cup1 > 1.0e20)
{ cup1 = PRESOLVE_INF ; }
cup[j1] = cup1 ;
if (sol)
{ sol[j1] += sol[j2] ; }
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic); } }
/*
Empty column j2.
*/
dcost[j2] = 0.0 ;
if (sol)
{ sol[j2] = clo2 ; }
CoinBigIndex k2cs = mcstrt[j2] ;
CoinBigIndex k2ce = k2cs + hincol[j2] ;
for (CoinBigIndex k = k2cs ; k < k2ce ; ++k)
{ presolve_delete_from_row(hrow[k],j2,mrstrt,hinrow,hcol,rowels) ; }
hincol[j2] = 0 ;
PRESOLVE_REMOVE_LINK(clink,j2) ;
continue ; }
/*
Unequal reduced costs. In this case, we may be able to fix one of the columns
or prove dual infeasibility. Given column a_k, duals y, objective
coefficient c_k, the reduced cost cbar_k = c_k - dot(y,a_k). Given
a_1 = a_2, substitute for dot(y,a_1) in the relation for cbar_2 to get
cbar_2 = (c_2 - c_1) + cbar_1
Independent elements here are variable bounds l_k, u_k, and difference
(c_2 - c_1). Infinite bounds for l_k, u_k will constrain the sign of cbar_k.
Assume minimization. If you do the case analysis, you find these cases of
interest:
l_1 u_1 l_2 u_2 cbar_1 c_2-c_1 cbar_2 result
A any finite -inf any <= 0 > 0 <= 0 x_1 -> NBUB
B -inf any any finite <= 0 < 0 < 0 x_2 -> NBUB
C finite any any +inf >= 0 < 0 >= 0 x_1 -> NBLB
D any +inf finite any >= 0 > 0 >= 0 x_2 -> NBLB
E -inf any any +inf <= 0 < 0 >= 0 dual infeas
F any inf -inf any >= 0 > 0 <= 0 dual infeas
G any finite finite any > 0 no inference
H finite any any finite < 0 no inference
The cases labelled dual infeasible are primal unbounded.
To keep the code compact, we'll always aim to take x_2 to bound. In the cases
where x_1 should go to bound, we'll swap. The implementation is boolean
algebra. Define bits for infinite bounds and (c_2 > c_1), then look for the
correct patterns.
*/
else
{ int minterm = 0 ;
bool swapped = false ;
#if PRESOLVE_DEBUG
printf("bounds %g %g\n",lowerBound,upperBound);
#endif
if (c2 > c1) minterm |= 1<<0 ;
if (cup2 >= PRESOLVE_INF) minterm |= 1<<1 ;
if (clo2 <= -PRESOLVE_INF) minterm |= 1<<2 ;
if (cup1 >= PRESOLVE_INF) minterm |= 1<<3 ;
if (clo1 <= -PRESOLVE_INF) minterm |= 1<<4 ;
// for now be careful - just one special case
if (!clo1&&!clo2) {
if (c2 > c1 && cup1 >= upperBound)
minterm |= 1<<3;
else if (c2 < c1 && cup2 >= upperBound)
minterm |= 1<<1;
}
/*
The most common case in a well-formed system should be no inference. We're
looking for x00x1 (case G) and 0xx00 (case H). This is where we have the
potential to miss inferences: If there are three or more columns with the
same sum, sort[tgt] == j1 will only be compared to the second in the
group.
*/
if ((minterm&0x0d) == 0x1 || (minterm&0x13) == 0)
{ tgt = jj ;
continue ; }
/*
Next remove the unbounded cases, 1xx10 and x11x1.
*/
if ((minterm&0x13) == 0x12 || (minterm&0x0d) == 0x0d)
{ prob->setStatus(2) ;
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) Unbounded\n",j1,j2)) ;
break ; }
/*
With the no inference and unbounded cases removed, all that's left are the
cases where we can push a variable to bound. Swap if necessary (x01x1 or
0xx10) so that we're always fixing index j2. This means that column
sort[tgt] = j1 will be fixed. Unswapped, we fix column sort[jj] = j2.
*/
if ((minterm&0x0d) == 0x05 || (minterm&0x13) == 0x02)
{ CoinSwap(j1, j2) ;
CoinSwap(clo1, clo2) ;
CoinSwap(cup1, cup2) ;
CoinSwap(c1, c2) ;
int tmp1 = minterm&0x18 ;
int tmp2 = minterm&0x06 ;
int tmp3 = minterm&0x01 ;
minterm = (tmp1>>2)|(tmp2<<2)|(tmp3^0x01) ;
swapped = true ; }
/*
Force x_2 to upper bound? (Case B, boolean 1X100, where X == don't care.)
*/
if ((minterm&0x13) == 0x10)
{ fixed_up[nfixed_up++] = j2 ;
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBUB\n",j1,j2,j2)) ;
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atUpperBound) ; }
if (sol)
{ double delta2 = cup2-sol[j2] ;
sol[j2] = cup2 ;
sol[j1] -= delta2 ; }
if (swapped)
{ tgt = jj ; }
continue ; }
/*
Force x_2 to lower bound? (Case C, boolean X1011.)
*/
if ((minterm&0x0d) == 0x09)
{ fixed_down[nfixed_down++] = j2 ;
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) %d -> NBLB\n",j1,j2,j2)) ;
if (prob->colstat_)
{ if (prob->getColumnStatus(j1) == CoinPrePostsolveMatrix::basic ||
prob->getColumnStatus(j2) == CoinPrePostsolveMatrix::basic)
{ prob->setColumnStatus(j1,CoinPrePostsolveMatrix::basic) ; }
prob->setColumnStatus(j2,CoinPrePostsolveMatrix::atLowerBound) ; }
if (sol)
{ double delta2 = clo2-sol[j2] ;
sol[j2] = clo2 ;
sol[j1] -= delta2 ; }
if (swapped)
{ tgt = jj ; }
continue ; } }
/*
We should never reach this point in the loop --- all cases force a new
iteration or loop termination. If we get here, something happened that we
didn't anticipate.
*/
PRESOLVE_STMT(printf("DUPCOL: (%d,%d) UNEXPECTED!\n",j1,j2)) ; }
/*
What's left? Deallocate vectors, and call make_fixed_action to handle any
variables that were fixed to bound.
*/
delete[] rowmul ;
delete[] colsum ;
delete[] sort ;
delete [] rhs;
# if PRESOLVE_SUMMARY || PRESOLVE_DEBUG
if (nactions+nfixed_down+nfixed_up > 0)
{ printf("DUPLICATE COLS: %d combined, %d lb, %d ub\n",
nactions,nfixed_down,nfixed_up) ; }
# endif
if (nactions)
{ next = new dupcol_action(nactions,CoinCopyOfArray(actions,nactions),next) ; }
deleteAction(actions,action*) ;
if (nfixed_down)
{ next =
make_fixed_action::presolve(prob,fixed_down,nfixed_down,true,next) ; }
delete[]fixed_down ;
if (nfixed_up)
{ next =
make_fixed_action::presolve(prob,fixed_up,nfixed_up,false,next) ; }
delete[]fixed_up ;
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveDupcol(128) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next) ; }
void dupcol_action::postsolve(CoinPostsolveMatrix *prob) const
{
const action *const actions = actions_;
const int nactions = nactions_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *sol = prob->sol_;
double *dcost = prob->cost_;
double *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *link = prob->link_;
double *rcosts = prob->rcosts_;
double tolerance = prob->ztolzb_;
for (const action *f = &actions[nactions-1]; actions<=f; f--) {
int icol = f->ithis; // was fixed
int icol2 = f->ilast; // was kept
dcost[icol] = dcost[icol2];
clo[icol] = f->thislo;
cup[icol] = f->thisup;
clo[icol2] = f->lastlo;
cup[icol2] = f->lastup;
create_col(icol,f->nincol,f->colels,mcstrt,colels,hrow,link,
&prob->free_list_) ;
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
// hincol[icol] = hincol[icol2]; // right? - no - has to match number in create_col
hincol[icol] = f->nincol;
double l_j = f->thislo;
double u_j = f->thisup;
double l_k = f->lastlo;
double u_k = f->lastup;
double x_k_sol = sol[icol2];
if (l_j>-PRESOLVE_INF&& x_k_sol-l_j>=l_k-tolerance&&x_k_sol-l_j<=u_k+tolerance) {
// j at lb, leave k
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atLowerBound);
sol[icol] = l_j;
sol[icol2] = x_k_sol - sol[icol];
} else if (u_j<PRESOLVE_INF&& x_k_sol-u_j>=l_k-tolerance&&x_k_sol-u_j<=u_k+tolerance) {
// j at ub, leave k
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::atUpperBound);
sol[icol] = u_j;
sol[icol2] = x_k_sol - sol[icol];
} else if (l_k>-PRESOLVE_INF&& x_k_sol-l_k>=l_j-tolerance&&x_k_sol-l_k<=u_j+tolerance) {
// k at lb make j basic
prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
sol[icol2] = l_k;
sol[icol] = x_k_sol - l_k;
prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atLowerBound);
} else if (u_k<PRESOLVE_INF&& x_k_sol-u_k>=l_j-tolerance&&x_k_sol-u_k<=u_j+tolerance) {
// k at ub make j basic
prob->setColumnStatus(icol,prob->getColumnStatus(icol2));
sol[icol2] = u_k;
sol[icol] = x_k_sol - u_k;
prob->setColumnStatus(icol2,CoinPrePostsolveMatrix::atUpperBound);
} else {
// both free! superbasic time
sol[icol] = 0.0; // doesn't matter
prob->setColumnStatus(icol,CoinPrePostsolveMatrix::isFree);
}
// row activity doesn't change
// dj of both variables is the same
rcosts[icol] = rcosts[icol2];
// leave until desctructor
// deleteAction(f->colels,double *);
# if PRESOLVE_DEBUG
const double ztolzb = prob->ztolzb_;
if (! (clo[icol] - ztolzb <= sol[icol] && sol[icol] <= cup[icol] + ztolzb))
printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol], sol[icol], cup[icol]);
if (! (clo[icol2] - ztolzb <= sol[icol2] && sol[icol2] <= cup[icol2] + ztolzb))
printf("BAD DUPCOL BOUNDS: %g %g %g\n", clo[icol2], sol[icol2], cup[icol2]);
# endif
}
// leave until desctructor
// deleteAction(actions_,action *);
}
dupcol_action::~dupcol_action()
{
for (int i = nactions_-1; i >= 0; --i) {
deleteAction(actions_[i].colels, double *);
}
deleteAction(actions_, action*);
}
/*
Routines for duplicate rows. This is definitely unfinished --- there's no
postsolve action.
*/
const char *duprow_action::name () const
{
return ("duprow_action");
}
// This is just ekkredc4, adapted into the new framework.
/*
I've made minimal changes for compatibility with dupcol: An initial scan to
accumulate rows of interest in sort.
-- lh, 040909 --
*/
const CoinPresolveAction
*duprow_action::presolve (CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
double startTime = 0.0;
int startEmptyRows=0;
int startEmptyColumns = 0;
if (prob->tuning_) {
startTime = CoinCpuTime();
startEmptyRows = prob->countEmptyRows();
startEmptyColumns = prob->countEmptyCols();
}
double *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
int ncols = prob->ncols_;
int nrows = prob->nrows_;
/*
Scan the rows for candidates, and write the indices into sort. We're not
interested in rows that are empty or prohibited.
Question: Should we exclude singletons, which are useful in other transforms?
Question: Why are we excluding integral columns?
*/
int *sort = new int[nrows] ;
int nlook = 0 ;
for (int i = 0 ; i < nrows ; i++)
{ if (hinrow[i] == 0) continue ;
if (prob->rowProhibited2(i)) continue ;
sort[nlook++] = i ; }
if (nlook == 0)
{ delete[] sort ;
return (next) ; }
double * workcol = new double[ncols+1];
double * workrow = new double[nrows+1];
init_random_vec(workcol, ncols);
compute_sums(nrows,hinrow,mrstrt,hcol,rowels,workcol,sort,workrow,nlook);
CoinSort_2(workrow,workrow+nlook,sort);
double *rlo = prob->rlo_;
double *rup = prob->rup_;
int nuseless_rows = 0;
bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
double tolerance = prob->feasibilityTolerance_;
double dval = workrow[0];
for (int jj = 1; jj < nlook; jj++) {
if (workrow[jj]==dval) {
int ithis=sort[jj];
int ilast=sort[jj-1];
CoinBigIndex krs = mrstrt[ithis];
CoinBigIndex kre = krs + hinrow[ithis];
if (hinrow[ithis] == hinrow[ilast]) {
int ishift = mrstrt[ilast] - krs;
CoinBigIndex k;
for (k=krs;k<kre;k++) {
if (hcol[k] != hcol[k+ishift] ||
rowels[k] != rowels[k+ishift]) {
break;
}
}
if (k == kre) {
/* now check rhs to see what is what */
double rlo1=rlo[ilast];
double rup1=rup[ilast];
double rlo2=rlo[ithis];
double rup2=rup[ithis];
int idelete=-1;
if (rlo1<=rlo2) {
if (rup2<=rup1) {
/* this is strictly tighter than last */
idelete=ilast;
} else if (fabs(rlo1-rlo2)<1.0e-12) {
/* last is strictly tighter than this */
idelete=ithis;
// swap so can carry on deleting
sort[jj-1]=ithis;
sort[jj]=ilast;
} else {
/* overlapping - could merge */
# if PRESOLVE_DEBUG
printf("overlapping duplicate row %g %g, %g %g\n",
rlo1,rup1,rlo2,rup2);
# endif
if (rup1<rlo2-tolerance&&!fixInfeasibility) {
// infeasible
prob->status_|= 1;
// wrong message - correct if works
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<ithis
<<rlo[ithis]
<<rup[ithis]
<<CoinMessageEol;
break;
}
}
} else {
// rlo2<=rlo1
if (rup1<=rup2) {
/* last is strictly tighter than this */
idelete=ithis;
// swap so can carry on deleting
sort[jj-1]=ithis;
sort[jj]=ilast;
} else {
/* overlapping - could merge */
# if PRESOLVE_DEBUG
printf("overlapping duplicate row %g %g, %g %g\n",
rlo1,rup1,rlo2,rup2);
# endif
if (rup2<rlo1-tolerance&&!fixInfeasibility) {
// infeasible
prob->status_|= 1;
// wrong message - correct if works
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<ithis
<<rlo[ithis]
<<rup[ithis]
<<CoinMessageEol;
break;
}
}
}
if (idelete>=0)
sort[nuseless_rows++]=idelete;
}
}
}
dval=workrow[jj];
}
delete[]workrow;
delete[]workcol;
if (nuseless_rows) {
# if PRESOLVE_SUMMARY
printf("DUPLICATE ROWS: %d\n", nuseless_rows);
# endif
next = useless_constraint_action::presolve(prob,
sort, nuseless_rows,
next);
}
delete[]sort;
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveDuprow(256) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next);
}
void duprow_action::postsolve(CoinPostsolveMatrix *prob) const
{
printf("STILL NO POSTSOLVE FOR DUPROW!\n");
abort();
}
syntax highlighted by Code2HTML, v. 0.9.1