// Copyright (C) 2003, International Business Machines
// Corporation and others. All Rights Reserved.
#include "CoinPragma.hpp"
#include "ClpNetworkBasis.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpSimplex.hpp"
#include "ClpMatrixBase.hpp"
#include "CoinIndexedVector.hpp"
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpNetworkBasis::ClpNetworkBasis ()
{
slackValue_=-1.0;
numberRows_=0;
numberColumns_=0;
parent_ = NULL;
descendant_ = NULL;
pivot_ = NULL;
rightSibling_ = NULL;
leftSibling_ = NULL;
sign_ = NULL;
stack_ = NULL;
permute_ = NULL;
permuteBack_ = NULL;
stack2_=NULL;
depth_ = NULL;
mark_ = NULL;
model_=NULL;
}
// Constructor from CoinFactorization
ClpNetworkBasis::ClpNetworkBasis(const ClpSimplex * model,
int numberRows, const double * pivotRegion,
const int * permuteBack,
const CoinBigIndex * startColumn,
const int * numberInColumn,
const int * indexRow, const double * element)
{
slackValue_=-1.0;
numberRows_=numberRows;
numberColumns_=numberRows;
parent_ = new int [ numberRows_+1];
descendant_ = new int [ numberRows_+1];
pivot_ = new int [ numberRows_+1];
rightSibling_ = new int [ numberRows_+1];
leftSibling_ = new int [ numberRows_+1];
sign_ = new double [ numberRows_+1];
stack_ = new int [ numberRows_+1];
stack2_ = new int[numberRows_+1];
depth_ = new int[numberRows_+1];
mark_ = new char[numberRows_+1];
permute_ = new int [numberRows_ + 1];
permuteBack_ = new int [numberRows_ + 1];
int i;
for (i=0;i<numberRows_+1;i++) {
parent_[i]=-1;
descendant_[i]=-1;
pivot_[i]=-1;
rightSibling_[i]=-1;
leftSibling_[i]=-1;
sign_[i]=-1.0;
stack_[i]=-1;
permute_[i]=i;
permuteBack_[i]=i;
stack2_[i]=-1;
depth_[i]=-1;
mark_[i]=0;
}
mark_[numberRows_]=1;
// pivotColumnBack gives order of pivoting into basis
// so pivotColumnback[0] is first slack in basis and
// it pivots on row permuteBack[0]
// a known root is given by permuteBack[numberRows_-1]
int lastPivot=numberRows_;
for (i=0;i<numberRows_;i++) {
int iPivot = permuteBack[i];
lastPivot=iPivot;
double sign;
if (pivotRegion[i]>0.0)
sign = 1.0;
else
sign =-1.0;
int other;
if (numberInColumn[i]>0) {
int iRow = indexRow[startColumn[i]];
other = permuteBack[iRow];
//assert (parent_[other]!=-1);
} else {
other = numberRows_;
}
sign_[iPivot] = sign;
int iParent = other;
parent_[iPivot] = other;
if (descendant_[iParent]>=0) {
// we have a sibling
int iRight = descendant_[iParent];
rightSibling_[iPivot]=iRight;
leftSibling_[iRight]=iPivot;
} else {
rightSibling_[iPivot]=-1;
}
descendant_[iParent] = iPivot;
leftSibling_[iPivot]=-1;
}
// do depth
int iPivot = numberRows_;
int nStack = 1;
stack_[0]=descendant_[numberRows_];
depth_[numberRows_]=-1; // root
while (nStack) {
// take off
int iNext = stack_[--nStack];
if (iNext>=0) {
depth_[iNext] = nStack;
iPivot = iNext;
int iRight = rightSibling_[iNext];
stack_[nStack++] = iRight;
if (descendant_[iNext]>=0)
stack_[nStack++] = descendant_[iNext];
}
}
model_=model;
check();
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpNetworkBasis::ClpNetworkBasis (const ClpNetworkBasis & rhs)
{
slackValue_=rhs.slackValue_;
numberRows_=rhs.numberRows_;
numberColumns_=rhs.numberColumns_;
if (rhs.parent_) {
parent_ = new int [numberRows_+1];
memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
} else {
parent_ = NULL;
}
if (rhs.descendant_) {
descendant_ = new int [numberRows_+1];
memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
} else {
descendant_ = NULL;
}
if (rhs.pivot_) {
pivot_ = new int [numberRows_+1];
memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
} else {
pivot_ = NULL;
}
if (rhs.rightSibling_) {
rightSibling_ = new int [numberRows_+1];
memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
} else {
rightSibling_ = NULL;
}
if (rhs.leftSibling_) {
leftSibling_ = new int [numberRows_+1];
memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
} else {
leftSibling_ = NULL;
}
if (rhs.sign_) {
sign_ = new double [numberRows_+1];
memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
} else {
sign_ = NULL;
}
if (rhs.stack_) {
stack_ = new int [numberRows_+1];
memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
} else {
stack_ = NULL;
}
if (rhs.permute_) {
permute_ = new int [numberRows_+1];
memcpy(permute_,rhs.permute_,(numberRows_+1)*sizeof(int));
} else {
permute_ = NULL;
}
if (rhs.permuteBack_) {
permuteBack_ = new int [numberRows_+1];
memcpy(permuteBack_,rhs.permuteBack_,(numberRows_+1)*sizeof(int));
} else {
permuteBack_ = NULL;
}
if (rhs.stack2_) {
stack2_ = new int [numberRows_+1];
memcpy(stack2_,rhs.stack2_,(numberRows_+1)*sizeof(int));
} else {
stack2_ = NULL;
}
if (rhs.depth_) {
depth_ = new int [numberRows_+1];
memcpy(depth_,rhs.depth_,(numberRows_+1)*sizeof(int));
} else {
depth_ = NULL;
}
if (rhs.mark_) {
mark_ = new char [numberRows_+1];
memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
} else {
mark_ = NULL;
}
model_=rhs.model_;
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpNetworkBasis::~ClpNetworkBasis ()
{
delete [] parent_;
delete [] descendant_;
delete [] pivot_;
delete [] rightSibling_;
delete [] leftSibling_;
delete [] sign_;
delete [] stack_;
delete [] permute_;
delete [] permuteBack_;
delete [] stack2_;
delete [] depth_;
delete [] mark_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpNetworkBasis &
ClpNetworkBasis::operator=(const ClpNetworkBasis& rhs)
{
if (this != &rhs) {
delete [] parent_;
delete [] descendant_;
delete [] pivot_;
delete [] rightSibling_;
delete [] leftSibling_;
delete [] sign_;
delete [] stack_;
delete [] permute_;
delete [] permuteBack_;
delete [] stack2_;
delete [] depth_;
delete [] mark_;
slackValue_=rhs.slackValue_;
numberRows_=rhs.numberRows_;
numberColumns_=rhs.numberColumns_;
if (rhs.parent_) {
parent_ = new int [numberRows_+1];
memcpy(parent_,rhs.parent_,(numberRows_+1)*sizeof(int));
} else {
parent_ = NULL;
}
if (rhs.descendant_) {
descendant_ = new int [numberRows_+1];
memcpy(descendant_,rhs.descendant_,(numberRows_+1)*sizeof(int));
} else {
descendant_ = NULL;
}
if (rhs.pivot_) {
pivot_ = new int [numberRows_+1];
memcpy(pivot_,rhs.pivot_,(numberRows_+1)*sizeof(int));
} else {
pivot_ = NULL;
}
if (rhs.rightSibling_) {
rightSibling_ = new int [numberRows_+1];
memcpy(rightSibling_,rhs.rightSibling_,(numberRows_+1)*sizeof(int));
} else {
rightSibling_ = NULL;
}
if (rhs.leftSibling_) {
leftSibling_ = new int [numberRows_+1];
memcpy(leftSibling_,rhs.leftSibling_,(numberRows_+1)*sizeof(int));
} else {
leftSibling_ = NULL;
}
if (rhs.sign_) {
sign_ = new double [numberRows_+1];
memcpy(sign_,rhs.sign_,(numberRows_+1)*sizeof(double));
} else {
sign_ = NULL;
}
if (rhs.stack_) {
stack_ = new int [numberRows_+1];
memcpy(stack_,rhs.stack_,(numberRows_+1)*sizeof(int));
} else {
stack_ = NULL;
}
if (rhs.permute_) {
permute_ = new int [numberRows_+1];
memcpy(permute_,rhs.permute_,(numberRows_+1)*sizeof(int));
} else {
permute_ = NULL;
}
if (rhs.permuteBack_) {
permuteBack_ = new int [numberRows_+1];
memcpy(permuteBack_,rhs.permuteBack_,(numberRows_+1)*sizeof(int));
} else {
permuteBack_ = NULL;
}
if (rhs.stack2_) {
stack2_ = new int [numberRows_+1];
memcpy(stack2_,rhs.stack2_,(numberRows_+1)*sizeof(int));
} else {
stack2_ = NULL;
}
if (rhs.depth_) {
depth_ = new int [numberRows_+1];
memcpy(depth_,rhs.depth_,(numberRows_+1)*sizeof(int));
} else {
depth_ = NULL;
}
if (rhs.mark_) {
mark_ = new char [numberRows_+1];
memcpy(mark_,rhs.mark_,(numberRows_+1)*sizeof(char));
} else {
mark_ = NULL;
}
}
return *this;
}
// checks looks okay
void ClpNetworkBasis::check()
{
// check depth
int iPivot = numberRows_;
int nStack = 1;
stack_[0]=descendant_[numberRows_];
depth_[numberRows_]=-1; // root
while (nStack) {
// take off
int iNext = stack_[--nStack];
if (iNext>=0) {
//assert (depth_[iNext]==nStack);
depth_[iNext] = nStack;
iPivot = iNext;
int iRight = rightSibling_[iNext];
stack_[nStack++] = iRight;
if (descendant_[iNext]>=0)
stack_[nStack++] = descendant_[iNext];
}
}
}
// prints
void ClpNetworkBasis::print()
{
int i;
printf(" parent descendant left right sign depth\n");
for (i=0;i<numberRows_+1;i++)
printf("%4d %7d %8d %7d %7d %5g %7d\n",
i,parent_[i],descendant_[i],leftSibling_[i],rightSibling_[i],
sign_[i],depth_[i]);
}
/* Replaces one Column to basis,
returns 0=OK
*/
int
ClpNetworkBasis::replaceColumn ( CoinIndexedVector * regionSparse,
int pivotRow)
{
// When things have settled down then redo this to make more elegant
// I am sure lots of loops can be combined
// regionSparse is empty
assert (!regionSparse->getNumElements());
model_->unpack(regionSparse, model_->sequenceIn());
// arc given by pivotRow is leaving basis
//int kParent = parent_[pivotRow];
// arc coming in has these two nodes
int * indices = regionSparse->getIndices();
int iRow0 = indices[0];
int iRow1;
if (regionSparse->getNumElements()==2)
iRow1 = indices[1];
else
iRow1 = numberRows_;
double sign = -regionSparse->denseVector()[iRow0];
regionSparse->clear();
// and outgoing
model_->unpack(regionSparse,model_->pivotVariable()[pivotRow]);
int jRow0 = indices[0];
int jRow1;
if (regionSparse->getNumElements()==2)
jRow1 = indices[1];
else
jRow1 = numberRows_;
regionSparse->clear();
// get correct pivotRow
//#define FULL_DEBUG
#ifdef FULL_DEBUG
printf ("irow %d %d, jrow %d %d\n",
iRow0,iRow1,jRow0,jRow1);
#endif
if (parent_[jRow0]==jRow1) {
int newPivot = jRow0;
if (newPivot!=pivotRow) {
#ifdef FULL_DEBUG
printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
#endif
pivotRow = newPivot;
}
} else {
//assert (parent_[jRow1]==jRow0);
int newPivot = jRow1;
if (newPivot!=pivotRow) {
#ifdef FULL_DEBUG
printf("pivot row of %d permuted to %d\n",pivotRow,newPivot);
#endif
pivotRow = newPivot;
}
}
bool extraPrint = (model_->numberIterations()>-3)&&
(model_->logLevel()>10);
if (extraPrint)
print();
#ifdef FULL_DEBUG
printf("In %d (region= %g, stored %g) %d (%g) pivoting on %d (%g)\n",
iRow1, sign, sign_[iRow1], iRow0, sign_[iRow0] ,pivotRow,sign_[pivotRow]);
#endif
// see which path outgoing pivot is on
int kRow = -1;
int jRow = iRow1;
while (jRow!=numberRows_) {
if (jRow==pivotRow) {
kRow = iRow1;
break;
} else {
jRow = parent_[jRow];
}
}
if (kRow<0) {
jRow = iRow0;
while (jRow!=numberRows_) {
if (jRow==pivotRow) {
kRow = iRow0;
break;
} else {
jRow = parent_[jRow];
}
}
}
//assert (kRow>=0);
if (iRow0==kRow) {
iRow0 = iRow1;
iRow1 = kRow;
sign = -sign;
}
// pivot row is on path from iRow1 back to root
// get stack of nodes to change
// Also get precursors for cleaning order
int nStack=1;
stack_[0]=iRow0;
while (kRow!=pivotRow) {
stack_[nStack++] = kRow;
if (sign*sign_[kRow]<0.0) {
sign_[kRow]= -sign_[kRow];
} else {
sign = -sign;
}
kRow = parent_[kRow];
//sign *= sign_[kRow];
}
stack_[nStack++]=pivotRow;
if (sign*sign_[pivotRow]<0.0) {
sign_[pivotRow]= -sign_[pivotRow];
} else {
sign = -sign;
}
int iParent = parent_[pivotRow];
while (nStack>1) {
int iLeft;
int iRight;
kRow = stack_[--nStack];
int newParent = stack_[nStack-1];
#ifdef FULL_DEBUG
printf("row %d, old parent %d, new parent %d, pivotrow %d\n",kRow,
iParent,newParent,pivotRow);
#endif
int i1 = permuteBack_[pivotRow];
int i2 = permuteBack_[kRow];
permuteBack_[pivotRow]=i2;
permuteBack_[kRow]=i1;
// do Btran permutation
permute_[i1]=kRow;
permute_[i2]=pivotRow;
pivotRow = kRow;
// Take out of old parent
iLeft = leftSibling_[kRow];
iRight = rightSibling_[kRow];
if (iLeft<0) {
// take out of tree
if (iRight>=0) {
leftSibling_[iRight]=iLeft;
descendant_[iParent]=iRight;
} else {
#ifdef FULL_DEBUG
printf("Saying %d (old parent of %d) has no descendants\n",
iParent, kRow);
#endif
descendant_[iParent]=-1;
}
} else {
// take out of tree
rightSibling_[iLeft] = iRight;
if (iRight>=0)
leftSibling_[iRight]=iLeft;
}
leftSibling_[kRow]=-1;
rightSibling_[kRow]=-1;
// now insert new one
// make this descendant of that
if (descendant_[newParent]>=0) {
// we will have a sibling
int jRight = descendant_[newParent];
rightSibling_[kRow]=jRight;
leftSibling_[jRight]=kRow;
} else {
rightSibling_[kRow]=-1;
}
descendant_[newParent] = kRow;
leftSibling_[kRow]=-1;
parent_[kRow]=newParent;
iParent = kRow;
}
// now redo all depths from stack_[1]
// This must be possible to combine - but later
{
int iPivot = stack_[1];
int iDepth=depth_[parent_[iPivot]]; //depth of parent
iDepth ++; //correct for this one
int nStack = 1;
stack_[0]=iPivot;
while (nStack) {
// take off
int iNext = stack_[--nStack];
if (iNext>=0) {
// add stack level
depth_[iNext]=nStack+iDepth;
stack_[nStack++] = rightSibling_[iNext];
if (descendant_[iNext]>=0)
stack_[nStack++] = descendant_[iNext];
}
}
}
if (extraPrint)
print();
//check();
return 0;
}
/* Updates one column (FTRAN) from region2 */
double
ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2,
int pivotRow)
{
regionSparse->clear ( );
double *region = regionSparse->denseVector ( );
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero = regionSparse2->getNumElements ( );
int *regionIndex = regionSparse->getIndices ( );
int i;
bool doTwo = (numberNonZero==2);
int i0=-1;
int i1=-1;
if (doTwo) {
i0 = regionIndex2[0];
i1 = regionIndex2[1];
}
double returnValue=0.0;
bool packed = regionSparse2->packedMode();
if (packed) {
if (doTwo&®ion2[0]*region2[1]<0.0) {
region[i0] = region2[0];
region2[0]=0.0;
region[i1] = region2[1];
region2[1]=0.0;
int iDepth0 = depth_[i0];
int iDepth1 = depth_[i1];
if (iDepth1>iDepth0) {
int temp = i0;
i0 = i1;
i1 = temp;
temp = iDepth0;
iDepth0 = iDepth1;
iDepth1 = temp;
}
numberNonZero=0;
if (pivotRow<0) {
while (iDepth0>iDepth1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
region2[numberNonZero] = pivotValue*sign_[i0];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
iDepth0--;
i0 = otherRow;
}
while (i0!=i1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
region2[numberNonZero] = pivotValue*sign_[i0];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
i0 = otherRow;
double pivotValue1 = region[i1];
// put back now ?
int iBack1 = permuteBack_[i1];
region2[numberNonZero] = pivotValue1*sign_[i1];
regionIndex2[numberNonZero++] = iBack1;
int otherRow1 = parent_[i1];
region[i1] =0.0;
region[otherRow1] += pivotValue1;
i1 = otherRow1;
}
} else {
while (iDepth0>iDepth1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
double value = pivotValue*sign_[i0];
region2[numberNonZero] = value;
regionIndex2[numberNonZero++] = iBack;
if (iBack==pivotRow)
returnValue = value;
int otherRow = parent_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
iDepth0--;
i0 = otherRow;
}
while (i0!=i1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
double value = pivotValue*sign_[i0];
region2[numberNonZero] = value;
regionIndex2[numberNonZero++] = iBack;
if (iBack==pivotRow)
returnValue = value;
int otherRow = parent_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
i0 = otherRow;
double pivotValue1 = region[i1];
// put back now ?
int iBack1 = permuteBack_[i1];
value = pivotValue1*sign_[i1];
region2[numberNonZero] = value;
regionIndex2[numberNonZero++] = iBack1;
if (iBack1==pivotRow)
returnValue = value;
int otherRow1 = parent_[i1];
region[i1] =0.0;
region[otherRow1] += pivotValue1;
i1 = otherRow1;
}
}
} else {
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
//mark_[numberRows_]=1;
for (i=0;i<numberNonZero;i++) {
int j = regionIndex2[i];
double value = region2[i];
region2[i]=0.0;
region[j]=value;
regionIndex[i]=j;
int iDepth = depth_[j];
if (iDepth>greatestDepth)
greatestDepth = iDepth;
// and back until marked
while (!mark_[j]) {
int iNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=iNext;
mark_[j]=1;
iDepth--;
j=parent_[j];
}
}
numberNonZero=0;
if (pivotRow<0) {
for (;greatestDepth>=0; greatestDepth--) {
int iPivot = stack2_[greatestDepth];
stack2_[greatestDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region[iPivot];
if (pivotValue) {
// put back now ?
int iBack = permuteBack_[iPivot];
region2[numberNonZero] = pivotValue*sign_[iPivot];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[iPivot];
region[iPivot] =0.0;
region[otherRow] += pivotValue;
}
iPivot = stack_[iPivot];
}
}
} else {
for (;greatestDepth>=0; greatestDepth--) {
int iPivot = stack2_[greatestDepth];
stack2_[greatestDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region[iPivot];
if (pivotValue) {
// put back now ?
int iBack = permuteBack_[iPivot];
double value = pivotValue*sign_[iPivot];
region2[numberNonZero] = value;
regionIndex2[numberNonZero++] = iBack;
if (iBack==pivotRow)
returnValue = value;
int otherRow = parent_[iPivot];
region[iPivot] =0.0;
region[otherRow] += pivotValue;
}
iPivot = stack_[iPivot];
}
}
}
}
} else {
if (doTwo&®ion2[i0]*region2[i1]<0.0) {
// If just +- 1 then could go backwards on depth until join
region[i0] = region2[i0];
region2[i0]=0.0;
region[i1] = region2[i1];
region2[i1]=0.0;
int iDepth0 = depth_[i0];
int iDepth1 = depth_[i1];
if (iDepth1>iDepth0) {
int temp = i0;
i0 = i1;
i1 = temp;
temp = iDepth0;
iDepth0 = iDepth1;
iDepth1 = temp;
}
numberNonZero=0;
while (iDepth0>iDepth1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[i0];
region2[iBack] = pivotValue*sign_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
iDepth0--;
i0 = otherRow;
}
while (i0!=i1) {
double pivotValue = region[i0];
// put back now ?
int iBack = permuteBack_[i0];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[i0];
region2[iBack] = pivotValue*sign_[i0];
region[i0] =0.0;
region[otherRow] += pivotValue;
i0 = otherRow;
double pivotValue1 = region[i1];
// put back now ?
int iBack1 = permuteBack_[i1];
regionIndex2[numberNonZero++] = iBack1;
int otherRow1 = parent_[i1];
region2[iBack1] = pivotValue1*sign_[i1];
region[i1] =0.0;
region[otherRow1] += pivotValue1;
i1 = otherRow1;
}
} else {
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
//mark_[numberRows_]=1;
for (i=0;i<numberNonZero;i++) {
int j = regionIndex2[i];
double value = region2[j];
region2[j]=0.0;
region[j]=value;
regionIndex[i]=j;
int iDepth = depth_[j];
if (iDepth>greatestDepth)
greatestDepth = iDepth;
// and back until marked
while (!mark_[j]) {
int iNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=iNext;
mark_[j]=1;
iDepth--;
j=parent_[j];
}
}
numberNonZero=0;
for (;greatestDepth>=0; greatestDepth--) {
int iPivot = stack2_[greatestDepth];
stack2_[greatestDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region[iPivot];
if (pivotValue) {
// put back now ?
int iBack = permuteBack_[iPivot];
regionIndex2[numberNonZero++] = iBack;
int otherRow = parent_[iPivot];
region2[iBack] = pivotValue*sign_[iPivot];
region[iPivot] =0.0;
region[otherRow] += pivotValue;
}
iPivot = stack_[iPivot];
}
}
}
if (pivotRow>=0)
returnValue = region2[pivotRow];
}
region[numberRows_]=0.0;
regionSparse2->setNumElements(numberNonZero);
#ifdef FULL_DEBUG
{
int i;
for (i=0;i<numberRows_;i++) {
assert(!mark_[i]);
assert (stack2_[i]==-1);
}
}
#endif
return returnValue;
}
/* Updates one column (FTRAN) to/from array
** For large problems you should ALWAYS know where the nonzeros
are, so please try and migrate to previous method after you
have got code working using this simple method - thank you!
(the only exception is if you know input is dense e.g. rhs) */
int
ClpNetworkBasis::updateColumn ( CoinIndexedVector * regionSparse,
double region2[] ) const
{
regionSparse->clear ( );
double *region = regionSparse->denseVector ( );
int numberNonZero = 0;
int *regionIndex = regionSparse->getIndices ( );
int i;
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
for (i=0;i<numberRows_;i++) {
double value = region2[i];
if (value) {
region2[i]=0.0;
region[i]=value;
regionIndex[numberNonZero++]=i;
int j=i;
int iDepth = depth_[j];
if (iDepth>greatestDepth)
greatestDepth = iDepth;
// and back until marked
while (!mark_[j]) {
int iNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=iNext;
mark_[j]=1;
iDepth--;
j=parent_[j];
}
}
}
numberNonZero=0;
for (;greatestDepth>=0; greatestDepth--) {
int iPivot = stack2_[greatestDepth];
stack2_[greatestDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region[iPivot];
if (pivotValue) {
// put back now ?
int iBack = permuteBack_[iPivot];
numberNonZero++;
int otherRow = parent_[iPivot];
region2[iBack] = pivotValue*sign_[iPivot];
region[iPivot] =0.0;
region[otherRow] += pivotValue;
}
iPivot = stack_[iPivot];
}
}
region[numberRows_]=0.0;
return numberNonZero;
}
/* Updates one column transpose (BTRAN)
For large problems you should ALWAYS know where the nonzeros
are, so please try and migrate to previous method after you
have got code working using this simple method - thank you!
(the only exception is if you know input is dense e.g. dense objective)
returns number of nonzeros */
int
ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
double region2[] ) const
{
// permute in after copying
// so will end up in right place
double *region = regionSparse->denseVector ( );
int *regionIndex = regionSparse->getIndices ( );
int i;
int numberNonZero=0;
memcpy(region,region2,numberRows_*sizeof(double));
for (i=0;i<numberRows_;i++) {
double value = region[i];
if (value) {
int k = permute_[i];
region[i]=0.0;
region2[k]=value;
regionIndex[numberNonZero++]=k;
mark_[k]=1;
}
}
// copy back
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
int smallestDepth=numberRows_;
for (i=0;i<numberNonZero;i++) {
int j = regionIndex[i];
// add in
int iDepth = depth_[j];
smallestDepth = CoinMin(iDepth,smallestDepth) ;
greatestDepth = CoinMax(iDepth,greatestDepth) ;
int jNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=jNext;
// and put all descendants on list
int iChild = descendant_[j];
while (iChild>=0) {
if (!mark_[iChild]) {
regionIndex[numberNonZero++] = iChild;
mark_[iChild]=1;
}
iChild = rightSibling_[iChild];
}
}
numberNonZero=0;
region2[numberRows_]=0.0;
int iDepth;
for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
int iPivot = stack2_[iDepth];
stack2_[iDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region2[iPivot];
int otherRow = parent_[iPivot];
double otherValue = region2[otherRow];
pivotValue = sign_[iPivot]*pivotValue+otherValue;
region2[iPivot]=pivotValue;
if (pivotValue)
numberNonZero++;
iPivot = stack_[iPivot];
}
}
return numberNonZero;
}
/* Updates one column (BTRAN) from region2 */
int
ClpNetworkBasis::updateColumnTranspose ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2) const
{
// permute in - presume small number so copy back
// so will end up in right place
regionSparse->clear ( );
double *region = regionSparse->denseVector ( );
double *region2 = regionSparse2->denseVector ( );
int *regionIndex2 = regionSparse2->getIndices ( );
int numberNonZero2 = regionSparse2->getNumElements ( );
int *regionIndex = regionSparse->getIndices ( );
int i;
int numberNonZero=0;
bool packed = regionSparse2->packedMode();
if (packed) {
for (i=0;i<numberNonZero2;i++) {
int k = regionIndex2[i];
int j = permute_[k];
double value = region2[i];
region2[i]=0.0;
region[j]=value;
mark_[j]=1;
regionIndex[numberNonZero++]=j;
}
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
int smallestDepth=numberRows_;
//mark_[numberRows_]=1;
for (i=0;i<numberNonZero2;i++) {
int j = regionIndex[i];
regionIndex2[i]=j;
// add in
int iDepth = depth_[j];
smallestDepth = CoinMin(iDepth,smallestDepth) ;
greatestDepth = CoinMax(iDepth,greatestDepth) ;
int jNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=jNext;
// and put all descendants on list
int iChild = descendant_[j];
while (iChild>=0) {
if (!mark_[iChild]) {
regionIndex2[numberNonZero++] = iChild;
mark_[iChild]=1;
}
iChild = rightSibling_[iChild];
}
}
for (;i<numberNonZero;i++) {
int j = regionIndex2[i];
// add in
int iDepth = depth_[j];
smallestDepth = CoinMin(iDepth,smallestDepth) ;
greatestDepth = CoinMax(iDepth,greatestDepth) ;
int jNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=jNext;
// and put all descendants on list
int iChild = descendant_[j];
while (iChild>=0) {
if (!mark_[iChild]) {
regionIndex2[numberNonZero++] = iChild;
mark_[iChild]=1;
}
iChild = rightSibling_[iChild];
}
}
numberNonZero2=0;
region[numberRows_]=0.0;
int iDepth;
for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
int iPivot = stack2_[iDepth];
stack2_[iDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region[iPivot];
int otherRow = parent_[iPivot];
double otherValue = region[otherRow];
pivotValue = sign_[iPivot]*pivotValue+otherValue;
region[iPivot]=pivotValue;
if (pivotValue) {
region2[numberNonZero2]=pivotValue;
regionIndex2[numberNonZero2++]=iPivot;
}
iPivot = stack_[iPivot];
}
}
// zero out
for (i=0;i<numberNonZero2;i++) {
int k = regionIndex2[i];
region[k]=0.0;
}
} else {
for (i=0;i<numberNonZero2;i++) {
int k = regionIndex2[i];
int j = permute_[k];
double value = region2[k];
region2[k]=0.0;
region[j]=value;
mark_[j]=1;
regionIndex[numberNonZero++]=j;
}
// copy back
// set up linked lists at each depth
// stack2 is start, stack is next
int greatestDepth=-1;
int smallestDepth=numberRows_;
//mark_[numberRows_]=1;
for (i=0;i<numberNonZero2;i++) {
int j = regionIndex[i];
double value = region[j];
region[j]=0.0;
region2[j]=value;
regionIndex2[i]=j;
// add in
int iDepth = depth_[j];
smallestDepth = CoinMin(iDepth,smallestDepth) ;
greatestDepth = CoinMax(iDepth,greatestDepth) ;
int jNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=jNext;
// and put all descendants on list
int iChild = descendant_[j];
while (iChild>=0) {
if (!mark_[iChild]) {
regionIndex2[numberNonZero++] = iChild;
mark_[iChild]=1;
}
iChild = rightSibling_[iChild];
}
}
for (;i<numberNonZero;i++) {
int j = regionIndex2[i];
// add in
int iDepth = depth_[j];
smallestDepth = CoinMin(iDepth,smallestDepth) ;
greatestDepth = CoinMax(iDepth,greatestDepth) ;
int jNext = stack2_[iDepth];
stack2_[iDepth]=j;
stack_[j]=jNext;
// and put all descendants on list
int iChild = descendant_[j];
while (iChild>=0) {
if (!mark_[iChild]) {
regionIndex2[numberNonZero++] = iChild;
mark_[iChild]=1;
}
iChild = rightSibling_[iChild];
}
}
numberNonZero2=0;
region2[numberRows_]=0.0;
int iDepth;
for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
int iPivot = stack2_[iDepth];
stack2_[iDepth]=-1;
while (iPivot>=0) {
mark_[iPivot]=0;
double pivotValue = region2[iPivot];
int otherRow = parent_[iPivot];
double otherValue = region2[otherRow];
pivotValue = sign_[iPivot]*pivotValue+otherValue;
region2[iPivot]=pivotValue;
if (pivotValue)
regionIndex2[numberNonZero2++]=iPivot;
iPivot = stack_[iPivot];
}
}
}
regionSparse2->setNumElements(numberNonZero2);
#ifdef FULL_DEBUG
{
int i;
for (i=0;i<numberRows_;i++) {
assert(!mark_[i]);
assert (stack2_[i]==-1);
}
}
#endif
return numberNonZero2;
}
syntax highlighted by Code2HTML, v. 0.9.1