// Copyright (C) 2003, International Business Machines
// Corporation and others. All Rights Reserved.
#include <stdio.h>
#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveEmpty.hpp" // for DROP_COL/DROP_ROW
#include "CoinPresolveZeros.hpp"
#include "CoinPresolveFixed.hpp"
#include "CoinPresolveTripleton.hpp"
#include "CoinPresolvePsdebug.hpp"
#include "CoinMessage.hpp"
#if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
#include "CoinPresolvePsdebug.hpp"
#endif
/*
* Substituting y away:
*
* y = (c - a x - d z) / b
*
* so adjust bounds by: c/b
* and x by: -a/b
* and z by: -d/b
*
* This affects both the row and col representations.
*
* mcstrt only modified if the column must be moved.
*
* for every row in icoly
* if icolx is also has an entry for row
* modify the icolx entry for row
* drop the icoly entry from row and modify the icolx entry
* else
* add a new entry to icolx column
* create a new icolx entry
* (this may require moving the column in memory)
* replace icoly entry from row and replace with icolx entry
*
* same for icolz
* The row and column reps are inconsistent during the routine,
* because icolx in the column rep is updated, and the entries corresponding
* to icolx in the row rep are updated, but nothing concerning icoly
* in the col rep is changed. icoly entries in the row rep are deleted,
* and icolx entries in both reps are consistent.
* At the end, we set the length of icoly to be zero, so the reps would
* be consistent if the row were deleted from the row rep.
* Both the row and icoly must be removed from both reps.
* In the col rep, icoly will be eliminated entirely, at the end of the routine;
* irow occurs in just two columns, one of which (icoly) is eliminated
* entirely, the other is icolx, which is not deleted here.
* In the row rep, irow will be eliminated entirely, but not here;
* icoly is removed from the rows it occurs in.
*/
static bool elim_tripleton(const char *msg,
CoinBigIndex *mcstrt,
double *rlo, double * acts, double *rup,
double *colels,
int *hrow, int *hcol,
int *hinrow, int *hincol,
presolvehlink *clink, int ncols,
presolvehlink *rlink, int nrows,
CoinBigIndex *mrstrt, double *rowels,
//double a, double b, double c,
double coeff_factorx,double coeff_factorz,
double bounds_factor,
int row0, int icolx, int icoly, int icolz)
{
CoinBigIndex kcs = mcstrt[icoly];
CoinBigIndex kce = kcs + hincol[icoly];
CoinBigIndex kcsx = mcstrt[icolx];
CoinBigIndex kcex = kcsx + hincol[icolx];
CoinBigIndex kcsz = mcstrt[icolz];
CoinBigIndex kcez = kcsz + hincol[icolz];
# if PRESOLVE_DEBUG
printf("%s %d x=%d y=%d z=%d cfx=%g cfz=%g nx=%d yrows=(", msg,
row0,icolx,icoly,icolz,coeff_factorx,coeff_factorz,hincol[icolx]) ;
# endif
for (CoinBigIndex kcoly=kcs; kcoly<kce; kcoly++) {
int row = hrow[kcoly];
// even though these values are updated, they remain consistent
PRESOLVEASSERT(kcex == kcsx + hincol[icolx]);
PRESOLVEASSERT(kcez == kcsz + hincol[icolz]);
// we don't need to update the row being eliminated
if (row != row0/* && hinrow[row] > 0*/) {
if (bounds_factor != 0.0) {
// (1)
if (-PRESOLVE_INF < rlo[row])
rlo[row] -= colels[kcoly] * bounds_factor;
// (2)
if (rup[row] < PRESOLVE_INF)
rup[row] -= colels[kcoly] * bounds_factor;
// and solution
if (acts)
{ acts[row] -= colels[kcoly] * bounds_factor; }
}
// see if row appears in colx
CoinBigIndex kcolx = presolve_find_row1(row, kcsx, kcex, hrow);
# if PRESOLVE_DEBUG
printf("%d%s ",row,(kcolx<kcex)?"x+":"") ;
# endif
// see if row appears in colz
CoinBigIndex kcolz = presolve_find_row1(row, kcsz, kcez, hrow);
# if PRESOLVE_DEBUG
printf("%d%s ",row,(kcolz<kcez)?"x+":"") ;
# endif
if (kcolx>=kcex&&kcolz<kcez) {
// swap
int iTemp;
iTemp=kcolx;
kcolx=kcolz;
kcolz=iTemp;
iTemp=kcsx;
kcsx=kcsz;
kcsz=iTemp;
iTemp=kcex;
kcex=kcez;
kcez=iTemp;
iTemp=icolx;
icolx=icolz;
icolz=iTemp;
double dTemp=coeff_factorx;
coeff_factorx=coeff_factorz;
coeff_factorz=dTemp;
}
if (kcolx<kcex) {
// before: both x and y are in the row
// after: only x is in the row
// so: number of elems in col x unchanged, and num elems in row is one less
// update col rep - just modify coefficent
// column y is deleted as a whole at the end of the loop
colels[kcolx] += colels[kcoly] * coeff_factorx;
// update row rep
// first, copy new value for col x into proper place in rowels
CoinBigIndex k2 = presolve_find_col(icolx, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
rowels[k2] = colels[kcolx];
if (kcolz<kcez) {
// before: both z and y are in the row
// after: only z is in the row
// so: number of elems in col z unchanged, and num elems in row is one less
// update col rep - just modify coefficent
// column y is deleted as a whole at the end of the loop
colels[kcolz] += colels[kcoly] * coeff_factorz;
// update row rep
// first, copy new value for col z into proper place in rowels
CoinBigIndex k2 = presolve_find_col(icolz, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
rowels[k2] = colels[kcolz];
// now delete col y from the row; this changes hinrow[row]
presolve_delete_from_row(row, icoly, mrstrt, hinrow, hcol, rowels);
} else {
// before: only y is in the row
// after: only z is in the row
// so: number of elems in col z is one greater, but num elems in row remains same
// update entry corresponding to icolz in row rep
// by just overwriting the icoly entry
{
CoinBigIndex k2 = presolve_find_col(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
hcol[k2] = icolz;
rowels[k2] = colels[kcoly] * coeff_factorz;
}
{
bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,
clink,ncols,icolz);
if (no_mem)
return (true);
// have to adjust various induction variables
kcoly = mcstrt[icoly] + (kcoly - kcs);
kcs = mcstrt[icoly]; // do this for ease of debugging
kce = mcstrt[icoly] + hincol[icoly];
kcolz = mcstrt[icolz] + (kcolz - kcs); // don't really need to do this
kcsz = mcstrt[icolz];
kcez = mcstrt[icolz] + hincol[icolz];
}
// there is now an unused entry in the memory after the column - use it
// mcstrt[ncols] == penultimate index of arrays hrow/colels
hrow[kcez] = row;
colels[kcez] = colels[kcoly] * coeff_factorz; // y factor is 0.0
hincol[icolz]++, kcez++; // expand the col
}
} else {
// before: only y is in the row
// after: only x and z are in the row
// update entry corresponding to icolx in row rep
// by just overwriting the icoly entry
{
CoinBigIndex k2 = presolve_find_col(icoly, mrstrt[row], mrstrt[row]+hinrow[row], hcol);
hcol[k2] = icolx;
rowels[k2] = colels[kcoly] * coeff_factorx;
}
presolve_expand_row(mrstrt,rowels,hcol,hinrow,rlink,nrows,row) ;
// there is now an unused entry in the memory after the column - use it
int krez = mrstrt[row]+hinrow[row];
hcol[krez] = icolz;
rowels[krez] = colels[kcoly] * coeff_factorz;
hinrow[row]++;
{
bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,
clink,ncols,icolx) ;
if (no_mem)
return (true);
// have to adjust various induction variables
kcoly = mcstrt[icoly] + (kcoly - kcs);
kcs = mcstrt[icoly]; // do this for ease of debugging
kce = mcstrt[icoly] + hincol[icoly];
kcolx = mcstrt[icolx] + (kcolx - kcs); // don't really need to do this
kcsx = mcstrt[icolx];
kcex = mcstrt[icolx] + hincol[icolx];
kcolz = mcstrt[icolz] + (kcolz - kcs); // don't really need to do this
kcsz = mcstrt[icolz];
kcez = mcstrt[icolz] + hincol[icolz];
}
// there is now an unused entry in the memory after the column - use it
hrow[kcex] = row;
colels[kcex] = colels[kcoly] * coeff_factorx; // y factor is 0.0
hincol[icolx]++, kcex++; // expand the col
{
bool no_mem = presolve_expand_col(mcstrt,colels,hrow,hincol,clink,
ncols,icolz);
if (no_mem)
return (true);
// have to adjust various induction variables
kcoly = mcstrt[icoly] + (kcoly - kcs);
kcs = mcstrt[icoly]; // do this for ease of debugging
kce = mcstrt[icoly] + hincol[icoly];
kcsx = mcstrt[icolx];
kcex = mcstrt[icolx] + hincol[icolx];
kcsz = mcstrt[icolz];
kcez = mcstrt[icolz] + hincol[icolz];
}
// there is now an unused entry in the memory after the column - use it
hrow[kcez] = row;
colels[kcez] = colels[kcoly] * coeff_factorz; // y factor is 0.0
hincol[icolz]++, kcez++; // expand the col
}
}
}
# if PRESOLVE_DEBUG
printf(")\n") ;
# endif
// delete the whole column
hincol[icoly] = 0;
return (false);
}
/*
*
* The col rep and row rep must be consistent.
*/
const CoinPresolveAction *tripleton_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 *colels = prob->colels_;
int *hrow = prob->hrow_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int ncols = prob->ncols_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *rowels = prob->rowels_;
int *hcol = prob->hcol_;
CoinBigIndex *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
int nrows = prob->nrows_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
presolvehlink *clink = prob->clink_;
presolvehlink *rlink = prob->rlink_;
const unsigned char *integerType = prob->integerType_;
double *cost = prob->cost_;
int numberLook = prob->numberRowsToDo_;
int iLook;
int * look = prob->rowsToDo_;
const double ztolzb = prob->ztolzb_;
action * actions = new action [nrows];
# ifdef ZEROFAULT
// initialise alignment padding bytes
memset(actions,0,nrows*sizeof(action)) ;
# endif
int nactions = 0;
int *zeros = new int[ncols];
memset(zeros,0,ncols*sizeof(int));
int nzeros = 0;
// If rowstat exists then all do
unsigned char *rowstat = prob->rowstat_;
double *acts = prob->acts_;
// unsigned char * colstat = prob->colstat_;
# if PRESOLVE_CONSISTENCY
presolve_links_ok(prob) ;
# endif
// wasfor (int irow=0; irow<nrows; irow++)
for (iLook=0;iLook<numberLook;iLook++) {
int irow = look[iLook];
if (hinrow[irow] == 3 &&
fabs(rup[irow] - rlo[irow]) <= ZTOLDP) {
double rhs = rlo[irow];
CoinBigIndex krs = mrstrt[irow];
CoinBigIndex kre = krs + hinrow[irow];
int icolx, icoly, icolz;
double coeffx, coeffy, coeffz;
CoinBigIndex k;
/* locate first column */
for (k=krs; k<kre; k++) {
if (hincol[hcol[k]] > 0) {
break;
}
}
PRESOLVEASSERT(k<kre);
coeffx = rowels[k];
if (fabs(coeffx) < ZTOLDP)
continue;
icolx = hcol[k];
/* locate second column */
for (k++; k<kre; k++) {
if (hincol[hcol[k]] > 0) {
break;
}
}
PRESOLVEASSERT(k<kre);
coeffy = rowels[k];
if (fabs(coeffy) < ZTOLDP)
continue;
icoly = hcol[k];
/* locate third column */
for (k++; k<kre; k++) {
if (hincol[hcol[k]] > 0) {
break;
}
}
PRESOLVEASSERT(k<kre);
coeffz = rowels[k];
if (fabs(coeffz) < ZTOLDP)
continue;
icolz = hcol[k];
// For now let's do obvious one
if (coeffx*coeffz>0.0) {
if(coeffx*coeffy>0.0)
continue;
} else if (coeffx*coeffy>0.0) {
int iTemp = icoly;
icoly=icolz;
icolz=iTemp;
double dTemp = coeffy;
coeffy=coeffz;
coeffz=dTemp;
} else {
int iTemp = icoly;
icoly=icolx;
icolx=iTemp;
double dTemp = coeffy;
coeffy=coeffx;
coeffx=dTemp;
}
// Not all same sign and y is odd one out
// don't bother with fixed variables
if (!(fabs(cup[icolx] - clo[icolx]) < ZTOLDP) &&
!(fabs(cup[icoly] - clo[icolx]) < ZTOLDP) &&
!(fabs(cup[icolz] - clo[icoly]) < ZTOLDP)) {
assert (coeffx*coeffz>0.0&&coeffx*coeffy<0.0);
/* don't do if y integer for now */
if (integerType[icoly])
continue;
// Only do if does not give implicit bounds on x and z
double cx = - coeffx/coeffy;
double cz = - coeffz/coeffy;
double rhsRatio = rhs/coeffy;
if (clo[icoly]>-1.0e30) {
if (clo[icolx]<-1.0e30||clo[icolz]<-1.0e30)
continue;
if (cx*clo[icolx]+cz*clo[icolz]+rhsRatio<clo[icoly]-ztolzb)
continue;
}
if (cup[icoly]<1.0e30) {
if (cup[icolx]>1.0e30||cup[icolz]>1.0e30)
continue;
if (cx*cup[icolx]+cz*cup[icolz]+rhsRatio>cup[icoly]+ztolzb)
continue;
}
CoinBigIndex krowx,krowy,krowz;
/* find this row in each of the columns and do counts */
bool singleton=false;
for (k=mcstrt[icoly]; k<mcstrt[icoly]+hincol[icoly]; k++) {
int jrow=hrow[k];
if (hinrow[jrow]==1)
singleton=true;
if (jrow == irow)
krowy=k;
else
prob->setRowUsed(jrow);
}
int nDuplicate=0;
for (k=mcstrt[icolx]; k<mcstrt[icolx]+hincol[icolx]; k++) {
int jrow=hrow[k];
if (jrow == irow)
krowx=k;
else if (prob->rowUsed(jrow))
nDuplicate++;;
}
for (k=mcstrt[icolz]; k<mcstrt[icolz]+hincol[icolz]; k++) {
int jrow=hrow[k];
if (jrow == irow)
krowz=k;
else if (prob->rowUsed(jrow))
nDuplicate++;;
}
int nAdded=hincol[icoly]-3-nDuplicate;
for (k=mcstrt[icoly]; k<mcstrt[icoly]+hincol[icoly]; k++) {
int jrow=hrow[k];
prob->unsetRowUsed(jrow);
}
// let singleton rows be taken care of first
if (singleton)
continue;
//if (nAdded<=1)
//printf("%d elements added, hincol %d , dups %d\n",nAdded,hincol[icoly],nDuplicate);
if (nAdded>1)
continue;
// it is possible that both x/z and y are singleton columns
// that can cause problems
if ((hincol[icolx] == 1 ||hincol[icolz] == 1) && hincol[icoly] == 1)
continue;
// common equations are of the form ax + by = 0, or x + y >= lo
{
action *s = &actions[nactions];
nactions++;
s->row = irow;
s->icolx = icolx;
s->icolz = icolz;
s->icoly = icoly;
s->cloy = clo[icoly];
s->cupy = cup[icoly];
s->costy = cost[icoly];
s->rlo = rlo[irow];
s->rup = rup[irow];
s->coeffx = coeffx;
s->coeffy = coeffy;
s->coeffz = coeffz;
s->ncoly = hincol[icoly];
s->colel = presolve_dupmajor(colels, hrow, hincol[icoly],
mcstrt[icoly]);
}
// costs
// the effect of maxmin cancels out
cost[icolx] += cost[icoly] * cx;
cost[icolz] += cost[icoly] * cz;
prob->change_bias(cost[icoly] * rhs / coeffy);
//if (cost[icoly]*rhs)
//printf("change %g col %d cost %g rhs %g coeff %g\n",cost[icoly]*rhs/coeffy,
// icoly,cost[icoly],rhs,coeffy);
if (rowstat) {
// update solution and basis
int numberBasic=0;
if (prob->columnIsBasic(icoly))
numberBasic++;
if (prob->rowIsBasic(irow))
numberBasic++;
if (numberBasic>1) {
if (!prob->columnIsBasic(icolx))
prob->setColumnStatus(icolx,CoinPrePostsolveMatrix::basic);
else
prob->setColumnStatus(icolz,CoinPrePostsolveMatrix::basic);
}
}
// Update next set of actions
{
prob->addCol(icolx);
int i,kcs,kce;
kcs = mcstrt[icoly];
kce = kcs + hincol[icoly];
for (i=kcs;i<kce;i++) {
int row = hrow[i];
prob->addRow(row);
}
kcs = mcstrt[icolx];
kce = kcs + hincol[icolx];
for (i=kcs;i<kce;i++) {
int row = hrow[i];
prob->addRow(row);
}
prob->addCol(icolz);
kcs = mcstrt[icolz];
kce = kcs + hincol[icolz];
for (i=kcs;i<kce;i++) {
int row = hrow[i];
prob->addRow(row);
}
}
/* transfer the colx factors to coly */
bool no_mem = elim_tripleton("ELIMT",
mcstrt, rlo, acts, rup, colels,
hrow, hcol, hinrow, hincol,
clink, ncols, rlink, nrows,
mrstrt, rowels,
cx,
cz,
rhs / coeffy,
irow, icolx, icoly,icolz);
if (no_mem)
throwCoinError("out of memory",
"tripleton_action::presolve");
// now remove irow from icolx and icolz in the col rep
// better if this were first.
presolve_delete_from_col(irow,icolx,mcstrt,hincol,hrow,colels) ;
presolve_delete_from_col(irow,icolz,mcstrt,hincol,hrow,colels) ;
// eliminate irow entirely from the row rep
hinrow[irow] = 0;
// eliminate irow entirely from the row rep
PRESOLVE_REMOVE_LINK(rlink, irow);
// eliminate coly entirely from the col rep
PRESOLVE_REMOVE_LINK(clink, icoly);
cost[icoly] = 0.0;
rlo[irow] = 0.0;
rup[irow] = 0.0;
zeros[icolx] = 1; // check for zeros
zeros[icolz] = 1; // check for zeros
nzeros++;
}
# if PRESOLVE_CONSISTENCY
presolve_links_ok(prob) ;
presolve_consistent(prob);
# endif
}
}
if (nactions) {
# if PRESOLVE_SUMMARY
printf("NTRIPLETONS: %d\n", nactions);
# endif
action *actions1 = new action[nactions];
CoinMemcpyN(actions, nactions, actions1);
next = new tripleton_action(nactions, actions1, next);
if (nzeros) {
int i;
nzeros=0;
for (i=0;i<ncols;i++)
if (zeros[i])
zeros[nzeros++]=i;
next = drop_zero_coefficients_action::presolve(prob, zeros, nzeros, next);
}
}
delete[]zeros;
deleteAction(actions,action*);
if (prob->tuning_) {
double thisTime=CoinCpuTime();
int droppedRows = prob->countEmptyRows() - startEmptyRows ;
int droppedColumns = prob->countEmptyCols() - startEmptyColumns;
printf("CoinPresolveTripleton(8) - %d rows, %d columns dropped in time %g, total %g\n",
droppedRows,droppedColumns,thisTime-startTime,thisTime-prob->startTime_);
}
return (next);
}
void tripleton_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_;
double *clo = prob->clo_;
double *cup = prob->cup_;
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_;
unsigned char *colstat = prob->colstat_;
unsigned char *rowstat = prob->rowstat_;
const double maxmin = prob->maxmin_;
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
char *cdone = prob->cdone_;
char *rdone = prob->rdone_;
# endif
CoinBigIndex &free_list = prob->free_list_;
const double ztolzb = prob->ztolzb_;
const double ztoldj = prob->ztoldj_;
// Space for accumulating two columns
int nrows = prob->nrows_;
int * index1 = new int[nrows];
double * element1 = new double[nrows];
memset(element1,0,nrows*sizeof(double));
int * index2 = new int[nrows];
double * element2 = new double[nrows];
memset(element2,0,nrows*sizeof(double));
for (const action *f = &actions[nactions-1]; actions<=f; f--) {
int irow = f->row;
// probably don't need this
double ylo0 = f->cloy;
double yup0 = f->cupy;
double coeffx = f->coeffx;
double coeffy = f->coeffy;
double coeffz = f->coeffz;
int jcolx = f->icolx;
int jcoly = f->icoly;
int jcolz = f->icolz;
// needed?
double rhs = f->rlo;
/* the column was in the reduced problem */
PRESOLVEASSERT(cdone[jcolx] && rdone[irow]==DROP_ROW&&cdone[jcolz]);
PRESOLVEASSERT(cdone[jcoly]==DROP_COL);
// probably don't need this
rlo[irow] = f->rlo;
rup[irow] = f->rup;
// probably don't need this
clo[jcoly] = ylo0;
cup[jcoly] = yup0;
dcost[jcoly] = f->costy;
dcost[jcolx] += f->costy*coeffx/coeffy;
dcost[jcolz] += f->costy*coeffz/coeffy;
// this is why we want coeffx < coeffy (55)
sol[jcoly] = (rhs - coeffx * sol[jcolx] - coeffz * sol[jcolz]) / coeffy;
// since this row is fixed
acts[irow] = rhs;
// acts[irow] always ok, since slack is fixed
if (rowstat)
prob->setRowStatus(irow,CoinPrePostsolveMatrix::atLowerBound);
// CLAIM:
// if the new pi value is chosen to keep the reduced cost
// of col x at its prior value, then the reduced cost of
// col y will be 0.
// also have to update row activities and bounds for rows affected by jcoly
//
// sol[jcolx] was found for coeffx that
// was += colels[kcoly] * coeff_factor;
// where coeff_factor == -coeffx / coeffy
//
// its contribution to activity was
// (colels[kcolx] + colels[kcoly] * coeff_factor) * sol[jcolx] (1)
//
// After adjustment, the two columns contribute:
// colels[kcoly] * sol[jcoly] + colels[kcolx] * sol[jcolx]
// == colels[kcoly] * ((rhs - coeffx * sol[jcolx]) / coeffy) + colels[kcolx] * sol[jcolx]
// == colels[kcoly] * rhs/coeffy + colels[kcoly] * coeff_factor * sol[jcolx] + colels[kcolx] * sol[jcolx]
// colels[kcoly] * rhs/coeffy + the expression (1)
//
// therefore, we must increase the row bounds by colels[kcoly] * rhs/coeffy,
// which is similar to the bias
double djy = maxmin * dcost[jcoly];
double djx = maxmin * dcost[jcolx];
double djz = maxmin * dcost[jcolz];
double bounds_factor = rhs/coeffy;
// need to reconstruct x and z
double multiplier1 = coeffx/coeffy;
double multiplier2 = coeffz/coeffy;
int * indy = reinterpret_cast<int *>(f->colel+f->ncoly);
int ystart = NO_LINK;
int nX=0,nZ=0;
int i,iRow;
for (i=0; i<f->ncoly; ++i) {
int iRow = indy[i];
double yValue = f->colel[i];
CoinBigIndex k = free_list;
assert(k >= 0 && k < prob->bulk0_) ;
free_list = link[free_list];
if (iRow != irow) {
// are these tests always true???
// undo elim_tripleton(1)
if (-PRESOLVE_INF < rlo[iRow])
rlo[iRow] += yValue * bounds_factor;
// undo elim_tripleton(2)
if (rup[iRow] < PRESOLVE_INF)
rup[iRow] += yValue * bounds_factor;
acts[iRow] += yValue * bounds_factor;
djy -= rowduals[iRow] * yValue;
}
hrow[k] = iRow;
PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
colels[k] = yValue;
link[k] = ystart;
ystart = k;
element1[iRow]=yValue*multiplier1;
index1[nX++]=iRow;
element2[iRow]=yValue*multiplier2;
index2[nZ++]=iRow;
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
mcstrt[jcoly] = ystart;
hincol[jcoly] = f->ncoly;
// find the tail
CoinBigIndex k=mcstrt[jcolx];
CoinBigIndex last = NO_LINK;
int numberInColumn = hincol[jcolx];
int numberToDo=numberInColumn;
for (i=0; i<numberToDo; ++i) {
iRow = hrow[k];
assert (iRow>=0&&iRow<nrows);
double value = colels[k]+element1[iRow];
element1[iRow]=0.0;
if (fabs(value)>=1.0e-15) {
colels[k]=value;
last=k;
k = link[k];
if (iRow != irow)
djx -= rowduals[iRow] * value;
} else {
numberInColumn--;
// add to free list
int nextk = link[k];
link[k]=free_list;
free_list=k;
assert (k>=0);
k=nextk;
if (last!=NO_LINK)
link[last]=k;
else
mcstrt[jcolx]=k;
}
}
for (i=0;i<nX;i++) {
int iRow = index1[i];
double xValue = element1[iRow];
element1[iRow]=0.0;
if (fabs(xValue)>=1.0e-15) {
if (iRow != irow)
djx -= rowduals[iRow] * xValue;
numberInColumn++;
CoinBigIndex k = free_list;
assert(k >= 0 && k < prob->bulk0_) ;
free_list = link[free_list];
hrow[k] = iRow;
PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
colels[k] = xValue;
if (last!=NO_LINK)
link[last]=k;
else
mcstrt[jcolx]=k;
last = k;
}
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
link[last]=NO_LINK;
assert(numberInColumn);
hincol[jcolx] = numberInColumn;
// find the tail
k=mcstrt[jcolz];
last = NO_LINK;
numberInColumn = hincol[jcolz];
numberToDo=numberInColumn;
for (i=0; i<numberToDo; ++i) {
iRow = hrow[k];
assert (iRow>=0&&iRow<nrows);
double value = colels[k]+element2[iRow];
element2[iRow]=0.0;
if (fabs(value)>=1.0e-15) {
colels[k]=value;
last=k;
k = link[k];
if (iRow != irow)
djz -= rowduals[iRow] * value;
} else {
numberInColumn--;
// add to free list
int nextk = link[k];
assert(free_list>=0);
link[k]=free_list;
free_list=k;
assert (k>=0);
k=nextk;
if (last!=NO_LINK)
link[last]=k;
else
mcstrt[jcolz]=k;
}
}
for (i=0;i<nZ;i++) {
int iRow = index2[i];
double zValue = element2[iRow];
element2[iRow]=0.0;
if (fabs(zValue)>=1.0e-15) {
if (iRow != irow)
djz -= rowduals[iRow] * zValue;
numberInColumn++;
CoinBigIndex k = free_list;
assert(k >= 0 && k < prob->bulk0_) ;
free_list = link[free_list];
hrow[k] = iRow;
PRESOLVEASSERT(rdone[hrow[k]] || hrow[k] == irow);
colels[k] = zValue;
if (last!=NO_LINK)
link[last]=k;
else
mcstrt[jcolz]=k;
last = k;
}
}
# if PRESOLVE_CONSISTENCY
presolve_check_free_list(prob) ;
# endif
link[last]=NO_LINK;
assert(numberInColumn);
hincol[jcolz] = numberInColumn;
// The only problem with keeping the reduced costs the way they were
// was that the variable's bound may have moved, requiring it
// to become basic.
//printf("djs x - %g (%g), y - %g (%g)\n",djx,coeffx,djy,coeffy);
if (colstat) {
if (prob->columnIsBasic(jcolx) ||
(fabs(clo[jcolx] - sol[jcolx]) < ztolzb && rcosts[jcolx] >= -ztoldj) ||
(fabs(cup[jcolx] - sol[jcolx]) < ztolzb && rcosts[jcolx] <= ztoldj) ||
(prob->getColumnStatus(jcolx) ==CoinPrePostsolveMatrix::isFree&&
fabs(rcosts[jcolx]) <= ztoldj)) {
// colx or y is fine as it is - make coly basic
prob->setColumnStatus(jcoly,CoinPrePostsolveMatrix::basic);
// 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[irow] = djy / coeffy;
rcosts[jcolx] = djx - rowduals[irow] * coeffx;
# if PRESOLVE_DEBUG
if (prob->columnIsBasic(jcolx)&&fabs(rcosts[jcolx])>1.0e-5)
printf("bad dj %d %g\n",jcolx,rcosts[jcolx]);
# endif
rcosts[jcolz] = djz - rowduals[irow] * coeffz;
//if (prob->columnIsBasic(jcolz))
//assert (fabs(rcosts[jcolz])<1.0e-5);
rcosts[jcoly] = 0.0;
} else {
prob->setColumnStatus(jcolx,CoinPrePostsolveMatrix::basic);
prob->setColumnStatusUsingValue(jcoly);
// change rowduals[jcolx] enough to cancel out rcosts[jcolx]
rowduals[irow] = djx / coeffx;
rcosts[jcolx] = 0.0;
// change rowduals[jcolx] enough to cancel out rcosts[jcolx]
//rowduals[irow] = djz / coeffz;
//rcosts[jcolz] = 0.0;
rcosts[jcolz] = djz - rowduals[irow] * coeffz;
rcosts[jcoly] = djy - rowduals[irow] * coeffy;
}
} else {
// No status array
// 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[irow] = djy / coeffy;
rcosts[jcoly] = 0.0;
}
// DEBUG CHECK
# if PRESOLVE_DEBUG
{
CoinBigIndex k = mcstrt[jcolx];
int nx = hincol[jcolx];
double dj = maxmin * dcost[jcolx];
for (int i=0; i<nx; ++i) {
int row = hrow[k];
double coeff = colels[k];
k = link[k];
dj -= rowduals[row] * coeff;
}
if (! (fabs(rcosts[jcolx] - dj) < 100*ZTOLDP))
printf("BAD DOUBLE X DJ: %d %d %g %g\n",
irow, jcolx, rcosts[jcolx], dj);
rcosts[jcolx]=dj;
}
{
CoinBigIndex k = mcstrt[jcoly];
int ny = hincol[jcoly];
double dj = maxmin * dcost[jcoly];
for (int i=0; i<ny; ++i) {
int row = hrow[k];
double coeff = colels[k];
k = link[k];
dj -= rowduals[row] * coeff;
//printf("b %d coeff %g dual %g dj %g\n",
// row,coeff,rowduals[row],dj);
}
if (! (fabs(rcosts[jcoly] - dj) < 100*ZTOLDP))
printf("BAD DOUBLE Y DJ: %d %d %g %g\n",
irow, jcoly, rcosts[jcoly], dj);
rcosts[jcoly]=dj;
//exit(0);
}
# endif
# if PRESOLVE_DEBUG || PRESOLVE_CONSISTENCY
cdone[jcoly] = TRIPLETON;
rdone[irow] = TRIPLETON;
# endif
}
delete [] index1;
delete [] element1;
delete [] index2;
delete [] element2;
}
tripleton_action::~tripleton_action()
{
for (int i=nactions_-1; i>=0; i--) {
delete[]actions_[i].colel;
}
deleteAction(actions_,action*);
}
static double *tripleton_mult;
static int *tripleton_id;
void check_tripletons(const CoinPresolveAction * paction)
{
const CoinPresolveAction * paction0 = paction;
if (paction) {
check_tripletons(paction->next);
if (strcmp(paction0->name(), "tripleton_action") == 0) {
const tripleton_action *daction = (const tripleton_action *)paction0;
for (int i=daction->nactions_-1; i>=0; --i) {
int icolx = daction->actions_[i].icolx;
int icoly = daction->actions_[i].icoly;
double coeffx = daction->actions_[i].coeffx;
double coeffy = daction->actions_[i].coeffy;
tripleton_mult[icoly] = -coeffx/coeffy;
tripleton_id[icoly] = icolx;
}
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1