// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include "CoinPragma.hpp"
#include "ClpFactorization.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif
#include "CoinHelperFunctions.hpp"
#include "CoinIndexedVector.hpp"
#include "ClpSimplex.hpp"
#include "ClpMatrixBase.hpp"
#ifndef SLIM_CLP
#include "ClpNetworkBasis.hpp"
#include "ClpNetworkMatrix.hpp"
//#define CHECK_NETWORK
#ifdef CHECK_NETWORK
const static bool doCheck=true;
#else
const static bool doCheck=false;
#endif
#endif
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization () :
CoinFactorization()
{
#ifndef SLIM_CLP
networkBasis_ = NULL;
#endif
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpFactorization::ClpFactorization (const ClpFactorization & rhs) :
CoinFactorization(rhs)
{
#ifndef SLIM_CLP
if (rhs.networkBasis_)
networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
else
networkBasis_=NULL;
#endif
}
ClpFactorization::ClpFactorization (const CoinFactorization & rhs) :
CoinFactorization(rhs)
{
#ifndef SLIM_CLP
networkBasis_=NULL;
#endif
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpFactorization::~ClpFactorization ()
{
#ifndef SLIM_CLP
delete networkBasis_;
#endif
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpFactorization &
ClpFactorization::operator=(const ClpFactorization& rhs)
{
if (this != &rhs) {
CoinFactorization::operator=(rhs);
#ifndef SLIM_CLP
delete networkBasis_;
if (rhs.networkBasis_)
networkBasis_ = new ClpNetworkBasis(*(rhs.networkBasis_));
else
networkBasis_=NULL;
#endif
}
return *this;
}
#if 0
static unsigned int saveList[10000];
int numberSave=-1;
inline bool isDense(int i) {
return ((saveList[i>>5]>>(i&31))&1)!=0;
}
inline void setDense(int i) {
unsigned int & value = saveList[i>>5];
int bit = i&31;
value |= (1<<bit);
}
#endif
int
ClpFactorization::factorize ( ClpSimplex * model,
int solveType, bool valuesPass)
{
ClpMatrixBase * matrix = model->clpMatrix();
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
if (!numberRows)
return 0;
// If too many compressions increase area
if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
areaFactor_ *= 1.1;
}
//int numberPivots=numberPivots_;
#if 0
if (model->algorithm()>0)
numberSave=-1;
#endif
#ifndef SLIM_CLP
if (!networkBasis_||doCheck) {
#endif
status_=-99;
int * pivotVariable = model->pivotVariable();
//returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
while (status_<-98) {
int i;
int numberBasic=0;
int numberRowBasic;
// Move pivot variables across if they look good
int * pivotTemp = model->rowArray(0)->getIndices();
assert (!model->rowArray(0)->getNumElements());
if (!matrix->rhsOffset(model)) {
#if 0
if (numberSave>0) {
int nStill=0;
int nAtBound=0;
int nZeroDual=0;
CoinIndexedVector * array = model->rowArray(3);
CoinIndexedVector * objArray = model->columnArray(1);
array->clear();
objArray->clear();
double * cost = model->costRegion();
double tolerance = model->primalTolerance();
double offset=0.0;
for (i=0;i<numberRows;i++) {
int iPivot = pivotVariable[i];
if (iPivot<numberColumns&&isDense(iPivot)) {
if (model->getColumnStatus(iPivot)==ClpSimplex::basic) {
nStill++;
double value=model->solutionRegion()[iPivot];
double dual = model->dualRowSolution()[i];
double lower=model->lowerRegion()[iPivot];
double upper=model->upperRegion()[iPivot];
ClpSimplex::Status status;
if (fabs(value-lower)<tolerance) {
status=ClpSimplex::atLowerBound;
nAtBound++;
} else if (fabs(value-upper)<tolerance) {
nAtBound++;
status=ClpSimplex::atUpperBound;
} else if (value>lower&&value<upper) {
status=ClpSimplex::superBasic;
} else {
status=ClpSimplex::basic;
}
if (status!=ClpSimplex::basic) {
if (model->getRowStatus(i)!=ClpSimplex::basic) {
model->setColumnStatus(iPivot,ClpSimplex::atLowerBound);
model->setRowStatus(i,ClpSimplex::basic);
pivotVariable[i]=i+numberColumns;
model->dualRowSolution()[i]=0.0;
model->djRegion(0)[i]=0.0;
array->add(i,dual);
offset += dual*model->solutionRegion(0)[i];
}
}
if (fabs(dual)<1.0e-5)
nZeroDual++;
}
}
}
printf("out of %d dense, %d still in basis, %d at bound, %d with zero dual - offset %g\n",
numberSave,nStill,nAtBound,nZeroDual,offset);
if (array->getNumElements()) {
// modify costs
model->clpMatrix()->transposeTimes(model,1.0,array,model->columnArray(0),
objArray);
array->clear();
int n=objArray->getNumElements();
int * indices = objArray->getIndices();
double * elements = objArray->denseVector();
for (i=0;i<n;i++) {
int iColumn = indices[i];
cost[iColumn] -= elements[iColumn];
elements[iColumn]=0.0;
}
objArray->setNumElements(0);
}
}
#endif
// Seems to prefer things in order so quickest
// way is to go though like this
for (i=0;i<numberRows;i++) {
if (model->getRowStatus(i) == ClpSimplex::basic)
pivotTemp[numberBasic++]=i;
}
numberRowBasic=numberBasic;
/* Put column basic variables into pivotVariable
This is done by ClpMatrixBase to allow override for gub
*/
matrix->generalExpanded(model,0,numberBasic);
} else {
// Long matrix - do a different way
bool fullSearch=false;
for (i=0;i<numberRows;i++) {
int iPivot = pivotVariable[i];
if (iPivot>=numberColumns) {
pivotTemp[numberBasic++]=iPivot-numberColumns;
}
}
numberRowBasic=numberBasic;
for (i=0;i<numberRows;i++) {
int iPivot = pivotVariable[i];
if (iPivot<numberColumns) {
if (iPivot>=0) {
pivotTemp[numberBasic++]=iPivot;
} else {
// not full basis
fullSearch=true;
break;
}
}
}
if (fullSearch) {
// do slow way
numberBasic=0;
for (i=0;i<numberRows;i++) {
if (model->getRowStatus(i) == ClpSimplex::basic)
pivotTemp[numberBasic++]=i;
}
numberRowBasic=numberBasic;
/* Put column basic variables into pivotVariable
This is done by ClpMatrixBase to allow override for gub
*/
matrix->generalExpanded(model,0,numberBasic);
}
}
if (numberBasic>model->maximumBasic()) {
#if 0 // ndef NDEBUG
printf("%d basic - should only be %d\n",
numberBasic,numberRows);
#endif
// Take out some
numberBasic=numberRowBasic;
for (int i=0;i<numberColumns;i++) {
if (model->getColumnStatus(i) == ClpSimplex::basic) {
if (numberBasic<numberRows)
numberBasic++;
else
model->setColumnStatus(i,ClpSimplex::superBasic);
}
}
numberBasic=numberRowBasic;
matrix->generalExpanded(model,0,numberBasic);
}
#ifndef SLIM_CLP
// see if matrix a network
#ifndef NO_RTTI
ClpNetworkMatrix* networkMatrix =
dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
#else
ClpNetworkMatrix* networkMatrix = NULL;
if (model->clpMatrix()->type()==11)
networkMatrix =
static_cast< ClpNetworkMatrix*>(model->clpMatrix());
#endif
// If network - still allow ordinary factorization first time for laziness
if (networkMatrix)
biasLU_=0; // All to U if network
//int saveMaximumPivots = maximumPivots();
delete networkBasis_;
networkBasis_ = NULL;
if (networkMatrix&&!doCheck)
maximumPivots(1);
#endif
//printf("L, U, R %d %d %d\n",numberElementsL(),numberElementsU(),numberElementsR());
while (status_==-99) {
// maybe for speed will be better to leave as many regions as possible
gutsOfDestructor();
gutsOfInitialize(2);
CoinBigIndex numberElements=numberRowBasic;
// compute how much in basis
int i;
// can change for gub
int numberColumnBasic = numberBasic-numberRowBasic;
numberElements +=matrix->countBasis(model,
pivotTemp+numberRowBasic,
numberRowBasic,
numberColumnBasic);
// and recompute as network side say different
if (model->numberIterations())
numberRowBasic = numberBasic - numberColumnBasic;
numberElements = 3 * numberBasic + 3 * numberElements + 10000;
#if 0
// If iteration not zero then may be compressed
getAreas ( !model->numberIterations() ? numberRows : numberBasic,
numberRowBasic+numberColumnBasic, numberElements,
2 * numberElements );
#else
getAreas ( numberRows,
numberRowBasic+numberColumnBasic, numberElements,
2 * numberElements );
#endif
//fill
// Fill in counts so we can skip part of preProcess
int * numberInRow = numberInRow_.array();
int * numberInColumn = numberInColumn_.array();
CoinZeroN ( numberInRow, numberRows_ + 1 );
CoinZeroN ( numberInColumn, maximumColumnsExtra_ + 1 );
double * elementU = elementU_.array();
int * indexRowU = indexRowU_.array();
CoinBigIndex * startColumnU = startColumnU_.array();
for (i=0;i<numberRowBasic;i++) {
int iRow = pivotTemp[i];
indexRowU[i]=iRow;
startColumnU[i]=i;
elementU[i]=slackValue_;
numberInRow[iRow]=1;
numberInColumn[i]=1;
}
startColumnU[numberRowBasic]=numberRowBasic;
// can change for gub so redo
numberColumnBasic = numberBasic-numberRowBasic;
matrix->fillBasis(model,
pivotTemp+numberRowBasic,
numberColumnBasic,
indexRowU_.array(),
startColumnU+numberRowBasic,
numberInRow,
numberInColumn+numberRowBasic,
elementU_.array());
#if 0
{
printf("%d row basic, %d column basic\n",numberRowBasic,numberColumnBasic);
for (int i=0;i<numberElements;i++)
printf("row %d col %d value %g\n",indexRowU_.array()[i],indexColumnU_[i],
elementU_.array()[i]);
}
#endif
// recompute number basic
numberBasic = numberRowBasic+numberColumnBasic;
if (numberBasic)
numberElements = startColumnU[numberBasic-1]
+numberInColumn[numberBasic-1];
else
numberElements=0;
lengthU_ = numberElements;
//saveFactorization("dump.d");
if (biasLU_>=3||numberRows_!=numberColumns_)
preProcess ( 2 );
else
preProcess ( 3 ); // no row copy
factor ( );
if (status_==-99) {
// get more memory
areaFactor(2.0*areaFactor());
} else if (status_==-1&&model->numberIterations()==0&&
denseThreshold_) {
// Round again without dense
denseThreshold_=0;
status_=-99;
}
}
// If we get here status is 0 or -1
if (status_ == 0) {
// We may need to tamper with order and redo - e.g. network with side
int useNumberRows = numberRows;
// **** we will also need to add test in dual steepest to do
// as we do for network
matrix->generalExpanded(model,12,useNumberRows);
const int * permuteBack = permuteBack_.array();
const int * back = pivotColumnBack_.array();
//int * pivotTemp = pivotColumn_.array();
//ClpDisjointCopyN ( pivotVariable, numberRows , pivotTemp );
// Redo pivot order
for (i=0;i<numberRowBasic;i++) {
int k = pivotTemp[i];
// so rowIsBasic[k] would be permuteBack[back[i]]
pivotVariable[permuteBack[back[i]]]=k+numberColumns;
}
for (;i<useNumberRows;i++) {
int k = pivotTemp[i];
// so rowIsBasic[k] would be permuteBack[back[i]]
pivotVariable[permuteBack[back[i]]]=k;
}
#if 0
if (numberSave>=0) {
numberSave=numberDense_;
memset(saveList,0,((numberRows_+31)>>5)*sizeof(int));
for (i=numberRows_-numberSave;i<numberRows_;i++) {
int k=pivotTemp[pivotColumn_.array()[i]];
setDense(k);
}
}
#endif
// Set up permutation vector
// these arrays start off as copies of permute
// (and we could use permute_ instead of pivotColumn (not back though))
ClpDisjointCopyN ( permute_.array(), useNumberRows , pivotColumn_.array() );
ClpDisjointCopyN ( permuteBack_.array(), useNumberRows , pivotColumnBack_.array() );
#ifndef SLIM_CLP
if (networkMatrix) {
maximumPivots(CoinMax(2000,maximumPivots()));
// redo arrays
for (int iRow=0;iRow<4;iRow++) {
int length =model->numberRows()+maximumPivots();
if (iRow==3||model->objectiveAsObject()->type()>1)
length += model->numberColumns();
model->rowArray(iRow)->reserve(length);
}
// create network factorization
if (doCheck)
delete networkBasis_; // temp
networkBasis_ = new ClpNetworkBasis(model,numberRows_,
pivotRegion_.array(),
permuteBack_.array(),
startColumnU_.array(),
numberInColumn_.array(),
indexRowU_.array(),
elementU_.array());
// kill off arrays in ordinary factorization
if (!doCheck) {
gutsOfDestructor();
// but make sure numberRows_ set
numberRows_ = model->numberRows();
status_=0;
#if 0
// but put back permute arrays so odd things will work
int numberRows = model->numberRows();
pivotColumnBack_ = new int [numberRows];
permute_ = new int [numberRows];
int i;
for (i=0;i<numberRows;i++) {
pivotColumnBack_[i]=i;
permute_[i]=i;
}
#endif
}
} else {
#endif
// See if worth going sparse and when
if (numberFtranCounts_>100) {
ftranCountInput_= CoinMax(ftranCountInput_,1.0);
ftranAverageAfterL_ = CoinMax(ftranCountAfterL_/ftranCountInput_,1.0);
ftranAverageAfterR_ = CoinMax(ftranCountAfterR_/ftranCountAfterL_,1.0);
ftranAverageAfterU_ = CoinMax(ftranCountAfterU_/ftranCountAfterR_,1.0);
if (btranCountInput_&&btranCountAfterU_&&btranCountAfterR_) {
btranAverageAfterU_ = CoinMax(btranCountAfterU_/btranCountInput_,1.0);
btranAverageAfterR_ = CoinMax(btranCountAfterR_/btranCountAfterU_,1.0);
btranAverageAfterL_ = CoinMax(btranCountAfterL_/btranCountAfterR_,1.0);
} else {
// we have not done any useful btrans (values pass?)
btranAverageAfterU_ = 1.0;
btranAverageAfterR_ = 1.0;
btranAverageAfterL_ = 1.0;
}
}
// scale back
ftranCountInput_ *= 0.8;
ftranCountAfterL_ *= 0.8;
ftranCountAfterR_ *= 0.8;
ftranCountAfterU_ *= 0.8;
btranCountInput_ *= 0.8;
btranCountAfterU_ *= 0.8;
btranCountAfterR_ *= 0.8;
btranCountAfterL_ *= 0.8;
#ifndef SLIM_CLP
}
#endif
} else if (status_ == -1&&(solveType==0||solveType==2)) {
// This needs redoing as it was merged coding - does not need array
int numberTotal = numberRows+numberColumns;
int * isBasic = new int [numberTotal];
int * rowIsBasic = isBasic+numberColumns;
int * columnIsBasic = isBasic;
for (i=0;i<numberTotal;i++)
isBasic[i]=-1;
for (i=0;i<numberRowBasic;i++) {
int iRow = pivotTemp[i];
rowIsBasic[iRow]=1;
}
for (;i<numberBasic;i++) {
int iColumn = pivotTemp[i];
columnIsBasic[iColumn]=1;
}
numberBasic=0;
for (i=0;i<numberRows;i++)
pivotVariable[i]=-1;
// mark as basic or non basic
const int * pivotColumn = pivotColumn_.array();
for (i=0;i<numberRows;i++) {
if (rowIsBasic[i]>=0) {
if (pivotColumn[numberBasic]>=0) {
rowIsBasic[i]=pivotColumn[numberBasic];
} else {
rowIsBasic[i]=-1;
model->setRowStatus(i,ClpSimplex::superBasic);
}
numberBasic++;
}
}
for (i=0;i<numberColumns;i++) {
if (columnIsBasic[i]>=0) {
if (pivotColumn[numberBasic]>=0)
columnIsBasic[i]=pivotColumn[numberBasic];
else
columnIsBasic[i]=-1;
numberBasic++;
}
}
// leave pivotVariable in useful form for cleaning basis
int * pivotVariable = model->pivotVariable();
for (i=0;i<numberRows;i++) {
pivotVariable[i]=-1;
}
for (i=0;i<numberRows;i++) {
if (model->getRowStatus(i) == ClpSimplex::basic) {
int iPivot = rowIsBasic[i];
if (iPivot>=0)
pivotVariable[iPivot]=i+numberColumns;
}
}
for (i=0;i<numberColumns;i++) {
if (model->getColumnStatus(i) == ClpSimplex::basic) {
int iPivot = columnIsBasic[i];
if (iPivot>=0)
pivotVariable[iPivot]=i;
}
}
delete [] isBasic;
double * columnLower = model->lowerRegion();
double * columnUpper = model->upperRegion();
double * columnActivity = model->solutionRegion();
double * rowLower = model->lowerRegion(0);
double * rowUpper = model->upperRegion(0);
double * rowActivity = model->solutionRegion(0);
//redo basis - first take ALL columns out
int iColumn;
double largeValue = model->largeValue();
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
// take out
if (!valuesPass) {
double lower=columnLower[iColumn];
double upper=columnUpper[iColumn];
double value=columnActivity[iColumn];
if (lower>-largeValue||upper<largeValue) {
if (fabs(value-lower)<fabs(value-upper)) {
model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
columnActivity[iColumn]=lower;
} else {
model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
columnActivity[iColumn]=upper;
}
} else {
model->setColumnStatus(iColumn,ClpSimplex::isFree);
}
} else {
model->setColumnStatus(iColumn,ClpSimplex::superBasic);
}
}
}
int iRow;
for (iRow=0;iRow<numberRows;iRow++) {
int iSequence=pivotVariable[iRow];
if (iSequence>=0) {
// basic
if (iSequence>=numberColumns) {
// slack in - leave
//assert (iSequence-numberColumns==iRow);
} else {
assert(model->getRowStatus(iRow)!=ClpSimplex::basic);
// put back structural
model->setColumnStatus(iSequence,ClpSimplex::basic);
}
} else {
// put in slack
model->setRowStatus(iRow,ClpSimplex::basic);
}
}
// Put back any key variables for gub (status_ not touched)
matrix->generalExpanded(model,1,status_);
// signal repeat
status_=-99;
// set fixed if they are
for (iRow=0;iRow<numberRows;iRow++) {
if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
if (rowLower[iRow]==rowUpper[iRow]) {
rowActivity[iRow]=rowLower[iRow];
model->setRowStatus(iRow,ClpSimplex::isFixed);
}
}
}
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
if (columnLower[iColumn]==columnUpper[iColumn]) {
columnActivity[iColumn]=columnLower[iColumn];
model->setColumnStatus(iColumn,ClpSimplex::isFixed);
}
}
}
}
}
#ifndef SLIM_CLP
} else {
// network - fake factorization - do nothing
status_=0;
numberPivots_ = 0;
}
#endif
#ifndef SLIM_CLP
if (!status_) {
// take out part if quadratic
if (model->algorithm()==2) {
ClpObjective * obj = model->objectiveAsObject();
ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
assert (quadraticObj);
CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
int numberXColumns = quadratic->getNumCols();
assert (numberXColumns<numberColumns);
int base = numberColumns-numberXColumns;
int * which = new int [numberXColumns];
int * pivotVariable = model->pivotVariable();
int * permute = pivotColumn();
int i;
int n=0;
for (i=0;i<numberRows;i++) {
int iSj = pivotVariable[i]-base;
if (iSj>=0&&iSj<numberXColumns)
which[n++]=permute[i];
}
if (n)
emptyRows(n,which);
delete [] which;
}
}
#endif
return status_;
}
/* Replaces one Column to basis,
returns 0=OK, 1=Probably OK, 2=singular, 3=no room
If checkBeforeModifying is true will do all accuracy checks
before modifying factorization. Whether to set this depends on
speed considerations. You could just do this on first iteration
after factorization and thereafter re-factorize
partial update already in U */
int
ClpFactorization::replaceColumn ( const ClpSimplex * model,
CoinIndexedVector * regionSparse,
CoinIndexedVector * tableauColumn,
int pivotRow,
double pivotCheck ,
bool checkBeforeModifying)
{
#ifndef SLIM_CLP
if (!networkBasis_) {
#endif
// see if FT
if (doForrestTomlin_) {
int returnCode= CoinFactorization::replaceColumn(regionSparse,
pivotRow,
pivotCheck,
checkBeforeModifying);
return returnCode;
} else {
return CoinFactorization::replaceColumnPFI(tableauColumn,
pivotRow,pivotCheck); // Note array
}
#ifndef SLIM_CLP
} else {
if (doCheck) {
int returnCode = CoinFactorization::replaceColumn(regionSparse,
pivotRow,
pivotCheck,
checkBeforeModifying);
networkBasis_->replaceColumn(regionSparse,
pivotRow);
return returnCode;
} else {
// increase number of pivots
numberPivots_++;
return networkBasis_->replaceColumn(regionSparse,
pivotRow);
}
}
#endif
}
/* Updates one column (FTRAN) from region2
number returned is negative if no room
region1 starts as zero and is zero at end */
int
ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2)
{
#ifdef CLP_DEBUG
regionSparse->checkClear();
#endif
if (!numberRows_)
return 0;
#ifndef SLIM_CLP
if (!networkBasis_) {
#endif
collectStatistics_ = true;
int returnValue= CoinFactorization::updateColumnFT(regionSparse,
regionSparse2);
collectStatistics_ = false;
return returnValue;
#ifndef SLIM_CLP
} else {
#ifdef CHECK_NETWORK
CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
double * check = new double[numberRows_];
int returnCode = CoinFactorization::updateColumnFT(regionSparse,
regionSparse2);
networkBasis_->updateColumn(regionSparse,save,-1);
int i;
double * array = regionSparse2->denseVector();
int * indices = regionSparse2->getIndices();
int n=regionSparse2->getNumElements();
memset(check,0,numberRows_*sizeof(double));
double * array2 = save->denseVector();
int * indices2 = save->getIndices();
int n2=save->getNumElements();
assert (n==n2);
if (save->packedMode()) {
for (i=0;i<n;i++) {
check[indices[i]]=array[i];
}
for (i=0;i<n;i++) {
double value2 = array2[i];
assert (check[indices2[i]]==value2);
}
} else {
for (i=0;i<numberRows_;i++) {
double value1 = array[i];
double value2 = array2[i];
assert (value1==value2);
}
}
delete save;
delete [] check;
return returnCode;
#else
networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
return 1;
#endif
}
#endif
}
/* Updates one column (FTRAN) from region2
number returned is negative if no room
region1 starts as zero and is zero at end */
int
ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2,
bool noPermute) const
{
#ifdef CLP_DEBUG
if (!noPermute)
regionSparse->checkClear();
#endif
if (!numberRows_)
return 0;
#ifndef SLIM_CLP
if (!networkBasis_) {
#endif
collectStatistics_ = true;
int returnValue= CoinFactorization::updateColumn(regionSparse,
regionSparse2,
noPermute);
collectStatistics_ = false;
return returnValue;
#ifndef SLIM_CLP
} else {
#ifdef CHECK_NETWORK
CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
double * check = new double[numberRows_];
int returnCode = CoinFactorization::updateColumn(regionSparse,
regionSparse2,
noPermute);
networkBasis_->updateColumn(regionSparse,save,-1);
int i;
double * array = regionSparse2->denseVector();
int * indices = regionSparse2->getIndices();
int n=regionSparse2->getNumElements();
memset(check,0,numberRows_*sizeof(double));
double * array2 = save->denseVector();
int * indices2 = save->getIndices();
int n2=save->getNumElements();
assert (n==n2);
if (save->packedMode()) {
for (i=0;i<n;i++) {
check[indices[i]]=array[i];
}
for (i=0;i<n;i++) {
double value2 = array2[i];
assert (check[indices2[i]]==value2);
}
} else {
for (i=0;i<numberRows_;i++) {
double value1 = array[i];
double value2 = array2[i];
assert (value1==value2);
}
}
delete save;
delete [] check;
return returnCode;
#else
networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
return 1;
#endif
}
#endif
}
/* Updates one column (FTRAN) from region2
number returned is negative if no room
region1 starts as zero and is zero at end */
int
ClpFactorization::updateColumnForDebug ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2,
bool noPermute) const
{
if (!noPermute)
regionSparse->checkClear();
if (!numberRows_)
return 0;
collectStatistics_ = false;
int returnValue= CoinFactorization::updateColumn(regionSparse,
regionSparse2,
noPermute);
return returnValue;
}
/* Updates one column (BTRAN) from region2
region1 starts as zero and is zero at end */
int
ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
CoinIndexedVector * regionSparse2) const
{
if (!numberRows_)
return 0;
#ifndef SLIM_CLP
if (!networkBasis_) {
#endif
collectStatistics_ = true;
return CoinFactorization::updateColumnTranspose(regionSparse,
regionSparse2);
collectStatistics_ = false;
#ifndef SLIM_CLP
} else {
#ifdef CHECK_NETWORK
CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
double * check = new double[numberRows_];
int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
regionSparse2);
networkBasis_->updateColumnTranspose(regionSparse,save);
int i;
double * array = regionSparse2->denseVector();
int * indices = regionSparse2->getIndices();
int n=regionSparse2->getNumElements();
memset(check,0,numberRows_*sizeof(double));
double * array2 = save->denseVector();
int * indices2 = save->getIndices();
int n2=save->getNumElements();
assert (n==n2);
if (save->packedMode()) {
for (i=0;i<n;i++) {
check[indices[i]]=array[i];
}
for (i=0;i<n;i++) {
double value2 = array2[i];
assert (check[indices2[i]]==value2);
}
} else {
for (i=0;i<numberRows_;i++) {
double value1 = array[i];
double value2 = array2[i];
assert (value1==value2);
}
}
delete save;
delete [] check;
return returnCode;
#else
return networkBasis_->updateColumnTranspose(regionSparse,regionSparse2);
#endif
}
#endif
}
/* makes a row copy of L for speed and to allow very sparse problems */
void
ClpFactorization::goSparse()
{
#ifndef SLIM_CLP
if (!networkBasis_)
#endif
CoinFactorization::goSparse();
}
// Cleans up i.e. gets rid of network basis
void
ClpFactorization::cleanUp()
{
#ifndef SLIM_CLP
delete networkBasis_;
networkBasis_=NULL;
#endif
resetStatistics();
}
/// Says whether to redo pivot order
bool
ClpFactorization::needToReorder() const
{
#ifdef CHECK_NETWORK
return true;
#endif
#ifndef SLIM_CLP
if (!networkBasis_)
#endif
return true;
#ifndef SLIM_CLP
else
return false;
#endif
}
// Get weighted row list
void
ClpFactorization::getWeights(int * weights) const
{
#ifndef SLIM_CLP
if (networkBasis_) {
// Network - just unit
for (int i=0;i<numberRows_;i++)
weights[i]=1;
return;
}
#endif
int * numberInRow = numberInRow_.array();
int * numberInColumn = numberInColumn_.array();
int * permuteBack = pivotColumnBack_.array();
int * indexRowU = indexRowU_.array();
const CoinBigIndex * startColumnU = startColumnU_.array();
const CoinBigIndex * startRowL = startRowL_.array();
if (!startRowL||!numberInRow_.array()) {
int * temp = new int[numberRows_];
memset(temp,0,numberRows_*sizeof(int));
int i;
for (i=0;i<numberRows_;i++) {
// one for pivot
temp[i]++;
CoinBigIndex j;
for (j=startColumnU[i];j<startColumnU[i]+numberInColumn[i];j++) {
int iRow=indexRowU[j];
temp[iRow]++;
}
}
CoinBigIndex * startColumnL = startColumnL_.array();
int * indexRowL = indexRowL_.array();
for (i=baseL_;i<baseL_+numberL_;i++) {
CoinBigIndex j;
for (j=startColumnL[i];j<startColumnL[i+1];j++) {
int iRow = indexRowL[j];
temp[iRow]++;
}
}
for (i=0;i<numberRows_;i++) {
int number = temp[i];
int iPermute = permuteBack[i];
weights[iPermute]=number;
}
delete [] temp;
} else {
int i;
for (i=0;i<numberRows_;i++) {
int number = startRowL[i+1]-startRowL[i]+numberInRow[i]+1;
//number = startRowL[i+1]-startRowL[i]+1;
//number = numberInRow[i]+1;
int iPermute = permuteBack[i];
weights[iPermute]=number;
}
}
}
syntax highlighted by Code2HTML, v. 0.9.1