// Copyright (C) 2003, International Business Machines
// Corporation and others. All Rights Reserved.
#include <cstdio>
#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinPackedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
// at end to get min/max!
#include "ClpPlusMinusOneMatrix.hpp"
#include "ClpMessage.hpp"
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix ()
: ClpMatrixBase()
{
setType(12);
matrix_ = NULL;
startPositive_ = NULL;
startNegative_ = NULL;
lengths_=NULL;
indices_=NULL;
numberRows_=0;
numberColumns_=0;
columnOrdered_=true;
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs)
: ClpMatrixBase(rhs)
{
matrix_ = NULL;
startPositive_ = NULL;
startNegative_ = NULL;
lengths_=NULL;
indices_=NULL;
numberRows_=rhs.numberRows_;
numberColumns_=rhs.numberColumns_;
columnOrdered_=rhs.columnOrdered_;
if (numberColumns_) {
CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
indices_ = new int [ numberElements];
memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
startPositive_ = new CoinBigIndex [ numberColumns_+1];
memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(CoinBigIndex));
startNegative_ = new CoinBigIndex [ numberColumns_];
memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(CoinBigIndex));
}
int numberRows = getNumRows();
if (rhs.rhsOffset_&&numberRows) {
rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
} else {
rhsOffset_=NULL;
}
}
// Constructor from arrays
ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(int numberRows, int numberColumns,
bool columnOrdered, const int * indices,
const CoinBigIndex * startPositive,
const CoinBigIndex * startNegative)
: ClpMatrixBase()
{
setType(12);
matrix_ = NULL;
lengths_=NULL;
numberRows_=numberRows;
numberColumns_=numberColumns;
columnOrdered_=columnOrdered;
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
int numberElements = startPositive[numberMajor];
startPositive_ = ClpCopyOfArray(startPositive,numberMajor+1);
startNegative_ = ClpCopyOfArray(startNegative,numberMajor);
indices_ = ClpCopyOfArray(indices,numberElements);
// Check valid
checkValid(false);
}
ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs)
: ClpMatrixBase()
{
setType(12);
matrix_ = NULL;
startPositive_ = NULL;
startNegative_ = NULL;
lengths_=NULL;
indices_=NULL;
int iColumn;
assert (rhs.isColOrdered());
// get matrix data pointers
const int * row = rhs.getIndices();
const CoinBigIndex * columnStart = rhs.getVectorStarts();
const int * columnLength = rhs.getVectorLengths();
const double * elementByColumn = rhs.getElements();
numberColumns_ = rhs.getNumCols();
numberRows_=-1;
indices_ = new int[rhs.getNumElements()];
startPositive_ = new CoinBigIndex [numberColumns_+1];
startNegative_ = new CoinBigIndex [numberColumns_];
int * temp = new int [rhs.getNumRows()];
CoinBigIndex j=0;
CoinBigIndex numberGoodP=0;
CoinBigIndex numberGoodM=0;
CoinBigIndex numberBad=0;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
CoinBigIndex k;
int iNeg=0;
startPositive_[iColumn]=j;
for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
k++) {
int iRow;
if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
iRow = row[k];
numberRows_ = CoinMax(numberRows_,iRow);
indices_[j++]=iRow;
numberGoodP++;
} else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
iRow = row[k];
numberRows_ = CoinMax(numberRows_,iRow);
temp[iNeg++]=iRow;
numberGoodM++;
} else {
numberBad++;
}
}
// move negative
startNegative_[iColumn]=j;
for (k=0;k<iNeg;k++) {
indices_[j++] = temp[k];
}
}
startPositive_[numberColumns_]=j;
delete [] temp;
if (numberBad) {
delete [] indices_;
indices_=NULL;
numberRows_=0;
numberColumns_=0;
delete [] startPositive_;
delete [] startNegative_;
// Put in statistics
startPositive_ = new CoinBigIndex [3];
startPositive_[0]=numberGoodP;
startPositive_[1]=numberGoodM;
startPositive_[2]=numberBad;
startNegative_ = NULL;
} else {
numberRows_ ++; // correct
// but number should be same as rhs
assert (numberRows_<=rhs.getNumRows());
numberRows_ = rhs.getNumRows();
columnOrdered_ = true;
}
// Check valid
if (!numberBad)
checkValid(false);
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
{
delete matrix_;
delete [] startPositive_;
delete [] startNegative_;
delete [] lengths_;
delete [] indices_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpPlusMinusOneMatrix &
ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
{
if (this != &rhs) {
ClpMatrixBase::operator=(rhs);
delete matrix_;
delete [] startPositive_;
delete [] startNegative_;
delete [] lengths_;
delete [] indices_;
matrix_ = NULL;
startPositive_ = NULL;
lengths_=NULL;
indices_=NULL;
numberRows_=rhs.numberRows_;
numberColumns_=rhs.numberColumns_;
columnOrdered_=rhs.columnOrdered_;
if (numberColumns_) {
CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
indices_ = new int [ numberElements];
memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
startPositive_ = new CoinBigIndex [ numberColumns_+1];
memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(CoinBigIndex));
startNegative_ = new CoinBigIndex [ numberColumns_];
memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(CoinBigIndex));
}
}
return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
{
return new ClpPlusMinusOneMatrix(*this);
}
/* Subset clone (without gaps). Duplicates are allowed
and order is as given */
ClpMatrixBase *
ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
int numberColumns,
const int * whichColumns) const
{
return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
numberColumns, whichColumns);
}
/* Subset constructor (without gaps). Duplicates are allowed
and order is as given */
ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
const ClpPlusMinusOneMatrix & rhs,
int numberRows, const int * whichRow,
int numberColumns, const int * whichColumn)
: ClpMatrixBase(rhs)
{
matrix_ = NULL;
startPositive_ = NULL;
startNegative_ = NULL;
lengths_=NULL;
indices_=NULL;
numberRows_=0;
numberColumns_=0;
columnOrdered_=rhs.columnOrdered_;
if (numberRows<=0||numberColumns<=0) {
startPositive_ = new CoinBigIndex [1];
startPositive_[0] = 0;
} else {
numberColumns_ = numberColumns;
numberRows_ = numberRows;
const int * index1 = rhs.indices_;
CoinBigIndex * startPositive1 = rhs.startPositive_;
int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
// Also swap incoming if not column ordered
if (!columnOrdered_) {
int temp1 = numberRows;
numberRows = numberColumns;
numberColumns = temp1;
const int * temp2;
temp2 = whichRow;
whichRow = whichColumn;
whichColumn = temp2;
}
// Throw exception if rhs empty
if (numberMajor1 <= 0 || numberMinor1 <= 0)
throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
// Array to say if an old row is in new copy
int * newRow = new int [numberMinor1];
int iRow;
for (iRow=0;iRow<numberMinor1;iRow++)
newRow[iRow] = -1;
// and array for duplicating rows
int * duplicateRow = new int [numberMinor];
int numberBad=0;
for (iRow=0;iRow<numberMinor;iRow++) {
duplicateRow[iRow] = -1;
int kRow = whichRow[iRow];
if (kRow>=0 && kRow < numberMinor1) {
if (newRow[kRow]<0) {
// first time
newRow[kRow]=iRow;
} else {
// duplicate
int lastRow = newRow[kRow];
newRow[kRow]=iRow;
duplicateRow[iRow] = lastRow;
}
} else {
// bad row
numberBad++;
}
}
if (numberBad)
throw CoinError("bad minor entries",
"subset constructor", "ClpPlusMinusOneMatrix");
// now get size and check columns
CoinBigIndex size = 0;
int iColumn;
numberBad=0;
for (iColumn=0;iColumn<numberMajor;iColumn++) {
int kColumn = whichColumn[iColumn];
if (kColumn>=0 && kColumn <numberMajor1) {
CoinBigIndex i;
for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
int kRow = index1[i];
kRow = newRow[kRow];
while (kRow>=0) {
size++;
kRow = duplicateRow[kRow];
}
}
} else {
// bad column
numberBad++;
printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
}
}
if (numberBad)
throw CoinError("bad major entries",
"subset constructor", "ClpPlusMinusOneMatrix");
// now create arrays
startPositive_ = new CoinBigIndex [numberMajor+1];
startNegative_ = new CoinBigIndex [numberMajor];
indices_ = new int[size];
// and fill them
size = 0;
startPositive_[0]=0;
CoinBigIndex * startNegative1 = rhs.startNegative_;
for (iColumn=0;iColumn<numberMajor;iColumn++) {
int kColumn = whichColumn[iColumn];
CoinBigIndex i;
for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
int kRow = index1[i];
kRow = newRow[kRow];
while (kRow>=0) {
indices_[size++] = kRow;
kRow = duplicateRow[kRow];
}
}
startNegative_[iColumn] = size;
for (;i<startPositive1[kColumn+1];i++) {
int kRow = index1[i];
kRow = newRow[kRow];
while (kRow>=0) {
indices_[size++] = kRow;
kRow = duplicateRow[kRow];
}
}
startPositive_[iColumn+1] = size;
}
delete [] newRow;
delete [] duplicateRow;
}
// Check valid
checkValid(false);
}
/* Returns a new matrix in reverse order without gaps */
ClpMatrixBase *
ClpPlusMinusOneMatrix::reverseOrderedCopy() const
{
int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
// count number in each row/column
CoinBigIndex * tempP = new CoinBigIndex [numberMinor];
CoinBigIndex * tempN = new CoinBigIndex [numberMinor];
memset(tempP,0,numberMinor*sizeof(CoinBigIndex));
memset(tempN,0,numberMinor*sizeof(CoinBigIndex));
CoinBigIndex j=0;
int i;
for (i=0;i<numberMajor;i++) {
for (;j<startNegative_[i];j++) {
int iRow = indices_[j];
tempP[iRow]++;
}
for (;j<startPositive_[i+1];j++) {
int iRow = indices_[j];
tempN[iRow]++;
}
}
int * newIndices = new int [startPositive_[numberMajor]];
CoinBigIndex * newP = new CoinBigIndex [numberMinor+1];
CoinBigIndex * newN = new CoinBigIndex[numberMinor];
int iRow;
j=0;
// do starts
for (iRow=0;iRow<numberMinor;iRow++) {
newP[iRow]=j;
j += tempP[iRow];
tempP[iRow]=newP[iRow];
newN[iRow] = j;
j += tempN[iRow];
tempN[iRow]=newN[iRow];
}
newP[numberMinor]=j;
j=0;
for (i=0;i<numberMajor;i++) {
for (;j<startNegative_[i];j++) {
int iRow = indices_[j];
CoinBigIndex put = tempP[iRow];
newIndices[put++] = i;
tempP[iRow] = put;
}
for (;j<startPositive_[i+1];j++) {
int iRow = indices_[j];
CoinBigIndex put = tempN[iRow];
newIndices[put++] = i;
tempN[iRow] = put;
}
}
delete [] tempP;
delete [] tempN;
ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
newCopy->passInCopy(numberMinor, numberMajor,
!columnOrdered_, newIndices, newP, newN);
return newCopy;
}
//unscaled versions
void
ClpPlusMinusOneMatrix::times(double scalar,
const double * x, double * y) const
{
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
int i;
CoinBigIndex j;
assert (columnOrdered_);
for (i=0;i<numberMajor;i++) {
double value = scalar*x[i];
if (value) {
for (j=startPositive_[i];j<startNegative_[i];j++) {
int iRow = indices_[j];
y[iRow] += value;
}
for (;j<startPositive_[i+1];j++) {
int iRow = indices_[j];
y[iRow] -= value;
}
}
}
}
void
ClpPlusMinusOneMatrix::transposeTimes(double scalar,
const double * x, double * y) const
{
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
int i;
CoinBigIndex j=0;
assert (columnOrdered_);
for (i=0;i<numberMajor;i++) {
double value = 0.0;
for (;j<startNegative_[i];j++) {
int iRow = indices_[j];
value += x[iRow];
}
for (;j<startPositive_[i+1];j++) {
int iRow = indices_[j];
value -= x[iRow];
}
y[i] += scalar*value;
}
}
void
ClpPlusMinusOneMatrix::times(double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale) const
{
// we know it is not scaled
times(scalar, x, y);
}
void
ClpPlusMinusOneMatrix::transposeTimes( double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale, double * spare) const
{
// we know it is not scaled
transposeTimes(scalar, x, y);
}
/* Return <code>x * A + y</code> in <code>z</code>.
Squashes small elements and knows about ClpSimplex */
void
ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
const CoinIndexedVector * rowArray,
CoinIndexedVector * y,
CoinIndexedVector * columnArray) const
{
// we know it is not scaled
columnArray->clear();
double * pi = rowArray->denseVector();
int numberNonZero=0;
int * index = columnArray->getIndices();
double * array = columnArray->denseVector();
int numberInRowArray = rowArray->getNumElements();
// maybe I need one in OsiSimplex
double zeroTolerance = model->factorization()->zeroTolerance();
int numberRows = model->numberRows();
bool packed = rowArray->packedMode();
#ifndef NO_RTTI
ClpPlusMinusOneMatrix* rowCopy =
dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
#else
ClpPlusMinusOneMatrix* rowCopy =
static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
#endif
double factor = 0.3;
// We may not want to do by row if there may be cache problems
int numberColumns = model->numberColumns();
// It would be nice to find L2 cache size - for moment 512K
// Be slightly optimistic
if (numberColumns*sizeof(double)>1000000) {
if (numberRows*10<numberColumns)
factor=0.1;
else if (numberRows*4<numberColumns)
factor=0.15;
else if (numberRows*2<numberColumns)
factor=0.2;
}
if (numberInRowArray>factor*numberRows||!rowCopy) {
assert (!y->getNumElements());
// do by column
// Need to expand if packed mode
int iColumn;
CoinBigIndex j=0;
assert (columnOrdered_);
if (packed) {
// need to expand pi into y
assert(y->capacity()>=numberRows);
double * piOld = pi;
pi = y->denseVector();
const int * whichRow = rowArray->getIndices();
int i;
// modify pi so can collapse to one loop
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=scalar*piOld[i];
}
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
double value = 0.0;
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
value += pi[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
value -= pi[iRow];
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=0.0;
}
} else {
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
double value = 0.0;
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
value += pi[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
value -= pi[iRow];
}
value *= scalar;
if (fabs(value)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn]=value;
}
}
}
columnArray->setNumElements(numberNonZero);
} else {
// do by row
rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
}
}
/* Return <code>x * A + y</code> in <code>z</code>.
Squashes small elements and knows about ClpSimplex */
void
ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
const CoinIndexedVector * rowArray,
CoinIndexedVector * y,
CoinIndexedVector * columnArray) const
{
columnArray->clear();
double * pi = rowArray->denseVector();
int numberNonZero=0;
int * index = columnArray->getIndices();
double * array = columnArray->denseVector();
int numberInRowArray = rowArray->getNumElements();
// maybe I need one in OsiSimplex
double zeroTolerance = model->factorization()->zeroTolerance();
const int * column = indices_;
const CoinBigIndex * startPositive = startPositive_;
const CoinBigIndex * startNegative = startNegative_;
const int * whichRow = rowArray->getIndices();
bool packed = rowArray->packedMode();
if (numberInRowArray>2) {
// do by rows
int iRow;
double * markVector = y->denseVector(); // probably empty .. but
int numberOriginal=0;
int i;
if (packed) {
numberNonZero=0;
// and set up mark as char array
char * marked = (char *) (index+columnArray->capacity());
double * array2 = y->denseVector();
#ifdef CLP_DEBUG
int numberColumns = model->numberColumns();
for (i=0;i<numberColumns;i++) {
assert(!marked[i]);
assert(!array2[i]);
}
#endif
for (i=0;i<numberInRowArray;i++) {
iRow = whichRow[i];
double value = pi[i]*scalar;
CoinBigIndex j;
for (j=startPositive[iRow];j<startNegative[iRow];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array2[iColumn] += value;
}
for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array2[iColumn] -= value;
}
}
// get rid of tiny values and zero out marked
numberOriginal=numberNonZero;
numberNonZero=0;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
if (marked[iColumn]) {
double value = array2[iColumn];
array2[iColumn]=0.0;
marked[iColumn]=0;
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
}
} else {
numberNonZero=0;
// and set up mark as char array
char * marked = (char *) markVector;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
marked[iColumn]=0;
}
for (i=0;i<numberInRowArray;i++) {
iRow = whichRow[i];
double value = pi[iRow]*scalar;
CoinBigIndex j;
for (j=startPositive[iRow];j<startNegative[iRow];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array[iColumn] += value;
}
for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array[iColumn] -= value;
}
}
// get rid of tiny values and zero out marked
numberOriginal=numberNonZero;
numberNonZero=0;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
marked[iColumn]=0;
if (fabs(array[iColumn])>zeroTolerance) {
index[numberNonZero++]=iColumn;
} else {
array[iColumn]=0.0;
}
}
}
} else if (numberInRowArray==2) {
/* do by rows when two rows (do longer first when not packed
and shorter first if packed */
int iRow0 = whichRow[0];
int iRow1 = whichRow[1];
CoinBigIndex j;
if (packed) {
double pi0 = pi[0];
double pi1 = pi[1];
if (startPositive[iRow0+1]-startPositive[iRow0]>
startPositive[iRow1+1]-startPositive[iRow1]) {
int temp = iRow0;
iRow0 = iRow1;
iRow1 = temp;
pi0=pi1;
pi1=pi[0];
}
// and set up mark as char array
char * marked = (char *) (index+columnArray->capacity());
int * lookup = y->getIndices();
double value = pi0*scalar;
for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
int iColumn = column[j];
array[numberNonZero] = value;
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
index[numberNonZero++]=iColumn;
}
for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
int iColumn = column[j];
array[numberNonZero] = -value;
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
index[numberNonZero++]=iColumn;
}
int numberOriginal = numberNonZero;
value = pi1*scalar;
for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
int iColumn = column[j];
if (marked[iColumn]) {
int iLookup = lookup[iColumn];
array[iLookup] += value;
} else {
if (fabs(value)>zeroTolerance) {
array[numberNonZero] = value;
index[numberNonZero++]=iColumn;
}
}
}
for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
int iColumn = column[j];
if (marked[iColumn]) {
int iLookup = lookup[iColumn];
array[iLookup] -= value;
} else {
if (fabs(value)>zeroTolerance) {
array[numberNonZero] = -value;
index[numberNonZero++]=iColumn;
}
}
}
// get rid of tiny values and zero out marked
int nDelete=0;
for (j=0;j<numberOriginal;j++) {
int iColumn = index[j];
marked[iColumn]=0;
if (fabs(array[j])<=zeroTolerance)
nDelete++;
}
if (nDelete) {
numberOriginal=numberNonZero;
numberNonZero=0;
for (j=0;j<numberOriginal;j++) {
int iColumn = index[j];
double value = array[j];
array[j]=0.0;
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
}
} else {
if (startPositive[iRow0+1]-startPositive[iRow0]<
startPositive[iRow1+1]-startPositive[iRow1]) {
int temp = iRow0;
iRow0 = iRow1;
iRow1 = temp;
}
int numberOriginal;
int i;
numberNonZero=0;
double value;
value = pi[iRow0]*scalar;
CoinBigIndex j;
for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
int iColumn = column[j];
index[numberNonZero++]=iColumn;
array[iColumn] = value;
}
for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
int iColumn = column[j];
index[numberNonZero++]=iColumn;
array[iColumn] = -value;
}
value = pi[iRow1]*scalar;
for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
int iColumn = column[j];
double value2= array[iColumn];
if (value2) {
value2 += value;
} else {
value2 = value;
index[numberNonZero++]=iColumn;
}
array[iColumn] = value2;
}
for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
int iColumn = column[j];
double value2= array[iColumn];
if (value2) {
value2 -= value;
} else {
value2 = -value;
index[numberNonZero++]=iColumn;
}
array[iColumn] = value2;
}
// get rid of tiny values and zero out marked
numberOriginal=numberNonZero;
numberNonZero=0;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
if (fabs(array[iColumn])>zeroTolerance) {
index[numberNonZero++]=iColumn;
} else {
array[iColumn]=0.0;
}
}
}
} else if (numberInRowArray==1) {
// Just one row
int iRow=rowArray->getIndices()[0];
numberNonZero=0;
double value;
iRow = whichRow[0];
CoinBigIndex j;
if (packed) {
value = pi[0]*scalar;
if (fabs(value)>zeroTolerance) {
for (j=startPositive[iRow];j<startNegative[iRow];j++) {
int iColumn = column[j];
array[numberNonZero] = value;
index[numberNonZero++]=iColumn;
}
for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
int iColumn = column[j];
array[numberNonZero] = -value;
index[numberNonZero++]=iColumn;
}
}
} else {
value = pi[iRow]*scalar;
if (fabs(value)>zeroTolerance) {
for (j=startPositive[iRow];j<startNegative[iRow];j++) {
int iColumn = column[j];
array[iColumn] = value;
index[numberNonZero++]=iColumn;
}
for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
int iColumn = column[j];
array[iColumn] = -value;
index[numberNonZero++]=iColumn;
}
}
}
}
columnArray->setNumElements(numberNonZero);
if (packed)
columnArray->setPacked();
y->setNumElements(0);
}
/* Return <code>x *A in <code>z</code> but
just for indices in y. */
void
ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * model,
const CoinIndexedVector * rowArray,
const CoinIndexedVector * y,
CoinIndexedVector * columnArray) const
{
columnArray->clear();
double * pi = rowArray->denseVector();
double * array = columnArray->denseVector();
int jColumn;
int numberToDo = y->getNumElements();
const int * which = y->getIndices();
assert (!rowArray->packedMode());
columnArray->setPacked();
for (jColumn=0;jColumn<numberToDo;jColumn++) {
int iColumn = which[jColumn];
double value = 0.0;
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
value += pi[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
value -= pi[iRow];
}
array[jColumn]=value;
}
}
/// returns number of elements in column part of basis,
CoinBigIndex
ClpPlusMinusOneMatrix::countBasis(ClpSimplex * model,
const int * whichColumn,
int numberBasic,
int & numberColumnBasic)
{
int i;
CoinBigIndex numberElements=0;
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
}
return numberElements;
}
void
ClpPlusMinusOneMatrix::fillBasis(ClpSimplex * model,
const int * whichColumn,
int & numberColumnBasic,
int * indexRowU, int * start,
int * rowCount, int * columnCount,
double * elementU)
{
int i;
CoinBigIndex numberElements=start[0];
assert (columnOrdered_);
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=1.0;
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=-1.0;
}
start[i+1]=numberElements;
columnCount[i]=numberElements-start[i];
}
}
/* Unpacks a column into an CoinIndexedvector
*/
void
ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
CoinIndexedVector * rowArray,
int iColumn) const
{
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
rowArray->add(iRow,1.0);
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
rowArray->add(iRow,-1.0);
}
}
/* Unpacks a column into an CoinIndexedvector
** in packed foramt
Note that model is NOT const. Bounds and objective could
be modified if doing column generation (just for this variable) */
void
ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
CoinIndexedVector * rowArray,
int iColumn) const
{
int * index = rowArray->getIndices();
double * array = rowArray->denseVector();
int number = 0;
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
array[number]=1.0;
index[number++]=iRow;
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
array[number]=-1.0;
index[number++]=iRow;
}
rowArray->setNumElements(number);
rowArray->setPackedMode(true);
}
/* Adds multiple of a column into an CoinIndexedvector
You can use quickAdd to add to vector */
void
ClpPlusMinusOneMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
int iColumn, double multiplier) const
{
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
rowArray->quickAdd(iRow,multiplier);
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
rowArray->quickAdd(iRow,-multiplier);
}
}
/* Adds multiple of a column into an array */
void
ClpPlusMinusOneMatrix::add(const ClpSimplex * model,double * array,
int iColumn, double multiplier) const
{
CoinBigIndex j=startPositive_[iColumn];
for (;j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
array[iRow] += multiplier;
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
array[iRow] -= multiplier;
}
}
// Return a complete CoinPackedMatrix
CoinPackedMatrix *
ClpPlusMinusOneMatrix::getPackedMatrix() const
{
if (!matrix_) {
int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
int numberElements = startPositive_[numberMajor];
double * elements = new double [numberElements];
CoinBigIndex j=0;
int i;
for (i=0;i<numberMajor;i++) {
for (;j<startNegative_[i];j++) {
elements[j]=1.0;
}
for (;j<startPositive_[i+1];j++) {
elements[j]=-1.0;
}
}
matrix_ = new CoinPackedMatrix(columnOrdered_,numberMinor,numberMajor,
getNumElements(),
elements,indices_,
startPositive_,getVectorLengths());
delete [] elements;
delete [] lengths_;
lengths_ = NULL;
}
return matrix_;
}
/* A vector containing the elements in the packed matrix. Note that there
might be gaps in this list, entries that do not belong to any
major-dimension vector. To get the actual elements one should look at
this vector together with vectorStarts and vectorLengths. */
const double *
ClpPlusMinusOneMatrix::getElements() const
{
if (!matrix_)
getPackedMatrix();
return matrix_->getElements();
}
const CoinBigIndex *
ClpPlusMinusOneMatrix::getVectorStarts() const
{
return startPositive_;
}
/* The lengths of the major-dimension vectors. */
const int *
ClpPlusMinusOneMatrix::getVectorLengths() const
{
if (!lengths_) {
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
lengths_ = new int [numberMajor];
int i;
for (i=0;i<numberMajor;i++) {
lengths_[i]=startPositive_[i+1]-startPositive_[i];
}
}
return lengths_;
}
/* Delete the columns whose indices are listed in <code>indDel</code>. */
void
ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
{
int iColumn;
CoinBigIndex newSize=startPositive_[numberColumns_];;
int numberBad=0;
// Use array to make sure we can have duplicates
int * which = new int[numberColumns_];
memset(which,0,numberColumns_*sizeof(int));
int nDuplicate=0;
for (iColumn=0;iColumn<numDel;iColumn++) {
int jColumn = indDel[iColumn];
if (jColumn<0||jColumn>=numberColumns_) {
numberBad++;
} else {
newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
if (which[jColumn])
nDuplicate++;
else
which[jColumn]=1;
}
}
if (numberBad)
throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
int newNumber = numberColumns_-numDel+nDuplicate;
// Get rid of temporary arrays
delete [] lengths_;
lengths_=NULL;
delete matrix_;
matrix_= NULL;
CoinBigIndex * newPositive = new CoinBigIndex [newNumber+1];
CoinBigIndex * newNegative = new CoinBigIndex [newNumber];
int * newIndices = new int [newSize];
newNumber=0;
newSize=0;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
if (!which[iColumn]) {
CoinBigIndex start,end;
CoinBigIndex i;
start = startPositive_[iColumn];
end=startNegative_[iColumn];
newPositive[newNumber]=newSize;
for (i=start;i<end;i++)
newIndices[newSize++]=indices_[i];
start = startNegative_[iColumn];
end=startPositive_[iColumn+1];
newNegative[newNumber++]=newSize;
for (i=start;i<end;i++)
newIndices[newSize++]=indices_[i];
}
}
newPositive[newNumber]=newSize;
delete [] which;
delete [] startPositive_;
startPositive_= newPositive;
delete [] startNegative_;
startNegative_= newNegative;
delete [] indices_;
indices_= newIndices;
numberColumns_ = newNumber;
}
/* Delete the rows whose indices are listed in <code>indDel</code>. */
void
ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
{
int iRow;
int numberBad=0;
// Use array to make sure we can have duplicates
int * which = new int[numberRows_];
memset(which,0,numberRows_*sizeof(int));
int nDuplicate=0;
for (iRow=0;iRow<numDel;iRow++) {
int jRow = indDel[iRow];
if (jRow<0||jRow>=numberRows_) {
numberBad++;
} else {
if (which[jRow])
nDuplicate++;
else
which[jRow]=1;
}
}
if (numberBad)
throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
CoinBigIndex iElement;
CoinBigIndex numberElements=startPositive_[numberColumns_];
CoinBigIndex newSize=0;
for (iElement=0;iElement<numberElements;iElement++) {
iRow = indices_[iElement];
if (!which[iRow])
newSize++;
}
int newNumber = numberRows_-numDel+nDuplicate;
// Get rid of temporary arrays
delete [] lengths_;
lengths_=NULL;
delete matrix_;
matrix_= NULL;
int * newIndices = new int [newSize];
newSize=0;
int iColumn;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
CoinBigIndex start,end;
CoinBigIndex i;
start = startPositive_[iColumn];
end=startNegative_[iColumn];
startPositive_[newNumber]=newSize;
for (i=start;i<end;i++) {
iRow = indices_[i];
if (!which[iRow])
newIndices[newSize++]=iRow;
}
start = startNegative_[iColumn];
end=startPositive_[iColumn+1];
startNegative_[newNumber]=newSize;
for (i=start;i<end;i++) {
iRow = indices_[i];
if (!which[iRow])
newIndices[newSize++]=iRow;
}
}
startPositive_[numberColumns_]=newSize;
delete [] which;
delete [] indices_;
indices_= newIndices;
numberRows_ = newNumber;
}
bool
ClpPlusMinusOneMatrix::isColOrdered() const
{
return columnOrdered_;
}
/* Number of entries in the packed matrix. */
CoinBigIndex
ClpPlusMinusOneMatrix::getNumElements() const
{
int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
if (startPositive_)
return startPositive_[numberMajor];
else
return 0;
}
// pass in copy (object takes ownership)
void
ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
bool columnOrdered, int * indices,
CoinBigIndex * startPositive, CoinBigIndex * startNegative)
{
columnOrdered_=columnOrdered;
startPositive_ = startPositive;
startNegative_ = startNegative;
indices_ = indices;
numberRows_=numberRows;
numberColumns_=numberColumns;
// Check valid
checkValid(false);
}
// Just checks matrix valid - will say if dimensions not quite right if detail
void
ClpPlusMinusOneMatrix::checkValid(bool detail) const
{
int maxIndex=-1;
int minIndex= columnOrdered_ ? numberRows_ : numberColumns_;
int number= !columnOrdered_ ? numberRows_ : numberColumns_;
int numberElements = getNumElements();
CoinBigIndex last=-1;
int bad=0;
for (int i=0;i<number;i++) {
if(startPositive_[i]<last)
bad++;
else
last = startPositive_[i];
if(startNegative_[i]<last)
bad++;
else
last = startNegative_[i];
}
if(startPositive_[number]<last)
bad++;
CoinAssertHint(!bad,"starts are not monotonic");
for (CoinBigIndex cbi=0;cbi<numberElements;cbi++) {
maxIndex = CoinMax(indices_[cbi],maxIndex);
minIndex = CoinMin(indices_[cbi],minIndex);
}
CoinAssert(maxIndex<(columnOrdered_ ? numberRows_ : numberColumns_));
CoinAssert(minIndex>=0);
if (detail) {
if (minIndex>0||maxIndex+1<(columnOrdered_ ? numberRows_ : numberColumns_))
printf("Not full range of indices - %d to %d\n",minIndex,maxIndex);
}
}
/* Given positive integer weights for each row fills in sum of weights
for each column (and slack).
Returns weights vector
*/
CoinBigIndex *
ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
{
int numberRows = model->numberRows();
int numberColumns =model->numberColumns();
int number = numberRows+numberColumns;
CoinBigIndex * weights = new CoinBigIndex[number];
int i;
for (i=0;i<numberColumns;i++) {
CoinBigIndex j;
CoinBigIndex count=0;
for (j=startPositive_[i];j<startPositive_[i+1];j++) {
int iRow=indices_[j];
count += inputWeights[iRow];
}
weights[i]=count;
}
for (i=0;i<numberRows;i++) {
weights[i+numberColumns]=inputWeights[i];
}
return weights;
}
// Append Columns
void
ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
{
int iColumn;
CoinBigIndex size=0;
int numberBad=0;
for (iColumn=0;iColumn<number;iColumn++) {
int n=columns[iColumn]->getNumElements();
const double * element = columns[iColumn]->getElements();
size += n;
int i;
for (i=0;i<n;i++) {
if (fabs(element[i])!=1.0)
numberBad++;
}
}
if (numberBad)
throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
// Get rid of temporary arrays
delete [] lengths_;
lengths_=NULL;
delete matrix_;
matrix_= NULL;
int numberNow = startPositive_[numberColumns_];
CoinBigIndex * temp;
temp = new CoinBigIndex [numberColumns_+1+number];
memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(CoinBigIndex ));
delete [] startPositive_;
startPositive_= temp;
temp = new CoinBigIndex [numberColumns_+number];
memcpy(temp,startNegative_,numberColumns_*sizeof(CoinBigIndex ));
delete [] startNegative_;
startNegative_= temp;
int * temp2 = new int [numberNow+size];
memcpy(temp2,indices_,numberNow*sizeof(int));
delete [] indices_;
indices_= temp2;
// now add
size=numberNow;
for (iColumn=0;iColumn<number;iColumn++) {
int n=columns[iColumn]->getNumElements();
const int * row = columns[iColumn]->getIndices();
const double * element = columns[iColumn]->getElements();
int i;
for (i=0;i<n;i++) {
if (element[i]==1.0)
indices_[size++]=row[i];
}
startNegative_[iColumn+numberColumns_]=size;
for (i=0;i<n;i++) {
if (element[i]==-1.0)
indices_[size++]=row[i];
}
startPositive_[iColumn+numberColumns_+1]=size;
}
numberColumns_ += number;
}
// Append Rows
void
ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
{
// Allocate arrays to use for counting
int * countPositive = new int [numberColumns_+1];
memset(countPositive,0,numberColumns_*sizeof(int));
int * countNegative = new int [numberColumns_];
memset(countNegative,0,numberColumns_*sizeof(int));
int iRow;
CoinBigIndex size=0;
int numberBad=0;
for (iRow=0;iRow<number;iRow++) {
int n=rows[iRow]->getNumElements();
const int * column = rows[iRow]->getIndices();
const double * element = rows[iRow]->getElements();
size += n;
int i;
for (i=0;i<n;i++) {
int iColumn = column[i];
if (element[i]==1.0)
countPositive[iColumn]++;
else if (element[i]==-1.0)
countNegative[iColumn]++;
else
numberBad++;
}
}
if (numberBad)
throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
// Get rid of temporary arrays
delete [] lengths_;
lengths_=NULL;
delete matrix_;
matrix_= NULL;
int numberNow = startPositive_[numberColumns_];
int * newIndices = new int [numberNow+size];
// Update starts and turn counts into positions
// also move current indices
int iColumn;
CoinBigIndex numberAdded=0;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
int n,move;
CoinBigIndex now;
now = startPositive_[iColumn];
move = startNegative_[iColumn]-now;
n = countPositive[iColumn];
startPositive_[iColumn] += numberAdded;
memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
countPositive[iColumn]= startNegative_[iColumn]+numberAdded;
numberAdded += n;
now = startNegative_[iColumn];
move = startPositive_[iColumn+1]-now;
n = countNegative[iColumn];
startNegative_[iColumn] += numberAdded;
memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
countNegative[iColumn]= startPositive_[iColumn+1]+numberAdded;
numberAdded += n;
}
delete [] indices_;
indices_ = newIndices;
startPositive_[numberColumns_] += numberAdded;
// Now put in
for (iRow=0;iRow<number;iRow++) {
int newRow = numberRows_+iRow;
int n=rows[iRow]->getNumElements();
const int * column = rows[iRow]->getIndices();
const double * element = rows[iRow]->getElements();
int i;
for (i=0;i<n;i++) {
int iColumn = column[i];
int put;
if (element[i]==1.0) {
put = countPositive[iColumn];
countPositive[iColumn] = put+1;
} else {
put = countNegative[iColumn];
countNegative[iColumn] = put+1;
}
indices_[put]=newRow;
}
}
delete [] countPositive;
delete [] countNegative;
numberRows_ += number;
}
/* Returns largest and smallest elements of both signs.
Largest refers to largest absolute value.
*/
void
ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
double & smallestPositive, double & largestPositive)
{
int iColumn;
bool plusOne=false;
bool minusOne=false;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
if (startNegative_[iColumn]>startPositive_[iColumn])
plusOne=true;
if (startPositive_[iColumn+1]>startNegative_[iColumn])
minusOne=true;
}
if (minusOne) {
smallestNegative=-1.0;
largestNegative=-1.0;
} else {
smallestNegative=0.0;
largestNegative=0.0;
}
if (plusOne) {
smallestPositive=1.0;
largestPositive=1.0;
} else {
smallestPositive=0.0;
largestPositive=0.0;
}
}
// Says whether it can do partial pricing
bool
ClpPlusMinusOneMatrix::canDoPartialPricing() const
{
return true;
}
// Partial pricing
void
ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
int & bestSequence, int & numberWanted)
{
numberWanted=currentWanted_;
int start = (int) (startFraction*numberColumns_);
int end = CoinMin((int) (endFraction*numberColumns_+1),numberColumns_);
CoinBigIndex j;
double tolerance=model->currentDualTolerance();
double * reducedCost = model->djRegion();
const double * duals = model->dualRowSolution();
const double * cost = model->costRegion();
double bestDj;
if (bestSequence>=0)
bestDj = fabs(reducedCost[bestSequence]);
else
bestDj=tolerance;
int sequenceOut = model->sequenceOut();
int saveSequence = bestSequence;
int iSequence;
for (iSequence=start;iSequence<end;iSequence++) {
if (iSequence!=sequenceOut) {
double value;
ClpSimplex::Status status = model->getStatus(iSequence);
switch(status) {
case ClpSimplex::basic:
case ClpSimplex::isFixed:
break;
case ClpSimplex::isFree:
case ClpSimplex::superBasic:
value=cost[iSequence];
j=startPositive_[iSequence];
for (;j<startNegative_[iSequence];j++) {
int iRow = indices_[j];
value -= duals[iRow];
}
for (;j<startPositive_[iSequence+1];j++) {
int iRow = indices_[j];
value += duals[iRow];
}
value = fabs(value);
if (value>FREE_ACCEPT*tolerance) {
numberWanted--;
// we are going to bias towards free (but only if reasonable)
value *= FREE_BIAS;
if (value>bestDj) {
// check flagged variable and correct dj
if (!model->flagged(iSequence)) {
bestDj=value;
bestSequence = iSequence;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
break;
case ClpSimplex::atUpperBound:
value=cost[iSequence];
j=startPositive_[iSequence];
for (;j<startNegative_[iSequence];j++) {
int iRow = indices_[j];
value -= duals[iRow];
}
for (;j<startPositive_[iSequence+1];j++) {
int iRow = indices_[j];
value += duals[iRow];
}
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!model->flagged(iSequence)) {
bestDj=value;
bestSequence = iSequence;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
break;
case ClpSimplex::atLowerBound:
value=cost[iSequence];
j=startPositive_[iSequence];
for (;j<startNegative_[iSequence];j++) {
int iRow = indices_[j];
value -= duals[iRow];
}
for (;j<startPositive_[iSequence+1];j++) {
int iRow = indices_[j];
value += duals[iRow];
}
value = -value;
if (value>tolerance) {
numberWanted--;
if (value>bestDj) {
// check flagged variable and correct dj
if (!model->flagged(iSequence)) {
bestDj=value;
bestSequence = iSequence;
} else {
// just to make sure we don't exit before got something
numberWanted++;
}
}
}
break;
}
}
if (!numberWanted)
break;
}
if (bestSequence!=saveSequence) {
// recompute dj
double value=cost[bestSequence];
j=startPositive_[bestSequence];
for (;j<startNegative_[bestSequence];j++) {
int iRow = indices_[j];
value -= duals[iRow];
}
for (;j<startPositive_[bestSequence+1];j++) {
int iRow = indices_[j];
value += duals[iRow];
}
reducedCost[bestSequence] = value;
savedBestSequence_ = bestSequence;
savedBestDj_ = reducedCost[savedBestSequence_];
}
currentWanted_=numberWanted;
}
// Allow any parts of a created CoinMatrix to be deleted
void
ClpPlusMinusOneMatrix::releasePackedMatrix() const
{
delete matrix_;
delete [] lengths_;
matrix_=NULL;
lengths_=NULL;
}
/* Returns true if can combine transposeTimes and subsetTransposeTimes
and if it would be faster */
bool
ClpPlusMinusOneMatrix::canCombine(const ClpSimplex * model,
const CoinIndexedVector * pi) const
{
int numberInRowArray = pi->getNumElements();
int numberRows = model->numberRows();
bool packed = pi->packedMode();
// factor should be smaller if doing both with two pi vectors
double factor = 0.27;
// We may not want to do by row if there may be cache problems
// It would be nice to find L2 cache size - for moment 512K
// Be slightly optimistic
if (numberColumns_*sizeof(double)>1000000) {
if (numberRows*10<numberColumns_)
factor *= 0.333333333;
else if (numberRows*4<numberColumns_)
factor *= 0.5;
else if (numberRows*2<numberColumns_)
factor *= 0.66666666667;
//if (model->numberIterations()%50==0)
//printf("%d nonzero\n",numberInRowArray);
}
// if not packed then bias a bit more towards by column
if (!packed)
factor *= 0.9;
return (numberInRowArray>factor*numberRows||!model->rowCopy());
}
// These have to match ClpPrimalColumnSteepest version
#define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
// Updates two arrays for steepest
void
ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model,
const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
CoinIndexedVector * spare,
double referenceIn, double devex,
// Array for exact devex to say what is in reference framework
unsigned int * reference,
double * weights, double scaleFactor)
{
// put row of tableau in dj1
double * pi = pi1->denseVector();
int numberNonZero=0;
int * index = dj1->getIndices();
double * array = dj1->denseVector();
int numberInRowArray = pi1->getNumElements();
double zeroTolerance = model->factorization()->zeroTolerance();
bool packed = pi1->packedMode();
// do by column
int iColumn;
assert (!spare->getNumElements());
double * piWeight = pi2->denseVector();
assert (!pi2->packedMode());
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=1.0;
// Note scale factor was -1.0
if (packed) {
// need to expand pi into y
assert(spare->capacity()>=model->numberRows());
double * piOld = pi;
pi = spare->denseVector();
const int * whichRow = pi1->getIndices();
int i;
// modify pi so can collapse to one loop
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=piOld[i];
}
CoinBigIndex j;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double value = 0.0;
for (j=startPositive_[iColumn];j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
value -= pi[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
value += pi[iRow];
}
if (fabs(value)>zeroTolerance) {
// and do other array
double modification = 0.0;
for (j=startPositive_[iColumn];j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
modification += piWeight[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
modification -= piWeight[iRow];
}
double thisWeight = weights[iColumn];
double pivot = value*scaleFactor;
double pivotSquared = pivot * pivot;
thisWeight += pivotSquared * devex + pivot * modification;
if (thisWeight<DEVEX_TRY_NORM) {
if (referenceIn<0.0) {
// steepest
thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
} else {
// exact
thisWeight = referenceIn*pivotSquared;
if (reference(iColumn))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[iColumn] = thisWeight;
if (!killDjs) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
}
// zero out
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=0.0;
}
} else {
CoinBigIndex j;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double value = 0.0;
for (j=startPositive_[iColumn];j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
value -= pi[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
value += pi[iRow];
}
if (fabs(value)>zeroTolerance) {
// and do other array
double modification = 0.0;
for (j=startPositive_[iColumn];j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
modification += piWeight[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
modification -= piWeight[iRow];
}
double thisWeight = weights[iColumn];
double pivot = value*scaleFactor;
double pivotSquared = pivot * pivot;
thisWeight += pivotSquared * devex + pivot * modification;
if (thisWeight<DEVEX_TRY_NORM) {
if (referenceIn<0.0) {
// steepest
thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
} else {
// exact
thisWeight = referenceIn*pivotSquared;
if (reference(iColumn))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[iColumn] = thisWeight;
if (!killDjs) {
array[iColumn]=value;
index[numberNonZero++]=iColumn;
}
}
}
}
dj1->setNumElements(numberNonZero);
spare->setNumElements(0);
if (packed)
dj1->setPackedMode(true);
}
// Updates second array for steepest and does devex weights
void
ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex * model,
CoinIndexedVector * dj1,
const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
double referenceIn, double devex,
// Array for exact devex to say what is in reference framework
unsigned int * reference,
double * weights, double scaleFactor)
{
int number = dj1->getNumElements();
const int * index = dj1->getIndices();
double * array = dj1->denseVector();
assert( dj1->packedMode());
double * piWeight = pi2->denseVector();
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=1.0;
for (int k=0;k<number;k++) {
int iColumn = index[k];
double pivot = array[k]*scaleFactor;
if (killDjs)
array[k]=0.0;
// and do other array
double modification = 0.0;
CoinBigIndex j;
for (j=startPositive_[iColumn];j<startNegative_[iColumn];j++) {
int iRow = indices_[j];
modification += piWeight[iRow];
}
for (;j<startPositive_[iColumn+1];j++) {
int iRow = indices_[j];
modification -= piWeight[iRow];
}
double thisWeight = weights[iColumn];
double pivotSquared = pivot * pivot;
thisWeight += pivotSquared * devex + pivot * modification;
if (thisWeight<DEVEX_TRY_NORM) {
if (referenceIn<0.0) {
// steepest
thisWeight = CoinMax(DEVEX_TRY_NORM,DEVEX_ADD_ONE+pivotSquared);
} else {
// exact
thisWeight = referenceIn*pivotSquared;
if (reference(iColumn))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[iColumn] = thisWeight;
}
}
/* Set the dimensions of the matrix. In effect, append new empty
columns/rows to the matrix. A negative number for either dimension
means that that dimension doesn't change. Otherwise the new dimensions
MUST be at least as large as the current ones otherwise an exception
is thrown. */
void
ClpPlusMinusOneMatrix::setDimensions(int newnumrows, int newnumcols)
{
if (newnumrows < 0)
newnumrows = numberRows_;
if (newnumrows < numberRows_)
throw CoinError("Bad new rownum (less than current)",
"setDimensions", "CoinPackedMatrix");
if (newnumcols < 0)
newnumcols = numberColumns_;
if (newnumcols < numberColumns_)
throw CoinError("Bad new colnum (less than current)",
"setDimensions", "CoinPackedMatrix");
int number = 0;
int length = 0;
if (columnOrdered_) {
length = numberColumns_;
numberColumns_ = newnumcols;
number = numberColumns_;
} else {
length = numberRows_;
numberRows_ = newnumrows;
number = numberRows_;
}
if (number > length) {
CoinBigIndex * temp;
int i;
CoinBigIndex end = startPositive_[length];
temp = new CoinBigIndex [number+1];
memcpy(temp,startPositive_,(length+1)*sizeof(CoinBigIndex));
delete [] startPositive_;
for (i=length+1;i<number+1;i++)
temp[i]= end;
startPositive_=temp;
temp = new CoinBigIndex [number];
memcpy(temp,startNegative_,length*sizeof(CoinBigIndex));
delete [] startNegative_;
for (i=length;i<number;i++)
temp[i]= end;
startNegative_=temp;
}
}
#ifndef SLIM_CLP
/* Append a set of rows/columns to the end of the matrix. Returns number of errors
i.e. if any of the new rows/columns contain an index that's larger than the
number of columns-1/rows-1 (if numberOther>0) or duplicates
If 0 then rows, 1 if columns */
int
ClpPlusMinusOneMatrix::appendMatrix(int number, int type,
const CoinBigIndex * starts, const int * index,
const double * element, int numberOther)
{
int numberErrors=0;
// make into CoinPackedVector
CoinPackedVectorBase ** vectors=
new CoinPackedVectorBase * [number];
int iVector;
for (iVector=0;iVector<number;iVector++) {
int iStart = starts[iVector];
vectors[iVector] =
new CoinPackedVector(starts[iVector+1]-iStart,
index+iStart,element+iStart);
}
if (type==0) {
// rows
appendRows(number,vectors);
} else {
// columns
appendCols(number,vectors);
}
for (iVector=0;iVector<number;iVector++)
delete vectors[iVector];
delete [] vectors;
return numberErrors;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1