// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <cstdio>
#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "ClpQuadraticObjective.hpp"
#include "ClpNonLinearCost.hpp"
// at end to get min/max!
#include "ClpGubDynamicMatrix.hpp"
#include "ClpMessage.hpp"
//#define CLP_DEBUG
//#define CLP_DEBUG_PRINT
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::ClpGubDynamicMatrix ()
: ClpGubMatrix(),
objectiveOffset_(0.0),
startColumn_(NULL),
row_(NULL),
element_(NULL),
cost_(NULL),
fullStart_(NULL),
id_(NULL),
dynamicStatus_(NULL),
lowerColumn_(NULL),
upperColumn_(NULL),
lowerSet_(NULL),
upperSet_(NULL),
numberGubColumns_(0),
firstAvailable_(0),
savedFirstAvailable_(0),
firstDynamic_(0),
lastDynamic_(0),
numberElements_(0)
{
setType(13);
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::ClpGubDynamicMatrix (const ClpGubDynamicMatrix & rhs)
: ClpGubMatrix(rhs)
{
objectiveOffset_ = rhs.objectiveOffset_;
numberGubColumns_ = rhs.numberGubColumns_;
firstAvailable_ = rhs.firstAvailable_;
savedFirstAvailable_ = rhs.savedFirstAvailable_;
firstDynamic_ = rhs.firstDynamic_;
lastDynamic_ = rhs.lastDynamic_;
numberElements_ = rhs.numberElements_;
startColumn_ = ClpCopyOfArray(rhs.startColumn_,numberGubColumns_+1);
CoinBigIndex numberElements = startColumn_[numberGubColumns_];
row_ = ClpCopyOfArray(rhs.row_,numberElements);;
element_ = ClpCopyOfArray(rhs.element_,numberElements);;
cost_ = ClpCopyOfArray(rhs.cost_,numberGubColumns_);
fullStart_ = ClpCopyOfArray(rhs.fullStart_,numberSets_+1);
id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_,numberGubColumns_);
upperColumn_ = ClpCopyOfArray(rhs.upperColumn_,numberGubColumns_);
dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,numberGubColumns_);
lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
}
/* This is the real constructor*/
ClpGubDynamicMatrix::ClpGubDynamicMatrix(ClpSimplex * model, int numberSets,
int numberGubColumns, const int * starts,
const double * lower, const double * upper,
const CoinBigIndex * startColumn, const int * row,
const double * element, const double * cost,
const double * lowerColumn, const double * upperColumn,
const unsigned char * status)
: ClpGubMatrix()
{
objectiveOffset_ = model->objectiveOffset();
model_ = model;
numberSets_ = numberSets;
numberGubColumns_ = numberGubColumns;
fullStart_ = ClpCopyOfArray(starts,numberSets_+1);
lower_ = ClpCopyOfArray(lower,numberSets_);
upper_ = ClpCopyOfArray(upper,numberSets_);
int numberColumns = model->numberColumns();
int numberRows = model->numberRows();
// Number of columns needed
int numberGubInSmall = numberSets_+numberRows + 2*model->factorizationFrequency()+2;
// for small problems this could be too big
//numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
int numberNeeded = numberGubInSmall + numberColumns;
firstAvailable_ = numberColumns;
savedFirstAvailable_ = numberColumns;
firstDynamic_ = numberColumns;
lastDynamic_ = numberNeeded;
startColumn_ = ClpCopyOfArray(startColumn,numberGubColumns_+1);
CoinBigIndex numberElements = startColumn_[numberGubColumns_];
row_ = ClpCopyOfArray(row,numberElements);
element_ = new double[numberElements];
CoinBigIndex i;
for (i=0;i<numberElements;i++)
element_[i]=element[i];
cost_ = new double[numberGubColumns_];
for (i=0;i<numberGubColumns_;i++) {
cost_[i]=cost[i];
// need sorted
CoinSort_2(row_+startColumn_[i],row_+startColumn_[i+1],element_+startColumn_[i]);
}
if (lowerColumn) {
lowerColumn_ = new double[numberGubColumns_];
for (i=0;i<numberGubColumns_;i++)
lowerColumn_[i]=lowerColumn[i];
} else {
lowerColumn_=NULL;
}
if (upperColumn) {
upperColumn_ = new double[numberGubColumns_];
for (i=0;i<numberGubColumns_;i++)
upperColumn_[i]=upperColumn[i];
} else {
upperColumn_=NULL;
}
if (upperColumn||lowerColumn) {
lowerSet_ = new double[numberSets_];
for (i=0;i<numberSets_;i++) {
if (lower[i]>-1.0e20)
lowerSet_[i]=lower[i];
else
lowerSet_[i] = -1.0e30;
}
upperSet_ = new double[numberSets_];
for (i=0;i<numberSets_;i++) {
if (upper[i]<1.0e20)
upperSet_[i]=upper[i];
else
upperSet_[i] = 1.0e30;
}
} else {
lowerSet_=NULL;
upperSet_=NULL;
}
start_=NULL;
end_=NULL;
dynamicStatus_=NULL;
id_ = new int[numberGubInSmall];
for (i=0;i<numberGubInSmall;i++)
id_[i]=-1;
ClpPackedMatrix* originalMatrixA =
dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
assert (originalMatrixA);
CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
originalMatrixA->setMatrixNull(); // so can be deleted safely
// guess how much space needed
double guess = originalMatrix->getNumElements()+10;
guess /= (double) numberColumns;
guess *= 2*numberGubColumns_;
numberElements_ = (int) CoinMin(guess,10000000.0);
numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
matrix_ = originalMatrix;
flags_ &= ~1;
// resize model (matrix stays same)
model->resize(numberRows,numberNeeded);
if (upperColumn_) {
// set all upper bounds so we have enough space
double * columnUpper = model->columnUpper();
for(i=firstDynamic_;i<lastDynamic_;i++)
columnUpper[i]=1.0e10;
}
// resize matrix
// extra 1 is so can keep number of elements handy
originalMatrix->reserve(numberNeeded,numberElements_,true);
originalMatrix->reserve(numberNeeded+1,numberElements_,false);
originalMatrix->getMutableVectorStarts()[numberColumns]=originalMatrix->getNumElements();
// redo number of columns
numberColumns = matrix_->getNumCols();
backward_ = new int[numberNeeded];
backToPivotRow_ = new int[numberNeeded];
// We know a bit better
delete [] changeCost_;
changeCost_ = new double [numberRows+numberSets_];
keyVariable_ = new int[numberSets_];
// signal to need new ordering
next_ = NULL;
for (int iColumn=0;iColumn<numberNeeded;iColumn++)
backward_[iColumn]=-1;
firstGub_=firstDynamic_;
lastGub_=lastDynamic_;
if (!lowerColumn_&&!upperColumn_)
gubType_=8;
if (status) {
status_ = ClpCopyOfArray(status,numberSets_);
} else {
status_= new unsigned char [numberSets_];
memset(status_,0,numberSets_);
int i;
for (i=0;i<numberSets_;i++) {
// make slack key
setStatus(i,ClpSimplex::basic);
}
}
saveStatus_= new unsigned char [numberSets_];
memset(saveStatus_,0,numberSets_);
savedKeyVariable_= new int [numberSets_];
memset(savedKeyVariable_,0,numberSets_*sizeof(int));
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpGubDynamicMatrix::~ClpGubDynamicMatrix ()
{
delete [] startColumn_;
delete [] row_;
delete [] element_;
delete [] cost_;
delete [] fullStart_;
delete [] id_;
delete [] dynamicStatus_;
delete [] lowerColumn_;
delete [] upperColumn_;
delete [] lowerSet_;
delete [] upperSet_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpGubDynamicMatrix &
ClpGubDynamicMatrix::operator=(const ClpGubDynamicMatrix& rhs)
{
if (this != &rhs) {
ClpGubMatrix::operator=(rhs);
delete [] startColumn_;
delete [] row_;
delete [] element_;
delete [] cost_;
delete [] fullStart_;
delete [] id_;
delete [] dynamicStatus_;
delete [] lowerColumn_;
delete [] upperColumn_;
delete [] lowerSet_;
delete [] upperSet_;
objectiveOffset_ = rhs.objectiveOffset_;
numberGubColumns_ = rhs.numberGubColumns_;
firstAvailable_ = rhs.firstAvailable_;
savedFirstAvailable_ = rhs.savedFirstAvailable_;
firstDynamic_ = rhs.firstDynamic_;
lastDynamic_ = rhs.lastDynamic_;
numberElements_ = rhs.numberElements_;
startColumn_ = ClpCopyOfArray(rhs.startColumn_,numberGubColumns_+1);
int numberElements = startColumn_[numberGubColumns_];
row_ = ClpCopyOfArray(rhs.row_,numberElements);;
element_ = ClpCopyOfArray(rhs.element_,numberElements);;
cost_ = ClpCopyOfArray(rhs.cost_,numberGubColumns_);
fullStart_ = ClpCopyOfArray(rhs.fullStart_,numberSets_+1);
id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
lowerColumn_ = ClpCopyOfArray(rhs.lowerColumn_,numberGubColumns_);
upperColumn_ = ClpCopyOfArray(rhs.upperColumn_,numberGubColumns_);
dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,numberGubColumns_);
lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
}
return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpGubDynamicMatrix::clone() const
{
return new ClpGubDynamicMatrix(*this);
}
// Partial pricing
void
ClpGubDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
int & bestSequence, int & numberWanted)
{
assert(!model->rowScale());
numberWanted=currentWanted_;
if (!numberSets_) {
// no gub
ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
return;
} else {
// and do some proportion of full set
int startG2 = (int) (startFraction*numberSets_);
int endG2 = (int) (endFraction*numberSets_+0.1);
endG2 = CoinMin(endG2,numberSets_);
//printf("gub price - set start %d end %d\n",
// startG2,endG2);
double tolerance=model->currentDualTolerance();
double * reducedCost = model->djRegion();
const double * duals = model->dualRowSolution();
double * cost = model->costRegion();
double bestDj;
int numberRows = model->numberRows();
int numberColumns = lastDynamic_;
// If nothing found yet can go all the way to end
int endAll = endG2;
if (bestSequence<0&&!startG2)
endAll = numberSets_;
if (bestSequence>=0)
bestDj = fabs(reducedCost[bestSequence]);
else
bestDj=tolerance;
int saveSequence = bestSequence;
double djMod=0.0;
double infeasibilityCost = model->infeasibilityCost();
double bestDjMod=0.0;
//printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
// startG2,endG2,numberWanted);
int bestType=-1;
int bestSet=-1;
const double * element =matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
int * length = matrix_->getMutableVectorLengths();
#if 0
// make sure first available is clean (in case last iteration rejected)
cost[firstAvailable_]=0.0;
length[firstAvailable_]=0;
model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
{
for (int i=firstAvailable_;i<lastDynamic_;i++)
assert(!cost[i]);
}
#endif
#ifdef CLP_DEBUG
{
for (int i=firstDynamic_;i<firstAvailable_;i++) {
assert (getDynamicStatus(id_[i-firstDynamic_])==inSmall);
}
}
#endif
int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
for (int iSet = startG2;iSet<endAll;iSet++) {
if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
// give up
numberWanted=0;
break;
} else if (iSet==endG2&&bestSequence>=0) {
break;
}
CoinBigIndex j;
int iBasic = keyVariable_[iSet];
if (iBasic>=numberColumns) {
djMod = - weight(iSet)*infeasibilityCost;
} else {
// get dj without
assert (model->getStatus(iBasic)==ClpSimplex::basic);
djMod=0.0;
for (j=startColumn[iBasic];
j<startColumn[iBasic]+length[iBasic];j++) {
int jRow=row[j];
djMod -= duals[jRow]*element[j];
}
djMod += cost[iBasic];
// See if gub slack possible - dj is djMod
if (getStatus(iSet)==ClpSimplex::atLowerBound) {
double value = -djMod;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSet)) {
bestDj=value;
bestSequence = numberRows+numberColumns+iSet;
bestDjMod = djMod;
bestType=0;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
abort();
}
}
}
} else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
double value = djMod;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSet)) {
bestDj=value;
bestSequence = numberRows+numberColumns+iSet;
bestDjMod = djMod;
bestType=0;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
abort();
}
}
}
}
}
for (int iSequence=fullStart_[iSet];iSequence<fullStart_[iSet+1];iSequence++) {
DynamicStatus status = getDynamicStatus(iSequence);
if (status!=inSmall) {
double value=cost_[iSequence]-djMod;
for (j=startColumn_[iSequence];
j<startColumn_[iSequence+1];j++) {
int jRow=row_[j];
value -= duals[jRow]*element_[j];
}
// change sign if at lower bound
if (status==atLowerBound)
value = -value;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!flagged(iSequence)) {
bestDj=value;
bestSequence = iSequence;
bestDjMod = djMod;
bestType=1;
bestSet = iSet;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
}
}
if (numberWanted<=0) {
numberWanted=0;
break;
}
}
// Do packed part before gub and small gub - but lightly
int saveMinNeg=minimumGoodReducedCosts_;
int saveSequence2 = bestSequence;
if (bestSequence>=0)
minimumGoodReducedCosts_=-2;
int saveLast=lastGub_;
lastGub_=firstAvailable_;
currentWanted_=numberWanted;
ClpGubMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
minimumGoodReducedCosts_=saveMinNeg;
lastGub_=saveLast;
if (bestSequence!=saveSequence2) {
bestType=-1; // in normal or small gub part
saveSequence=bestSequence;
}
if (bestSequence!=saveSequence||bestType>=0) {
double * lowerColumn = model->lowerRegion();
double * upperColumn = model->upperRegion();
double * solution = model->solutionRegion();
if (bestType>0) {
// recompute dj and create
double value=cost_[bestSequence]-bestDjMod;
for (CoinBigIndex jBigIndex=startColumn_[bestSequence];
jBigIndex<startColumn_[bestSequence+1];jBigIndex++) {
int jRow=row_[jBigIndex];
value -= duals[jRow]*element_[jBigIndex];
}
double * element = matrix_->getMutableElements();
int * row = matrix_->getMutableIndices();
CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
int * length = matrix_->getMutableVectorLengths();
CoinBigIndex numberElements = startColumn[firstAvailable_];
int numberThis = startColumn_[bestSequence+1]-startColumn_[bestSequence];
if (numberElements+numberThis>numberElements_) {
// need to redo
numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
matrix_->reserve(numberColumns,numberElements_);
element = matrix_->getMutableElements();
row = matrix_->getMutableIndices();
// these probably okay but be safe
startColumn = matrix_->getMutableVectorStarts();
length = matrix_->getMutableVectorLengths();
}
// already set startColumn[firstAvailable_]=numberElements;
length[firstAvailable_]=numberThis;
model->costRegion()[firstAvailable_]=cost_[bestSequence];
CoinBigIndex base = startColumn_[bestSequence];
for (int j=0;j<numberThis;j++) {
row[numberElements]=row_[base+j];
element[numberElements++]=element_[base+j];
}
id_[firstAvailable_-firstDynamic_]=bestSequence;
//printf("best %d\n",bestSequence);
backward_[firstAvailable_]=bestSet;
model->solutionRegion()[firstAvailable_]=0.0;
if (!lowerColumn_&&!upperColumn_) {
model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
lowerColumn[firstAvailable_] = 0.0;
upperColumn[firstAvailable_] = COIN_DBL_MAX;
} else {
DynamicStatus status = getDynamicStatus(bestSequence);
if (lowerColumn_)
lowerColumn[firstAvailable_] = lowerColumn_[bestSequence];
else
lowerColumn[firstAvailable_] = 0.0;
if (upperColumn_)
upperColumn[firstAvailable_] = upperColumn_[bestSequence];
else
upperColumn[firstAvailable_] = COIN_DBL_MAX;
if (status==atLowerBound) {
solution[firstAvailable_]=lowerColumn[firstAvailable_];
model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
} else {
solution[firstAvailable_]=upperColumn[firstAvailable_];
model->setStatus(firstAvailable_,ClpSimplex::atUpperBound);
}
}
model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
lowerColumn[firstAvailable_],
upperColumn[firstAvailable_],cost_[bestSequence]);
bestSequence=firstAvailable_;
// firstAvailable_ only updated if good pivot (in updatePivot)
startColumn[firstAvailable_+1]=numberElements;
//printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod);
reducedCost[bestSequence] = value;
gubSlackIn_=-1;
} else {
// slack - make last column
gubSlackIn_= bestSequence - numberRows - numberColumns;
bestSequence = numberColumns + 2*numberRows;
reducedCost[bestSequence] = bestDjMod;
//printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod);
model->setStatus(bestSequence,getStatus(gubSlackIn_));
if (getStatus(gubSlackIn_)==ClpSimplex::atUpperBound)
solution[bestSequence] = upper_[gubSlackIn_];
else
solution[bestSequence] = lower_[gubSlackIn_];
lowerColumn[bestSequence] = lower_[gubSlackIn_];
upperColumn[bestSequence] = upper_[gubSlackIn_];
model->costRegion()[bestSequence] = 0.0;
model->nonLinearCost()->setOne(bestSequence,solution[bestSequence],lowerColumn[bestSequence],
upperColumn[bestSequence],0.0);
}
savedBestSequence_ = bestSequence;
savedBestDj_ = reducedCost[savedBestSequence_];
}
// See if may be finished
if (!startG2&&bestSequence<0)
infeasibilityWeight_=model_->infeasibilityCost();
else if (bestSequence>=0)
infeasibilityWeight_=-1.0;
}
currentWanted_=numberWanted;
}
// This is local to Gub to allow synchronization when status is good
int
ClpGubDynamicMatrix::synchronize(ClpSimplex * model,int mode)
{
int returnNumber=0;
switch (mode) {
case 0:
{
#ifdef CLP_DEBUG
{
for (int i=0;i<numberSets_;i++)
assert(toIndex_[i]==-1);
}
#endif
// lookup array
int * lookup = new int[lastDynamic_];
int iColumn;
int numberColumns = model->numberColumns();
double * element = matrix_->getMutableElements();
int * row = matrix_->getMutableIndices();
CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
int * length = matrix_->getMutableVectorLengths();
double * cost = model->costRegion();
double * lowerColumn = model->lowerRegion();
double * upperColumn = model->upperRegion();
int * pivotVariable = model->pivotVariable();
CoinBigIndex numberElements=startColumn[firstDynamic_];
// first just do lookup and basic stuff
int currentNumber = firstAvailable_;
firstAvailable_ = firstDynamic_;
int numberToDo=0;
double objectiveChange=0.0;
double * solution = model->solutionRegion();
for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
int iSet = backward_[iColumn];
if (toIndex_[iSet]<0) {
toIndex_[iSet]=0;
fromIndex_[numberToDo++]=iSet;
}
if (model->getStatus(iColumn)==ClpSimplex::basic||iColumn==keyVariable_[iSet]) {
lookup[iColumn]=firstAvailable_;
if (iColumn!=keyVariable_[iSet]) {
int iPivot = backToPivotRow_[iColumn];
backToPivotRow_[firstAvailable_]=iPivot;
pivotVariable[iPivot]=firstAvailable_;
}
firstAvailable_++;
} else {
int jColumn = id_[iColumn-firstDynamic_];
setDynamicStatus(jColumn,atLowerBound);
if (lowerColumn_||upperColumn_) {
if (model->getStatus(iColumn)==ClpSimplex::atUpperBound)
setDynamicStatus(jColumn,atUpperBound);
// treat solution as if exactly at a bound
double value = solution[iColumn];
if (fabs(value-lowerColumn[iColumn])<fabs(value-upperColumn[iColumn]))
value = lowerColumn[iColumn];
else
value = upperColumn[iColumn];
objectiveChange += cost[iColumn]*value;
// redo lower and upper on sets
double shift=value;
if (lowerSet_[iSet]>-1.0e20)
lower_[iSet] = lowerSet_[iSet]-shift;
if (upperSet_[iSet]<1.0e20)
upper_[iSet] = upperSet_[iSet]-shift;
}
lookup[iColumn]=-1;
}
}
model->setObjectiveOffset(model->objectiveOffset()+objectiveChange);
firstAvailable_ = firstDynamic_;
for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
if (lookup[iColumn]>=0) {
// move
int jColumn = id_[iColumn-firstDynamic_];
id_[firstAvailable_-firstDynamic_] = jColumn;
int numberThis = startColumn_[jColumn+1]-startColumn_[jColumn];
length[firstAvailable_]=numberThis;
cost[firstAvailable_]=cost[iColumn];
lowerColumn[firstAvailable_]=lowerColumn[iColumn];
upperColumn[firstAvailable_]=upperColumn[iColumn];
double originalLower = lowerColumn_ ? lowerColumn_[jColumn] : 0.0;
double originalUpper = upperColumn_ ? upperColumn_[jColumn] : COIN_DBL_MAX;
if (originalUpper>1.0e30)
originalUpper = COIN_DBL_MAX;
model->nonLinearCost()->setOne(firstAvailable_,solution[iColumn],
originalLower,originalUpper,
cost_[jColumn]);
CoinBigIndex base = startColumn_[jColumn];
for (int j=0;j<numberThis;j++) {
row[numberElements]=row_[base+j];
element[numberElements++]=element_[base+j];
}
model->setStatus(firstAvailable_,model->getStatus(iColumn));
backward_[firstAvailable_] = backward_[iColumn];
solution[firstAvailable_] = solution[iColumn];
firstAvailable_++;
startColumn[firstAvailable_]=numberElements;
}
}
// clean up next_
int * temp = new int [firstAvailable_];
for (int jSet=0;jSet<numberToDo;jSet++) {
int iSet = fromIndex_[jSet];
toIndex_[iSet]=-1;
int last = keyVariable_[iSet];
int j = next_[last];
bool setTemp=true;
if (last<lastDynamic_) {
last = lookup[last];
assert (last>=0);
keyVariable_[iSet]=last;
} else if (j>=0) {
int newJ = lookup[j];
assert (newJ>=0);
j=next_[j];
next_[last]=newJ;
last = newJ;
} else {
next_[last] = -(iSet+numberColumns+1);
setTemp=false;
}
while (j>=0) {
int newJ = lookup[j];
assert (newJ>=0);
temp[last]=newJ;
last = newJ;
j=next_[j];
}
if (setTemp)
temp[last]=-(keyVariable_[iSet]+1);
if (lowerSet_) {
// we only need to get lower_ and upper_ correct
double shift=0.0;
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
if (getDynamicStatus(j)==atUpperBound)
shift += upperColumn_[j];
else if (getDynamicStatus(j)==atLowerBound&&lowerColumn_)
shift += lowerColumn_[j];
if (lowerSet_[iSet]>-1.0e20)
lower_[iSet] = lowerSet_[iSet]-shift;
if (upperSet_[iSet]<1.0e20)
upper_[iSet] = upperSet_[iSet]-shift;
}
}
// move to next_
memcpy(next_+firstDynamic_,temp+firstDynamic_,(firstAvailable_-firstDynamic_)*sizeof(int));
// if odd iterations may be one out so adjust currentNumber
currentNumber=CoinMin(currentNumber+1,lastDynamic_);
// zero solution
CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
// zero cost
CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
// zero lengths
CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
model->setStatus(iColumn,ClpSimplex::atLowerBound);
backward_[iColumn]=-1;
}
delete [] lookup;
delete [] temp;
// make sure fromIndex clean
fromIndex_[0]=-1;
//#define CLP_DEBUG
#ifdef CLP_DEBUG
// debug
{
int i;
int numberRows=model->numberRows();
char * xxxx = new char[numberColumns];
memset(xxxx,0,numberColumns);
for (i=0;i<numberRows;i++) {
int iPivot = pivotVariable[i];
assert (model->getStatus(iPivot)==ClpSimplex::basic);
if (iPivot<numberColumns&&backward_[iPivot]>=0)
xxxx[iPivot]=1;
}
for (i=0;i<numberSets_;i++) {
int key=keyVariable_[i];
int iColumn =next_[key];
int k=0;
while(iColumn>=0) {
k++;
assert (k<100);
assert (backward_[iColumn]==i);
iColumn=next_[iColumn];
}
int stop = -(key+1);
while (iColumn!=stop) {
assert (iColumn<0);
iColumn = -iColumn-1;
k++;
assert (k<100);
assert (backward_[iColumn]==i);
iColumn=next_[iColumn];
}
iColumn =next_[key];
while (iColumn>=0) {
assert (xxxx[iColumn]);
xxxx[iColumn]=0;
iColumn=next_[iColumn];
}
}
for (i=0;i<numberColumns;i++) {
if (i<numberColumns&&backward_[i]>=0) {
assert (!xxxx[i]||i==keyVariable_[backward_[i]]);
}
}
delete [] xxxx;
}
{
for (int i=0;i<numberSets_;i++)
assert(toIndex_[i]==-1);
}
#endif
savedFirstAvailable_ = firstAvailable_;
}
break;
// flag a variable
case 1:
{
// id will be sitting at firstAvailable
int sequence = id_[firstAvailable_-firstDynamic_];
assert (!flagged(sequence));
setFlagged(sequence);
model->clearFlagged(firstAvailable_);
}
break;
// unflag all variables
case 2:
{
for (int i=0;i<numberGubColumns_;i++) {
if (flagged(i)) {
unsetFlagged(i);
returnNumber++;
}
}
}
break;
// just reset costs and bounds (primal)
case 3:
{
double * cost = model->costRegion();
double * solution = model->solutionRegion();
double * lowerColumn = model->columnLower();
double * upperColumn = model->columnUpper();
for (int i=firstDynamic_;i<firstAvailable_;i++) {
int jColumn = id_[i-firstDynamic_];
cost[i]=cost_[jColumn];
if (!lowerColumn_&&!upperColumn_) {
lowerColumn[i] = 0.0;
upperColumn[i] = COIN_DBL_MAX;
} else {
if (lowerColumn_)
lowerColumn[i] = lowerColumn_[jColumn];
else
lowerColumn[i] = 0.0;
if (upperColumn_)
upperColumn[i] = upperColumn_[jColumn];
else
upperColumn[i] = COIN_DBL_MAX;
}
if (model->nonLinearCost())
model->nonLinearCost()->setOne(i,solution[i],
lowerColumn[i],
upperColumn[i],cost_[jColumn]);
}
if (!model->numberIterations()&&rhsOffset_) {
lastRefresh_ = - refreshFrequency_; // force refresh
}
}
break;
// and get statistics for column generation
case 4:
{
// In theory we should subtract out ones we have done but ....
// If key slack then dual 0.0
// If not then slack could be dual infeasible
// dj for key is zero so that defines dual on set
int i;
int numberColumns = model->numberColumns();
double * dual = model->dualRowSolution();
double infeasibilityCost = model->infeasibilityCost();
double dualTolerance = model->dualTolerance();
double relaxedTolerance=dualTolerance;
// we can't really trust infeasibilities if there is dual error
double error = CoinMin(1.0e-2,model->largestDualError());
// allow tolerance at least slightly bigger than standard
relaxedTolerance = relaxedTolerance + error;
// but we will be using difference
relaxedTolerance -= dualTolerance;
double objectiveOffset = 0.0;
for (i=0;i<numberSets_;i++) {
int kColumn = keyVariable_[i];
double value=0.0;
if (kColumn<numberColumns) {
kColumn = id_[kColumn-firstDynamic_];
// dj without set
value = cost_[kColumn];
for (CoinBigIndex j=startColumn_[kColumn];
j<startColumn_[kColumn+1];j++) {
int iRow = row_[j];
value -= dual[iRow]*element_[j];
}
double infeasibility = 0.0;
if (getStatus(i)==ClpSimplex::atLowerBound) {
if (-value>dualTolerance)
infeasibility=-value-dualTolerance;
} else if (getStatus(i)==ClpSimplex::atUpperBound) {
if (value>dualTolerance)
infeasibility=value-dualTolerance;
}
if (infeasibility>0.0) {
sumDualInfeasibilities_ += infeasibility;
if (infeasibility>relaxedTolerance)
sumOfRelaxedDualInfeasibilities_ += infeasibility;
numberDualInfeasibilities_ ++;
}
} else {
// slack key - may not be feasible
assert (getStatus(i)==ClpSimplex::basic);
// negative as -1.0 for slack
value=-weight(i)*infeasibilityCost;
}
// Now subtract out from all
for (CoinBigIndex k= fullStart_[i];k<fullStart_[i+1];k++) {
if (getDynamicStatus(k)!=inSmall) {
double djValue = cost_[k]-value;
for (CoinBigIndex j=startColumn_[k];
j<startColumn_[k+1];j++) {
int iRow = row_[j];
djValue -= dual[iRow]*element_[j];
}
double infeasibility=0.0;
double shift=0.0;
if (getDynamicStatus(k)==atLowerBound) {
if (lowerColumn_)
shift= lowerColumn_[k];
if (djValue<-dualTolerance)
infeasibility=-djValue-dualTolerance;
} else {
// at upper bound
shift= upperColumn_[k];
if (djValue>dualTolerance)
infeasibility=djValue-dualTolerance;
}
objectiveOffset += shift*cost_[k];
if (infeasibility>0.0) {
sumDualInfeasibilities_ += infeasibility;
if (infeasibility>relaxedTolerance)
sumOfRelaxedDualInfeasibilities_ += infeasibility;
numberDualInfeasibilities_ ++;
}
}
}
}
model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
}
break;
// see if time to re-factorize
case 5:
{
if (firstAvailable_>numberSets_+model->numberRows()+ model->factorizationFrequency())
returnNumber=4;
}
break;
// return 1 if there may be changing bounds on variable (column generation)
case 6:
{
returnNumber = (lowerColumn_!=NULL||upperColumn_!=NULL) ? 1 : 0;
#if 0
if (!returnNumber) {
// may be gub slacks
for (int i=0;i<numberSets_;i++) {
if (upper_[i]>lower_[i]) {
returnNumber=1;
break;
}
}
}
#endif
}
break;
// restore firstAvailable_
case 7:
{
int iColumn;
int * length = matrix_->getMutableVectorLengths();
double * cost = model->costRegion();
double * solution = model->solutionRegion();
int currentNumber=firstAvailable_;
firstAvailable_ = savedFirstAvailable_;
// zero solution
CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
// zero cost
CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
// zero lengths
CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
model->setStatus(iColumn,ClpSimplex::atLowerBound);
backward_[iColumn]=-1;
}
}
break;
// make sure set is clean
case 8:
{
int sequenceIn = model->sequenceIn();
if (sequenceIn<model->numberColumns()) {
int iSet = backward_[sequenceIn];
if (iSet>=0&&lowerSet_) {
// we only need to get lower_ and upper_ correct
double shift=0.0;
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
if (getDynamicStatus(j)==atUpperBound)
shift += upperColumn_[j];
else if (getDynamicStatus(j)==atLowerBound&&lowerColumn_)
shift += lowerColumn_[j];
if (lowerSet_[iSet]>-1.0e20)
lower_[iSet] = lowerSet_[iSet]-shift;
if (upperSet_[iSet]<1.0e20)
upper_[iSet] = upperSet_[iSet]-shift;
}
if (sequenceIn==firstAvailable_) {
// not really in small problem
int iBig=id_[sequenceIn-firstDynamic_];
if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound)
setDynamicStatus(iBig,atLowerBound);
else
setDynamicStatus(iBig,atUpperBound);
}
}
}
break;
// adjust lower,upper
case 9:
{
int sequenceIn = model->sequenceIn();
if (sequenceIn>=firstDynamic_&&sequenceIn<lastDynamic_&&lowerSet_) {
int iSet = backward_[sequenceIn];
assert (iSet>=0);
int inBig = id_[sequenceIn-firstDynamic_];
const double * solution = model->solutionRegion();
setDynamicStatus(inBig,inSmall);
if (lowerSet_[iSet]>-1.0e20)
lower_[iSet] += solution[sequenceIn];
if (upperSet_[iSet]<1.0e20)
upper_[iSet] += solution[sequenceIn];
model->setObjectiveOffset(model->objectiveOffset()-
solution[sequenceIn]*cost_[inBig]);
}
}
}
return returnNumber;
}
// Add a new variable to a set
void
ClpGubDynamicMatrix::insertNonBasic(int sequence, int iSet)
{
int last = keyVariable_[iSet];
int j=next_[last];
while (j>=0) {
last=j;
j=next_[j];
}
next_[last]=-(sequence+1);
next_[sequence]=j;
}
// Sets up an effective RHS and does gub crash if needed
void
ClpGubDynamicMatrix::useEffectiveRhs(ClpSimplex * model, bool cheapest)
{
// Do basis - cheapest or slack if feasible (unless cheapest set)
int longestSet=0;
int iSet;
for (iSet=0;iSet<numberSets_;iSet++)
longestSet = CoinMax(longestSet,fullStart_[iSet+1]-fullStart_[iSet]);
double * upper = new double[longestSet+1];
double * cost = new double[longestSet+1];
double * lower = new double[longestSet+1];
double * solution = new double[longestSet+1];
assert (!next_);
delete [] next_;
int numberColumns = model->numberColumns();
next_ = new int[numberColumns+numberSets_+CoinMax(2*longestSet,lastDynamic_-firstDynamic_)];
char * mark = new char[numberColumns];
memset(mark,0,numberColumns);
for (int iColumn=0;iColumn<numberColumns;iColumn++)
next_[iColumn]=COIN_INT_MAX;
int i;
int * keys = new int[numberSets_];
int * back = new int[numberGubColumns_];
CoinFillN(back,numberGubColumns_,-1);
for (i=0;i<numberSets_;i++)
keys[i]=COIN_INT_MAX;
delete [] dynamicStatus_;
dynamicStatus_ = new unsigned char [numberGubColumns_];
memset(dynamicStatus_,0,numberGubColumns_); // for clarity
for (i=0;i<numberGubColumns_;i++)
setDynamicStatus(i,atLowerBound);
// set up chains
for (i=firstDynamic_;i<lastDynamic_;i++){
if (id_[i-firstDynamic_]>=0) {
if (model->getStatus(i)==ClpSimplex::basic)
mark[i]=1;
int iSet = backward_[i];
assert (iSet>=0);
int iNext = keys[iSet];
next_[i]=iNext;
keys[iSet]=i;
back[id_[i-firstDynamic_]]=i;
} else {
model->setStatus(i,ClpSimplex::atLowerBound);
backward_[i]=-1;
}
}
double * columnSolution = model->solutionRegion();
int numberRows = getNumRows();
toIndex_ = new int[numberSets_];
for (iSet=0;iSet<numberSets_;iSet++)
toIndex_[iSet]=-1;
fromIndex_ = new int [numberRows+numberSets_];
double tolerance = model->primalTolerance();
double * element = matrix_->getMutableElements();
int * row = matrix_->getMutableIndices();
CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
int * length = matrix_->getMutableVectorLengths();
double objectiveOffset = 0.0;
for (iSet=0;iSet<numberSets_;iSet++) {
int j;
int numberBasic=0;
int iBasic=-1;
int iStart = fullStart_[iSet];
int iEnd=fullStart_[iSet+1];
// find one with smallest length
int smallest=numberRows+1;
double value=0.0;
j=keys[iSet];
while (j!=COIN_INT_MAX) {
if (model->getStatus(j)== ClpSimplex::basic) {
if (length[j]<smallest) {
smallest=length[j];
iBasic=j;
}
numberBasic++;
}
value += columnSolution[j];
j = next_[j];
}
bool done=false;
if (numberBasic>1||(numberBasic==1&&getStatus(iSet)==ClpSimplex::basic)) {
if (getStatus(iSet)==ClpSimplex::basic)
iBasic = iSet+numberColumns;// slack key - use
done=true;
} else if (numberBasic==1) {
// see if can be key
double thisSolution = columnSolution[iBasic];
if (thisSolution<0.0) {
value -= thisSolution;
thisSolution = 0.0;
columnSolution[iBasic]=thisSolution;
}
// try setting slack to a bound
assert (upper_[iSet]<1.0e20||lower_[iSet]>-1.0e20);
double cost1 = COIN_DBL_MAX;
int whichBound=-1;
if (upper_[iSet]<1.0e20) {
// try slack at ub
double newBasic = thisSolution +upper_[iSet]-value;
if (newBasic>=-tolerance) {
// can go
whichBound=1;
cost1 = newBasic*cost_[iBasic];
// But if exact then may be good solution
if (fabs(upper_[iSet]-value)<tolerance)
cost1=-COIN_DBL_MAX;
}
}
if (lower_[iSet]>-1.0e20) {
// try slack at lb
double newBasic = thisSolution +lower_[iSet]-value;
if (newBasic>=-tolerance) {
// can go but is it cheaper
double cost2 = newBasic*cost_[iBasic];
// But if exact then may be good solution
if (fabs(lower_[iSet]-value)<tolerance)
cost2=-COIN_DBL_MAX;
if (cost2<cost1)
whichBound=0;
}
}
if (whichBound!=-1) {
// key
done=true;
if (whichBound) {
// slack to upper
columnSolution[iBasic]=thisSolution + upper_[iSet]-value;
setStatus(iSet,ClpSimplex::atUpperBound);
} else {
// slack to lower
columnSolution[iBasic]=thisSolution + lower_[iSet]-value;
setStatus(iSet,ClpSimplex::atLowerBound);
}
}
}
if (!done) {
if (!cheapest) {
// see if slack can be key
if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
done=true;
setStatus(iSet,ClpSimplex::basic);
iBasic=iSet+numberColumns;
}
}
if (!done) {
// set non basic if there was one
if (iBasic>=0)
model->setStatus(iBasic,ClpSimplex::atLowerBound);
// find cheapest
int numberInSet = iEnd-iStart;
if (!lowerColumn_) {
CoinZeroN(lower,numberInSet);
} else {
for (int j=0;j<numberInSet;j++)
lower[j]= lowerColumn_[j+iStart];
}
if (!upperColumn_) {
CoinFillN(upper,numberInSet,COIN_DBL_MAX);
} else {
for (int j=0;j<numberInSet;j++)
upper[j]= upperColumn_[j+iStart];
}
CoinFillN(solution,numberInSet,0.0);
// and slack
iBasic=numberInSet;
solution[iBasic]=-value;
lower[iBasic]=-upper_[iSet];
upper[iBasic]=-lower_[iSet];
int kphase;
if (value>=lower_[iSet]-tolerance&&value<=upper_[iSet]+tolerance) {
// feasible
kphase=1;
cost[iBasic]=0.0;
for (int j=0;j<numberInSet;j++)
cost[j]= cost_[j+iStart];
} else {
// infeasible
kphase=0;
// remember bounds are flipped so opposite to natural
if (value<lower_[iSet]-tolerance)
cost[iBasic]=1.0;
else
cost[iBasic]=-1.0;
CoinZeroN(cost,numberInSet);
}
double dualTolerance =model->dualTolerance();
for (int iphase =kphase;iphase<2;iphase++) {
if (iphase) {
cost[numberInSet]=0.0;
for (int j=0;j<numberInSet;j++)
cost[j]= cost_[j+iStart];
}
// now do one row lp
bool improve=true;
while (improve) {
improve=false;
double dual = cost[iBasic];
int chosen =-1;
double best=dualTolerance;
int way=0;
for (int i=0;i<=numberInSet;i++) {
double dj = cost[i]-dual;
double improvement =0.0;
double distance=0.0;
if (iphase||i<numberInSet)
assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
if (dj>dualTolerance)
improvement = dj*(solution[i]-lower[i]);
else if (dj<-dualTolerance)
improvement = dj*(solution[i]-upper[i]);
if (improvement>best) {
best=improvement;
chosen=i;
if (dj<0.0) {
way = 1;
distance = upper[i]-solution[i];
} else {
way = -1;
distance = solution[i]-lower[i];
}
}
}
if (chosen>=0) {
improve=true;
// now see how far
if (way>0) {
// incoming increasing so basic decreasing
// if phase 0 then go to nearest bound
double distance=upper[chosen]-solution[chosen];
double basicDistance;
if (!iphase) {
assert (iBasic==numberInSet);
assert (solution[iBasic]>upper[iBasic]);
basicDistance = solution[iBasic]-upper[iBasic];
} else {
basicDistance = solution[iBasic]-lower[iBasic];
}
// need extra coding for unbounded
assert (CoinMin(distance,basicDistance)<1.0e20);
if (distance>basicDistance) {
// incoming becomes basic
solution[chosen] += basicDistance;
if (!iphase)
solution[iBasic]=upper[iBasic];
else
solution[iBasic]=lower[iBasic];
iBasic = chosen;
} else {
// flip
solution[chosen]=upper[chosen];
solution[iBasic] -= distance;
}
} else {
// incoming decreasing so basic increasing
// if phase 0 then go to nearest bound
double distance=solution[chosen]-lower[chosen];
double basicDistance;
if (!iphase) {
assert (iBasic==numberInSet);
assert (solution[iBasic]<lower[iBasic]);
basicDistance = lower[iBasic]-solution[iBasic];
} else {
basicDistance = upper[iBasic]-solution[iBasic];
}
// need extra coding for unbounded - for now just exit
if (CoinMin(distance,basicDistance)>1.0e20) {
printf("unbounded on set %d\n",iSet);
iphase=1;
iBasic=numberInSet;
break;
}
if (distance>basicDistance) {
// incoming becomes basic
solution[chosen] -= basicDistance;
if (!iphase)
solution[iBasic]=lower[iBasic];
else
solution[iBasic]=upper[iBasic];
iBasic = chosen;
} else {
// flip
solution[chosen]=lower[chosen];
solution[iBasic] += distance;
}
}
if (!iphase) {
if(iBasic<numberInSet)
break; // feasible
else if (solution[iBasic]>=lower[iBasic]&&
solution[iBasic]<=upper[iBasic])
break; // feasible (on flip)
}
}
}
}
// do solution i.e. bounds
if (lowerColumn_||upperColumn_) {
for (int j=0;j<numberInSet;j++) {
if (j!=iBasic) {
objectiveOffset += solution[j]*cost[j];
if (lowerColumn_&&upperColumn_) {
if (fabs(solution[j]-lowerColumn_[j+iStart])>
fabs(solution[j]-upperColumn_[j+iStart]))
setDynamicStatus(j+iStart,atUpperBound);
} else if (upperColumn_&&solution[j]>0.0) {
setDynamicStatus(j+iStart,atUpperBound);
} else {
setDynamicStatus(j+iStart,atLowerBound);
}
}
}
}
// convert iBasic back and do bounds
if (iBasic==numberInSet) {
// slack basic
setStatus(iSet,ClpSimplex::basic);
iBasic=iSet+numberColumns;
} else {
iBasic += fullStart_[iSet];
if (back[iBasic]>=0) {
// exists
iBasic = back[iBasic];
} else {
// create
CoinBigIndex numberElements = startColumn[firstAvailable_];
int numberThis = startColumn_[iBasic+1]-startColumn_[iBasic];
if (numberElements+numberThis>numberElements_) {
// need to redo
numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
matrix_->reserve(numberColumns,numberElements_);
element = matrix_->getMutableElements();
row = matrix_->getMutableIndices();
// these probably okay but be safe
startColumn = matrix_->getMutableVectorStarts();
length = matrix_->getMutableVectorLengths();
}
length[firstAvailable_]=numberThis;
model->costRegion()[firstAvailable_]=cost_[iBasic];
if (lowerColumn_)
model->lowerRegion()[firstAvailable_] = lowerColumn_[iBasic];
else
model->lowerRegion()[firstAvailable_] = 0.0;
if (upperColumn_)
model->upperRegion()[firstAvailable_] = upperColumn_[iBasic];
else
model->upperRegion()[firstAvailable_] = COIN_DBL_MAX;
columnSolution[firstAvailable_]=solution[iBasic-fullStart_[iSet]];
CoinBigIndex base = startColumn_[iBasic];
for (int j=0;j<numberThis;j++) {
row[numberElements]=row_[base+j];
element[numberElements++]=element_[base+j];
}
// already set startColumn[firstAvailable_]=numberElements;
id_[firstAvailable_-firstDynamic_]=iBasic;
setDynamicStatus(iBasic,inSmall);
backward_[firstAvailable_]=iSet;
iBasic=firstAvailable_;
firstAvailable_++;
startColumn[firstAvailable_]=numberElements;
}
model->setStatus(iBasic,ClpSimplex::basic);
// remember bounds flipped
if (upper[numberInSet]==lower[numberInSet])
setStatus(iSet,ClpSimplex::isFixed);
else if (solution[numberInSet]==upper[numberInSet])
setStatus(iSet,ClpSimplex::atLowerBound);
else if (solution[numberInSet]==lower[numberInSet])
setStatus(iSet,ClpSimplex::atUpperBound);
else
abort();
}
for (j=iStart;j<iEnd;j++) {
int iBack = back[j];
if (iBack>=0) {
if (model->getStatus(iBack)!=ClpSimplex::basic) {
int inSet=j-iStart;
columnSolution[iBack]=solution[inSet];
if (upper[inSet]==lower[inSet])
model->setStatus(iBack,ClpSimplex::isFixed);
else if (solution[inSet]==upper[inSet])
model->setStatus(iBack,ClpSimplex::atUpperBound);
else if (solution[inSet]==lower[inSet])
model->setStatus(iBack,ClpSimplex::atLowerBound);
}
}
}
}
}
keyVariable_[iSet]=iBasic;
}
model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
delete [] lower;
delete [] solution;
delete [] upper;
delete [] cost;
// make sure matrix is in good shape
matrix_->orderMatrix();
// create effective rhs
delete [] rhsOffset_;
rhsOffset_ = new double[numberRows];
// and redo chains
memset(mark,0,numberColumns);
for (int iColumnX=0;iColumnX<firstAvailable_;iColumnX++)
next_[iColumnX]=COIN_INT_MAX;
for (i=0;i<numberSets_;i++) {
keys[i]=COIN_INT_MAX;
int iKey = keyVariable_[i];
if (iKey<numberColumns)
model->setStatus(iKey,ClpSimplex::basic);
}
// set up chains
for (i=0;i<firstAvailable_;i++){
if (model->getStatus(i)==ClpSimplex::basic)
mark[i]=1;
int iSet = backward_[i];
if (iSet>=0) {
int iNext = keys[iSet];
next_[i]=iNext;
keys[iSet]=i;
}
}
for (i=0;i<numberSets_;i++) {
if (keys[i]!=COIN_INT_MAX) {
// something in set
int j;
if (getStatus(i)!=ClpSimplex::basic) {
// make sure fixed if it is
if (upper_[i]==lower_[i])
setStatus(i,ClpSimplex::isFixed);
// slack not key - choose one with smallest length
int smallest=numberRows+1;
int key=-1;
j = keys[i];
while (1) {
if (mark[j]&&length[j]<smallest) {
key=j;
smallest=length[j];
}
if (next_[j]!=COIN_INT_MAX) {
j = next_[j];
} else {
// correct end
next_[j]=-(keys[i]+1);
break;
}
}
if (key>=0) {
keyVariable_[i]=key;
} else {
// nothing basic - make slack key
//((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic);
// fudge to avoid const problem
status_[i]=1;
}
} else {
// slack key
keyVariable_[i]=numberColumns+i;
int j;
double sol=0.0;
j = keys[i];
while (1) {
sol += columnSolution[j];
if (next_[j]!=COIN_INT_MAX) {
j = next_[j];
} else {
// correct end
next_[j]=-(keys[i]+1);
break;
}
}
if (sol>upper_[i]+tolerance) {
setAbove(i);
} else if (sol<lower_[i]-tolerance) {
setBelow(i);
} else {
setFeasible(i);
}
}
// Create next_
int key = keyVariable_[i];
redoSet(model,key,keys[i],i);
} else {
// nothing in set!
next_[i+numberColumns]=-(i+numberColumns+1);
keyVariable_[i]=numberColumns+i;
double sol=0.0;
if (sol>upper_[i]+tolerance) {
setAbove(i);
} else if (sol<lower_[i]-tolerance) {
setBelow(i);
} else {
setFeasible(i);
}
}
}
delete [] keys;
delete [] mark;
delete [] back;
rhsOffset(model,true);
}
/* Returns effective RHS if it is being used. This is used for long problems
or big gub or anywhere where going through full columns is
expensive. This may re-compute */
double *
ClpGubDynamicMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,
bool check)
{
//forceRefresh=true;
//check=false;
#ifdef CLP_DEBUG
double * saveE=NULL;
if (rhsOffset_&&check) {
int numberRows = model->numberRows();
saveE = new double[numberRows];
}
#endif
if (rhsOffset_) {
#ifdef CLP_DEBUG
if (check) {
// no need - but check anyway
int numberRows = model->numberRows();
double * rhs = new double[numberRows];
int numberColumns = model->numberColumns();
int iRow;
CoinZeroN(rhs,numberRows);
// do ones at bounds before gub
const double * smallSolution = model->solutionRegion();
const double * element =matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
int iColumn;
for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
if (model->getStatus(iColumn)!=ClpSimplex::basic) {
double value = smallSolution[iColumn];
for (CoinBigIndex j=startColumn[iColumn];
j<startColumn[iColumn]+length[iColumn];j++) {
int jRow=row[j];
rhs[jRow] -= value*element[j];
}
}
}
if (lowerColumn_||upperColumn_) {
double * solution = new double [numberGubColumns_];
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value=0.0;
if(getDynamicStatus(iColumn)==atUpperBound)
value = upperColumn_[iColumn];
else if (lowerColumn_)
value = lowerColumn_[iColumn];
solution[iColumn]=value;
}
// ones at bounds in small and gub
for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
int jFull = id_[iColumn-firstDynamic_];
solution[jFull]=smallSolution[iColumn];
}
// zero all basic in small model
int * pivotVariable = model->pivotVariable();
for (iRow=0;iRow<numberRows;iRow++) {
int iColumn = pivotVariable[iRow];
if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=0.0;
}
}
// and now compute value to use for key
ClpSimplex::Status iStatus;
for (int iSet=0;iSet<numberSets_;iSet++) {
iColumn = keyVariable_[iSet];
if (iColumn<numberColumns) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=0.0;
double b=0.0;
// key is structural - where is slack
iStatus = getStatus(iSet);
assert (iStatus!=ClpSimplex::basic);
if (iStatus==ClpSimplex::atLowerBound)
b=lowerSet_[iSet];
else
b=upperSet_[iSet];
// subtract out others at bounds
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
b -= solution[j];
solution[iSequence]=b;
}
}
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value = solution[iColumn];
if (value) {
for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
int iRow = row_[j];
rhs[iRow] -= element_[j]*value;
}
}
}
// now do lower and upper bounds on sets
for (int iSet=0;iSet<numberSets_;iSet++) {
iColumn = keyVariable_[iSet];
double shift=0.0;
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++) {
if (getDynamicStatus(j)!=inSmall&&j!=iColumn) {
if (getDynamicStatus(j)==atLowerBound) {
if (lowerColumn_)
shift += lowerColumn_[j];
} else {
shift += upperColumn_[j];
}
}
}
if (lowerSet_[iSet]>-1.0e20)
assert(fabs(lower_[iSet] - (lowerSet_[iSet]-shift))<1.0e-3);
if (upperSet_[iSet]<1.0e20)
assert(fabs(upper_[iSet] -( upperSet_[iSet]-shift))<1.0e-3);
}
delete [] solution;
} else {
// no bounds
ClpSimplex::Status iStatus;
for (int iSet=0;iSet<numberSets_;iSet++) {
int iColumn = keyVariable_[iSet];
if (iColumn<numberColumns) {
int iSequence = id_[iColumn-firstDynamic_];
double b=0.0;
// key is structural - where is slack
iStatus = getStatus(iSet);
assert (iStatus!=ClpSimplex::basic);
if (iStatus==ClpSimplex::atLowerBound)
b=lower_[iSet];
else
b=upper_[iSet];
if (b) {
for (CoinBigIndex j= startColumn_[iSequence];j<startColumn_[iSequence+1];j++) {
int iRow = row_[j];
rhs[iRow] -= element_[j]*b;
}
}
}
}
}
for (iRow=0;iRow<numberRows;iRow++) {
if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
}
memcpy(saveE,rhs,numberRows*sizeof(double));
delete [] rhs;
}
#endif
if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
lastRefresh_+refreshFrequency_)) {
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
int iRow;
CoinZeroN(rhsOffset_,numberRows);
// do ones at bounds before gub
const double * smallSolution = model->solutionRegion();
const double * element =matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
int iColumn;
for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
if (model->getStatus(iColumn)!=ClpSimplex::basic) {
double value = smallSolution[iColumn];
for (CoinBigIndex j=startColumn[iColumn];
j<startColumn[iColumn]+length[iColumn];j++) {
int jRow=row[j];
rhsOffset_[jRow] -= value*element[j];
}
}
}
if (lowerColumn_||upperColumn_) {
double * solution = new double [numberGubColumns_];
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value=0.0;
if(getDynamicStatus(iColumn)==atUpperBound)
value = upperColumn_[iColumn];
else if (lowerColumn_)
value = lowerColumn_[iColumn];
solution[iColumn]=value;
}
// ones in gub and in small problem
for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
int jFull = id_[iColumn-firstDynamic_];
solution[jFull]=smallSolution[iColumn];
}
// zero all basic in small model
int * pivotVariable = model->pivotVariable();
for (iRow=0;iRow<numberRows;iRow++) {
int iColumn = pivotVariable[iRow];
if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=0.0;
}
}
// and now compute value to use for key
ClpSimplex::Status iStatus;
int iSet;
for ( iSet=0;iSet<numberSets_;iSet++) {
iColumn = keyVariable_[iSet];
if (iColumn<numberColumns) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=0.0;
double b=0.0;
// key is structural - where is slack
iStatus = getStatus(iSet);
assert (iStatus!=ClpSimplex::basic);
if (iStatus==ClpSimplex::atLowerBound)
b=lowerSet_[iSet];
else
b=upperSet_[iSet];
// subtract out others at bounds
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
b -= solution[j];
solution[iSequence]=b;
}
}
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value = solution[iColumn];
if (value) {
for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
int iRow = row_[j];
rhsOffset_[iRow] -= element_[j]*value;
}
}
}
// now do lower and upper bounds on sets
// and offset
double objectiveOffset = 0.0;
for ( iSet=0;iSet<numberSets_;iSet++) {
iColumn = keyVariable_[iSet];
double shift=0.0;
for (CoinBigIndex j=fullStart_[iSet];j<fullStart_[iSet+1];j++) {
if (getDynamicStatus(j)!=inSmall) {
double value=0.0;
if (getDynamicStatus(j)==atLowerBound) {
if (lowerColumn_)
value= lowerColumn_[j];
} else {
value= upperColumn_[j];
}
if (j!=iColumn)
shift += value;
objectiveOffset += value*cost_[j];
}
}
if (lowerSet_[iSet]>-1.0e20)
lower_[iSet] = lowerSet_[iSet]-shift;
if (upperSet_[iSet]<1.0e20)
upper_[iSet] = upperSet_[iSet]-shift;
}
delete [] solution;
model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
} else {
// no bounds
ClpSimplex::Status iStatus;
for (int iSet=0;iSet<numberSets_;iSet++) {
int iColumn = keyVariable_[iSet];
if (iColumn<numberColumns) {
int iSequence = id_[iColumn-firstDynamic_];
double b=0.0;
// key is structural - where is slack
iStatus = getStatus(iSet);
assert (iStatus!=ClpSimplex::basic);
if (iStatus==ClpSimplex::atLowerBound)
b=lower_[iSet];
else
b=upper_[iSet];
if (b) {
for (CoinBigIndex j= startColumn_[iSequence];j<startColumn_[iSequence+1];j++) {
int iRow = row_[j];
rhsOffset_[iRow] -= element_[j]*b;
}
}
}
}
}
#ifdef CLP_DEBUG
if (saveE) {
for (iRow=0;iRow<numberRows;iRow++) {
if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3)
printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]);
}
delete [] saveE;
}
#endif
lastRefresh_ = model->numberIterations();
}
}
return rhsOffset_;
}
/*
update information for a pivot (and effective rhs)
*/
int
ClpGubDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
{
// now update working model
int sequenceIn = model->sequenceIn();
int sequenceOut = model->sequenceOut();
bool doPrinting = (model->messageHandler()->logLevel()==63);
bool print=false;
int iSet;
int trueIn=-1;
int trueOut=-1;
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
if (sequenceIn==firstAvailable_) {
if (doPrinting)
printf("New variable ");
if (sequenceIn!=sequenceOut) {
insertNonBasic(firstAvailable_,backward_[firstAvailable_]);
setDynamicStatus(id_[sequenceIn-firstDynamic_],inSmall);
firstAvailable_++;
} else {
int bigSequence = id_[sequenceIn-firstDynamic_];
if (model->getStatus(sequenceIn)==ClpSimplex::atUpperBound)
setDynamicStatus(bigSequence,atUpperBound);
else
setDynamicStatus(bigSequence,atLowerBound);
}
synchronize(model,8);
}
if (sequenceIn<lastDynamic_) {
iSet = backward_[sequenceIn];
if (iSet>=0) {
int bigSequence = id_[sequenceIn-firstDynamic_];
trueIn=bigSequence+numberRows+numberColumns+numberSets_;
if (doPrinting)
printf(" incoming set %d big seq %d",iSet,bigSequence);
print=true;
}
} else if (sequenceIn>=numberRows+numberColumns) {
trueIn = numberRows+numberColumns+gubSlackIn_;
}
if (sequenceOut<lastDynamic_) {
iSet = backward_[sequenceOut];
if (iSet>=0) {
int bigSequence = id_[sequenceOut-firstDynamic_];
trueOut=bigSequence+firstDynamic_;
if (getDynamicStatus(bigSequence)!=inSmall) {
if (model->getStatus(sequenceOut)==ClpSimplex::atUpperBound)
setDynamicStatus(bigSequence,atUpperBound);
else
setDynamicStatus(bigSequence,atLowerBound);
}
if (doPrinting)
printf(" ,outgoing set %d big seq %d,",iSet,bigSequence);
print=true;
model->setSequenceIn(sequenceOut);
synchronize(model,8);
model->setSequenceIn(sequenceIn);
}
}
if (print&&doPrinting)
printf("\n");
ClpGubMatrix::updatePivot(model,oldInValue,oldOutValue);
// Redo true in and out
if (trueIn>=0)
trueSequenceIn_=trueIn;
if (trueOut>=0)
trueSequenceOut_=trueOut;
if (doPrinting&&0) {
for (int i=0;i<numberSets_;i++) {
printf("set %d key %d lower %g upper %g\n",i,keyVariable_[i],lower_[i],upper_[i]);
for (int j=fullStart_[i];j<fullStart_[i+1];j++)
if (getDynamicStatus(j)==atUpperBound) {
bool print=true;
for (int k=firstDynamic_;k<firstAvailable_;k++) {
if (id_[k-firstDynamic_]==j)
print=false;
if (id_[k-firstDynamic_]==j)
assert(getDynamicStatus(j)==inSmall);
}
if (print)
printf("variable %d at ub\n",j);
}
}
}
#ifdef CLP_DEBUG
char * inSmall = new char [numberGubColumns_];
memset(inSmall,0,numberGubColumns_);
for (int i=0;i<numberGubColumns_;i++)
if (getDynamicStatus(i)==ClpGubDynamicMatrix::inSmall)
inSmall[i]=1;
for (int i=firstDynamic_;i<firstAvailable_;i++) {
int k=id_[i-firstDynamic_];
inSmall[k]=0;
}
for (int i=0;i<numberGubColumns_;i++)
assert (!inSmall[i]);
delete [] inSmall;
#endif
return 0;
}
void
ClpGubDynamicMatrix::times(double scalar,
const double * x, double * y) const
{
if (model_->specialOptions()!=16) {
ClpPackedMatrix::times(scalar,x,y);
} else {
int iRow;
int numberColumns = model_->numberColumns();
int numberRows = model_->numberRows();
const double * element = matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
int * pivotVariable = model_->pivotVariable();
int numberToDo=0;
for (iRow=0;iRow<numberRows;iRow++) {
y[iRow] -= scalar*rhsOffset_[iRow];
int iColumn = pivotVariable[iRow];
if (iColumn<numberColumns) {
int iSet = backward_[iColumn];
if (iSet>=0&&toIndex_[iSet]<0) {
toIndex_[iSet]=0;
fromIndex_[numberToDo++]=iSet;
}
CoinBigIndex j;
double value = scalar*x[iColumn];
if (value) {
for (j=startColumn[iColumn];
j<startColumn[iColumn]+length[iColumn];j++) {
int jRow=row[j];
y[jRow] += value*element[j];
}
}
}
}
// and gubs which are interacting
for (int jSet=0;jSet<numberToDo;jSet++) {
int iSet = fromIndex_[jSet];
toIndex_[iSet]=-1;
int iKey=keyVariable_[iSet];
if (iKey<numberColumns) {
double valueKey;
if (getStatus(iSet)==ClpSimplex::atLowerBound)
valueKey = lower_[iSet];
else
valueKey = upper_[iSet];
double value = scalar*(x[iKey]-valueKey);
if (value) {
for (CoinBigIndex j=startColumn[iKey];
j<startColumn[iKey]+length[iKey];j++) {
int jRow=row[j];
y[jRow] += value*element[j];
}
}
}
}
}
}
/* Just for debug - may be extended to other matrix types later.
Returns number and sum of primal infeasibilities.
*/
int
ClpGubDynamicMatrix::checkFeasible(ClpSimplex * model,double & sum) const
{
int numberRows = model_->numberRows();
double * rhs = new double[numberRows];
int numberColumns = model_->numberColumns();
int iRow;
CoinZeroN(rhs,numberRows);
// do ones at bounds before gub
const double * smallSolution = model_->solutionRegion();
const double * element =matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
int iColumn;
int numberInfeasible=0;
const double * rowLower = model_->rowLower();
const double * rowUpper = model_->rowUpper();
sum=0.0;
for (iRow=0;iRow<numberRows;iRow++) {
double value=smallSolution[numberColumns+iRow];
if (value<rowLower[iRow]-1.0e-5||
value>rowUpper[iRow]+1.0e-5) {
//printf("row %d %g %g %g\n",
// iRow,rowLower[iRow],value,rowUpper[iRow]);
numberInfeasible++;
sum += CoinMax(rowLower[iRow]-value,value-rowUpper[iRow]);
}
rhs[iRow]=value;
}
const double * columnLower = model_->columnLower();
const double * columnUpper = model_->columnUpper();
for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
double value=smallSolution[iColumn];
if (value<columnLower[iColumn]-1.0e-5||
value>columnUpper[iColumn]+1.0e-5) {
//printf("column %d %g %g %g\n",
// iColumn,columnLower[iColumn],value,columnUpper[iColumn]);
numberInfeasible++;
sum += CoinMax(columnLower[iColumn]-value,value-columnUpper[iColumn]);
}
for (CoinBigIndex j=startColumn[iColumn];
j<startColumn[iColumn]+length[iColumn];j++) {
int jRow=row[j];
rhs[jRow] -= value*element[j];
}
}
double * solution = new double [numberGubColumns_];
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value=0.0;
if(getDynamicStatus(iColumn)==atUpperBound)
value = upperColumn_[iColumn];
else if (lowerColumn_)
value = lowerColumn_[iColumn];
solution[iColumn]=value;
}
// ones in small and gub
for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
int jFull = id_[iColumn-firstDynamic_];
solution[jFull]=smallSolution[iColumn];
}
// fill in all basic in small model
int * pivotVariable = model_->pivotVariable();
for (iRow=0;iRow<numberRows;iRow++) {
int iColumn = pivotVariable[iRow];
if (iColumn>=firstDynamic_&&iColumn<lastDynamic_) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=smallSolution[iColumn];
}
}
// and now compute value to use for key
ClpSimplex::Status iStatus;
for (int iSet=0;iSet<numberSets_;iSet++) {
iColumn = keyVariable_[iSet];
if (iColumn<numberColumns) {
int iSequence = id_[iColumn-firstDynamic_];
solution[iSequence]=0.0;
double b=0.0;
// key is structural - where is slack
iStatus = getStatus(iSet);
assert (iStatus!=ClpSimplex::basic);
if (iStatus==ClpSimplex::atLowerBound)
b=lower_[iSet];
else
b=upper_[iSet];
// subtract out others at bounds
for (int j=fullStart_[iSet];j<fullStart_[iSet+1];j++)
b -= solution[j];
solution[iSequence]=b;
}
}
for (iColumn=0;iColumn<numberGubColumns_;iColumn++) {
double value = solution[iColumn];
if ((lowerColumn_&&value<lowerColumn_[iColumn]-1.0e-5)||
(!lowerColumn_&&value<-1.0e-5)||
(upperColumn_&&value>upperColumn_[iColumn]+1.0e-5)) {
//printf("column %d %g %g %g\n",
// iColumn,lowerColumn_[iColumn],value,upperColumn_[iColumn]);
numberInfeasible++;
}
if (value) {
for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
int iRow = row_[j];
rhs[iRow] -= element_[j]*value;
}
}
}
for (iRow=0;iRow<numberRows;iRow++) {
if (fabs(rhs[iRow])>1.0e-5)
printf("rhs mismatch %d %g\n",iRow,rhs[iRow]);
}
delete [] solution;
delete [] rhs;
return numberInfeasible;
}
// Cleans data after setWarmStart
void
ClpGubDynamicMatrix::cleanData(ClpSimplex * model)
{
// and redo chains
int numberColumns = model->numberColumns();
int iColumn;
// do backward
int * mark = new int [numberGubColumns_];
for (iColumn=0;iColumn<numberGubColumns_;iColumn++)
mark[iColumn]=-1;
int i;
for (i=0;i<firstDynamic_;i++) {
assert (backward_[i]==-1);
next_[i]=-1;
}
for (i=firstDynamic_;i<firstAvailable_;i++) {
iColumn = id_[i-firstDynamic_];
mark[iColumn]=i;
}
for (i=0;i<numberSets_;i++) {
int iKey=keyVariable_[i];
int lastNext = -1;
int firstNext = -1;
for (CoinBigIndex k= fullStart_[i];k<fullStart_[i+1];k++) {
iColumn = mark[k];
if (iColumn>=0) {
if (iColumn!=iKey) {
if (lastNext>=0)
next_[lastNext]=iColumn;
else
firstNext = iColumn;
lastNext=iColumn;
}
backward_[iColumn]=i;
}
}
setFeasible(i);
if (firstNext>=0) {
// others
next_[iKey]=firstNext;
next_[lastNext]=-(iKey+1);
} else if (iKey<numberColumns) {
next_[iKey]=-(iKey+1);
}
}
delete [] mark;
// fill matrix
double * element = matrix_->getMutableElements();
int * row = matrix_->getMutableIndices();
CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
int * length = matrix_->getMutableVectorLengths();
CoinBigIndex numberElements = startColumn[firstDynamic_];
for (i=firstDynamic_;i<firstAvailable_;i++) {
int iColumn = id_[i-firstDynamic_];
int numberThis = startColumn_[iColumn+1]-startColumn_[iColumn];
length[i]=numberThis;
for (CoinBigIndex jBigIndex=startColumn_[iColumn];
jBigIndex<startColumn_[iColumn+1];jBigIndex++) {
row[numberElements] = row_[jBigIndex];
element[numberElements++] = element_[jBigIndex];
}
startColumn[i+1]=numberElements;
}
}
syntax highlighted by Code2HTML, v. 0.9.1