// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <stdio.h>
#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "CoinPresolveMatrix.hpp"
#include "CoinPresolveEmpty.hpp"
#include "CoinMessage.hpp"
/* \file
Routines to remove/reinsert empty columns and rows.
*/
/*
Physically remove empty columns, compressing mcstrt and hincol. The major
side effect is that columns are renumbered, thus clink_ is no longer valid
and must be rebuilt.
It's necessary to rebuild clink_ in order to do direct conversion of a
CoinPresolveMatrix to a CoinPostsolveMatrix by transferring the data arrays.
Without clink_, it's impractical to build link_ to match the transferred bulk
storage.
*/
const CoinPresolveAction
*drop_empty_cols_action::presolve (CoinPresolveMatrix *prob,
int *ecols,
int necols,
const CoinPresolveAction *next)
{
int ncols = prob->ncols_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
presolvehlink *clink = prob->clink_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *dcost = prob->cost_;
const double ztoldj = prob->ztoldj_;
unsigned char * integerType = prob->integerType_;
int * originalColumn = prob->originalColumn_;
const double maxmin = prob->maxmin_;
double * sol = prob->sol_;
unsigned char * colstat = prob->colstat_;
action *actions = new action[necols];
int * colmapping = new int [ncols+1];
bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
CoinZeroN(colmapping,ncols);
int i;
for (i=necols-1; i>=0; i--) {
int jcol = ecols[i];
colmapping[jcol]=-1;
action &e = actions[i];
e.jcol = jcol;
e.clo = clo[jcol];
e.cup = cup[jcol];
// adjust if integer
if (integerType[jcol]) {
e.clo = ceil(e.clo-1.0e-9);
e.cup = floor(e.cup+1.0e-9);
if (e.clo>e.cup&&!fixInfeasibility) {
// infeasible
prob->status_|= 1;
prob->messageHandler()->message(COIN_PRESOLVE_COLINFEAS,
prob->messages())
<<jcol
<<e.clo
<<e.cup
<<CoinMessageEol;
}
}
e.cost = dcost[jcol];
// there are no more constraints on this variable,
// so we had better be able to compute the answer now
if (fabs(dcost[jcol])<ztoldj)
dcost[jcol]=0.0;
if (dcost[jcol] * maxmin == 0.0) {
// hopefully, we can make this non-basic
// what does OSL currently do in this case???
e.sol = (-PRESOLVE_INF < e.clo
? e.clo
: e.cup < PRESOLVE_INF
? e.cup
: 0.0);
} else if (dcost[jcol] * maxmin > 0.0) {
if (-PRESOLVE_INF < e.clo)
e.sol = e.clo;
else {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDB,
prob->messages())
<<jcol
<<CoinMessageEol;
prob->status_ |= 2;
break;
}
} else {
if (e.cup < PRESOLVE_INF)
e.sol = e.cup;
else {
prob->messageHandler()->message(COIN_PRESOLVE_COLUMNBOUNDA,
prob->messages())
<<jcol
<<CoinMessageEol;
prob->status_ |= 2;
break;
}
}
#if PRESOLVE_DEBUG
if (e.sol * dcost[jcol]) {
//printf("\a>>>NON-ZERO COST FOR EMPTY COL %d: %g\n", jcol, dcost[jcol]);
}
#endif
prob->change_bias(e.sol * dcost[jcol]);
}
int ncols2=0;
// now move remaining ones down
for (i=0;i<ncols;i++) {
if (!colmapping[i]) {
mcstrt[ncols2] = mcstrt[i];
hincol[ncols2] = hincol[i];
clo[ncols2] = clo[i];
cup[ncols2] = cup[i];
dcost[ncols2] = dcost[i];
if (sol) {
sol[ncols2] = sol[i];
colstat[ncols2] = colstat[i];
}
integerType[ncols2] = integerType[i];
originalColumn[ncols2] = originalColumn[i];
colmapping[i] = ncols2++;
}
}
mcstrt[ncols2] = mcstrt[ncols] ;
colmapping[ncols] = ncols2 ;
/*
Rebuild clink_. At this point, all empty columns are linked out, so the
only columns left are columns that are to be saved, hence available in
colmapping. All we need to do is walk clink_ and write the new entries
into a new array.
*/
{ presolvehlink *newclink = new presolvehlink [ncols2+1] ;
for (int oldj = ncols ; oldj >= 0 ; oldj = clink[oldj].pre)
{ presolvehlink &oldlnk = clink[oldj] ;
int newj = colmapping[oldj] ;
assert(newj >= 0 && newj <= ncols2) ;
presolvehlink &newlnk = newclink[newj] ;
if (oldlnk.suc >= 0)
{ newlnk.suc = colmapping[oldlnk.suc] ; }
else
{ newlnk.suc = NO_LINK ; }
if (oldlnk.pre >= 0)
{ newlnk.pre = colmapping[oldlnk.pre] ; }
else
{ newlnk.pre = NO_LINK ; } }
delete [] clink ;
prob->clink_ = newclink ; }
delete [] colmapping;
prob->ncols_ = ncols2;
return (new drop_empty_cols_action(necols, actions, next));
}
const CoinPresolveAction *drop_empty_cols_action::presolve(CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
const int *hincol = prob->hincol_;
// const double *clo = prob->clo_;
// const double *cup = prob->cup_;
int ncols = prob->ncols_;
int i;
int nempty = 0;
int * empty = new int [ncols];
CoinBigIndex nelems2 = 0 ;
// count empty cols
for (i=0; i<ncols; i++)
{ nelems2 += hincol[i] ;
if (hincol[i] == 0) {
#if PRESOLVE_DEBUG
if (nempty==0)
printf("UNUSED COLS: ");
else
if (i < 100 && nempty%25 == 0)
printf("\n") ;
else
if (i >= 100 && i < 1000 && nempty%19 == 0)
printf("\n") ;
else
if (i >= 1000 && nempty%15 == 0)
printf("\n") ;
printf("%d ", i);
#endif
empty[nempty++] = i;
}
}
prob->nelems_ = nelems2 ;
if (nempty) {
#if PRESOLVE_DEBUG
printf("\ndropped %d cols\n", nempty);
#endif
next = drop_empty_cols_action::presolve(prob, empty, nempty, next);
}
delete [] empty;
return (next);
}
void drop_empty_cols_action::postsolve(CoinPostsolveMatrix *prob) const
{
const int nactions = nactions_;
const action *const actions = actions_;
int ncols = prob->ncols_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
// int *hrow = prob->hrow_;
double *clo = prob->clo_;
double *cup = prob->cup_;
double *sol = prob->sol_;
double *cost = prob->cost_;
double *rcosts = prob->rcosts_;
unsigned char *colstat = prob->colstat_;
const double maxmin = prob->maxmin_;
int ncols2 = ncols+nactions;
int * colmapping = new int [ncols2];
CoinZeroN(colmapping,ncols2);
# if PRESOLVE_DEBUG
char *cdone = prob->cdone_;
# endif
int action_i;
for (action_i = 0; action_i < nactions; action_i++) {
const action *e = &actions[action_i];
int jcol = e->jcol;
colmapping[jcol]=-1;
}
int i;
// now move remaining ones up
for (i=ncols2-1;i>=0;i--) {
if (!colmapping[i]) {
ncols--;
mcstrt[i] = mcstrt[ncols];
hincol[i] = hincol[ncols];
clo[i] = clo[ncols];
cup[i] = cup[ncols];
cost[i] = cost[ncols];
if (sol)
sol[i] = sol[ncols];
if (rcosts)
rcosts[i] = rcosts[ncols];
if (colstat)
colstat[i] = colstat[ncols];
# if PRESOLVE_DEBUG
cdone[i] = cdone[ncols];
# endif
}
}
assert (!ncols);
delete [] colmapping;
for (action_i = 0; action_i < nactions; action_i++) {
const action *e = &actions[action_i];
int jcol = e->jcol;
// now recreate jcol
clo[jcol] = e->clo;
cup[jcol] = e->cup;
if (sol)
sol[jcol] = e->sol;
cost[jcol] = e->cost;
if (rcosts)
rcosts[jcol] = maxmin*cost[jcol];
hincol[jcol] = 0;
mcstrt[jcol] = NO_LINK ;
# if PRESOLVE_DEBUG
cdone[jcol] = DROP_COL;
# endif
if (colstat)
prob->setColumnStatusUsingValue(jcol);
}
prob->ncols_ += nactions;
}
const CoinPresolveAction *drop_empty_rows_action::presolve(CoinPresolveMatrix *prob,
const CoinPresolveAction *next)
{
int ncols = prob->ncols_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *hrow = prob->hrow_;
int nrows = prob->nrows_;
// This is done after row copy needed
//int *mrstrt = prob->mrstrt_;
int *hinrow = prob->hinrow_;
//int *hcol = prob->hcol_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
unsigned char *rowstat = prob->rowstat_;
double *acts = prob->acts_;
int * originalRow = prob->originalRow_;
//presolvehlink *rlink = prob->rlink_;
bool fixInfeasibility = (prob->presolveOptions_&16384)!=0;
int i;
int nactions = 0;
for (i=0; i<nrows; i++)
if (hinrow[i] == 0)
nactions++;
if (nactions == 0)
return next;
else {
action *actions = new action[nactions];
int * rowmapping = new int [nrows];
nactions = 0;
int nrows2=0;
for (i=0; i<nrows; i++) {
if (hinrow[i] == 0) {
action &e = actions[nactions];
# if PRESOLVE_DEBUG
if (nactions == 0)
printf("UNUSED ROWS: ");
else
if (i < 100 && nactions%25 == 0)
printf("\n") ;
else
if (i >= 100 && i < 1000 && nactions%19 == 0)
printf("\n") ;
else
if (i >= 1000 && nactions%15 == 0)
printf("\n") ;
printf("%d ", i);
# endif
nactions++;
if (rlo[i] > 0.0 || rup[i] < 0.0) {
if ((rlo[i]<=prob->feasibilityTolerance_ &&
rup[i]>=-prob->feasibilityTolerance_)||fixInfeasibility) {
rlo[i]=0.0;
rup[i]=0.0;
} else {
prob->status_|= 1;
prob->messageHandler()->message(COIN_PRESOLVE_ROWINFEAS,
prob->messages())
<<i
<<rlo[i]
<<rup[i]
<<CoinMessageEol;
break;
}
}
e.row = i;
e.rlo = rlo[i];
e.rup = rup[i];
rowmapping[i]=-1;
} else {
// move down - we want to preserve order
rlo[nrows2]=rlo[i];
rup[nrows2]=rup[i];
originalRow[nrows2]=i;
if (acts) {
acts[nrows2]=acts[i];
rowstat[nrows2]=rowstat[i];
}
rowmapping[i]=nrows2++;
}
}
// remap matrix
for (i=0;i<ncols;i++) {
int j;
for (j=mcstrt[i];j<mcstrt[i]+hincol[i];j++)
hrow[j] = rowmapping[hrow[j]];
}
delete [] rowmapping;
prob->nrows_ = nrows2;
#if PRESOLVE_DEBUG
if (nactions)
printf("\ndropped %d rows\n", nactions);
#endif
return (new drop_empty_rows_action(nactions, actions, next));
}
}
void drop_empty_rows_action::postsolve(CoinPostsolveMatrix *prob) const
{
const int nactions = nactions_;
const action *const actions = actions_;
int ncols = prob->ncols_;
CoinBigIndex *mcstrt = prob->mcstrt_;
int *hincol = prob->hincol_;
int *hrow = prob->hrow_;
double *rlo = prob->rlo_;
double *rup = prob->rup_;
unsigned char *rowstat = prob->rowstat_;
double *rowduals = prob->rowduals_;
double *acts = prob->acts_;
# if PRESOLVE_DEBUG
char *rdone = prob->rdone_;
# endif
int nrows0 = prob->nrows0_;
int nrows = prob->nrows_;
int * rowmapping = new int [nrows0];
CoinZeroN(rowmapping,nrows0) ;
int i, action_i;
for (action_i = 0; action_i<nactions; action_i++) {
const action *e = &actions[action_i];
int hole = e->row;
rowmapping[hole]=-1;
}
// move data
for (i=nrows0-1; i>=0; i--) {
if (!rowmapping[i]) {
// not a hole
nrows--;
rlo[i]=rlo[nrows];
rup[i]=rup[nrows];
acts[i]=acts[nrows];
rowduals[i]=rowduals[nrows];
if (rowstat)
rowstat[i] = rowstat[nrows];
# if PRESOLVE_DEBUG
rdone[i] = rdone[nrows] ;
# endif
}
}
assert (!nrows);
// set up mapping for matrix
for (i=0;i<nrows0;i++) {
if (!rowmapping[i])
rowmapping[nrows++]=i;
}
for (int j=0; j<ncols; j++) {
CoinBigIndex start = mcstrt[j];
CoinBigIndex end = start + hincol[j];
for (CoinBigIndex k=start; k<end; ++k) {
hrow[k] = rowmapping[hrow[k]];
}
}
delete [] rowmapping;
for (action_i = 0; action_i < nactions; action_i++) {
const action *e = &actions[action_i];
int irow = e->row;
// Now recreate irow
rlo[irow] = e->rlo;
rup[irow] = e->rup;
if (rowstat)
prob->setRowStatus(irow,CoinPrePostsolveMatrix::basic);
rowduals[irow] = 0.0; // ???
acts[irow] = 0.0;
# if PRESOLVE_DEBUG
rdone[irow] = DROP_ROW;
# endif
}
prob->nrows_ = prob->nrows_+nactions;
}
syntax highlighted by Code2HTML, v. 0.9.1