// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include <cstdio>
#include "CoinPragma.hpp"
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
//#define THREAD
#include "ClpSimplex.hpp"
#include "ClpSimplexDual.hpp"
#include "ClpFactorization.hpp"
#ifndef SLIM_CLP
#include "ClpQuadraticObjective.hpp"
#endif
// at end to get min/max!
#include "ClpPackedMatrix.hpp"
#include "ClpMessage.hpp"
//#define COIN_PREFETCH
#ifdef COIN_PREFETCH
#if 1
#define coin_prefetch(mem) \
__asm__ __volatile__ ("prefetchnta %0" : : "m" (*((char *)(mem))))
#else
#define coin_prefetch(mem) \
__asm__ __volatile__ ("prefetch %0" : : "m" (*((char *)(mem))))
#endif
#else
// dummy
#define coin_prefetch(mem)
#endif
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpPackedMatrix::ClpPackedMatrix ()
: ClpMatrixBase(),
matrix_(NULL),
numberActiveColumns_(0),
flags_(2),
rowCopy_(NULL),
columnCopy_(NULL)
{
setType(1);
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs)
: ClpMatrixBase(rhs)
{
matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
numberActiveColumns_ = rhs.numberActiveColumns_;
flags_ = rhs.flags_;
int numberRows = getNumRows();
if (rhs.rhsOffset_&&numberRows) {
rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
} else {
rhsOffset_=NULL;
}
if (rhs.rowCopy_) {
assert ((flags_&4)!=0);
rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
} else {
rowCopy_ = NULL;
}
if (rhs.columnCopy_) {
assert ((flags_&8)!=0);
columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
} else {
columnCopy_ = NULL;
}
}
//-------------------------------------------------------------------
// assign matrix (for space reasons)
//-------------------------------------------------------------------
ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs)
: ClpMatrixBase()
{
matrix_ = rhs;
flags_ = 2;
numberActiveColumns_ = matrix_->getNumCols();
rowCopy_ = NULL;
columnCopy_=NULL;
setType(1);
}
ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs)
: ClpMatrixBase()
{
matrix_ = new CoinPackedMatrix(rhs);
numberActiveColumns_ = matrix_->getNumCols();
rowCopy_ = NULL;
flags_ = 2;
columnCopy_=NULL;
setType(1);
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpPackedMatrix::~ClpPackedMatrix ()
{
delete matrix_;
delete rowCopy_;
delete columnCopy_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpPackedMatrix &
ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
{
if (this != &rhs) {
ClpMatrixBase::operator=(rhs);
delete matrix_;
matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
numberActiveColumns_ = rhs.numberActiveColumns_;
flags_ = rhs.flags_;
delete rowCopy_;
delete columnCopy_;
if (rhs.rowCopy_) {
assert ((flags_&4)!=0);
rowCopy_ = new ClpPackedMatrix2(*rhs.rowCopy_);
} else {
rowCopy_ = NULL;
}
if (rhs.columnCopy_) {
assert ((flags_&8)!=0);
columnCopy_ = new ClpPackedMatrix3(*rhs.columnCopy_);
} else {
columnCopy_ = NULL;
}
}
return *this;
}
//-------------------------------------------------------------------
// Clone
//-------------------------------------------------------------------
ClpMatrixBase * ClpPackedMatrix::clone() const
{
return new ClpPackedMatrix(*this);
}
/* Subset clone (without gaps). Duplicates are allowed
and order is as given */
ClpMatrixBase *
ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
int numberColumns,
const int * whichColumns) const
{
return new ClpPackedMatrix(*this, numberRows, whichRows,
numberColumns, whichColumns);
}
/* Subset constructor (without gaps). Duplicates are allowed
and order is as given */
ClpPackedMatrix::ClpPackedMatrix (
const ClpPackedMatrix & rhs,
int numberRows, const int * whichRows,
int numberColumns, const int * whichColumns)
: ClpMatrixBase(rhs)
{
matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
numberColumns,whichColumns);
numberActiveColumns_ = matrix_->getNumCols();
rowCopy_ = NULL;
flags_ = rhs.flags_;
columnCopy_=NULL;
}
ClpPackedMatrix::ClpPackedMatrix (
const CoinPackedMatrix & rhs,
int numberRows, const int * whichRows,
int numberColumns, const int * whichColumns)
: ClpMatrixBase()
{
matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
numberColumns,whichColumns);
numberActiveColumns_ = matrix_->getNumCols();
rowCopy_ = NULL;
flags_ = 2;
columnCopy_=NULL;
setType(1);
}
/* Returns a new matrix in reverse order without gaps */
ClpMatrixBase *
ClpPackedMatrix::reverseOrderedCopy() const
{
ClpPackedMatrix * copy = new ClpPackedMatrix();
copy->matrix_= new CoinPackedMatrix();
copy->matrix_->setExtraGap(0.0);
copy->matrix_->setExtraMajor(0.0);
copy->matrix_->reverseOrderedCopyOf(*matrix_);
//copy->matrix_->removeGaps();
copy->numberActiveColumns_ = copy->matrix_->getNumCols();
copy->flags_ = flags_&(~2); // no gaps
return copy;
}
//unscaled versions
void
ClpPackedMatrix::times(double scalar,
const double * x, double * y) const
{
int iRow,iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
//memset(y,0,matrix_->getNumRows()*sizeof(double));
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
double value = scalar*x[iColumn];
if (value) {
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
y[iRow] += value*elementByColumn[j];
}
}
}
}
void
ClpPackedMatrix::transposeTimes(double scalar,
const double * x, double * y) const
{
int iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
if (!(flags_&2)) {
if (scalar==1.0) {
CoinBigIndex start=columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
CoinBigIndex next=columnStart[iColumn+1];
double value=y[iColumn];
for (j=start;j<next;j++) {
int jRow=row[j];
value += x[jRow]*elementByColumn[j];
}
start=next;
y[iColumn] = value;
}
} else if (scalar==-1.0) {
CoinBigIndex start=columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
CoinBigIndex next=columnStart[iColumn+1];
double value=y[iColumn];
for (j=start;j<next;j++) {
int jRow=row[j];
value -= x[jRow]*elementByColumn[j];
}
start=next;
y[iColumn] = value;
}
} else {
CoinBigIndex start=columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
CoinBigIndex next=columnStart[iColumn+1];
double value=0.0;
for (j=start;j<next;j++) {
int jRow=row[j];
value += x[jRow]*elementByColumn[j];
}
start=next;
y[iColumn] += value*scalar;
}
}
} else {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
double value=0.0;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int jRow=row[j];
value += x[jRow]*elementByColumn[j];
}
y[iColumn] += value*scalar;
}
}
}
void
ClpPackedMatrix::times(double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale) const
{
if (rowScale) {
int iRow,iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
double value = x[iColumn];
if (value) {
// scaled
value *= scalar*columnScale[iColumn];
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
y[iRow] += value*elementByColumn[j]*rowScale[iRow];
}
}
}
} else {
times(scalar,x,y);
}
}
void
ClpPackedMatrix::transposeTimes( double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale,
double * spare) const
{
if (rowScale) {
int iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
if (!spare) {
if (!(flags_&2)) {
CoinBigIndex start=columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
CoinBigIndex next=columnStart[iColumn+1];
double value=0.0;
// scaled
for (j=start;j<next;j++) {
int jRow=row[j];
value += x[jRow]*elementByColumn[j]*rowScale[jRow];
}
start=next;
y[iColumn] += value*scalar*columnScale[iColumn];
}
} else {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
double value=0.0;
// scaled
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int jRow=row[j];
value += x[jRow]*elementByColumn[j]*rowScale[jRow];
}
y[iColumn] += value*scalar*columnScale[iColumn];
}
}
} else {
// can use spare region
int iRow;
int numberRows = getNumRows();
for (iRow=0;iRow<numberRows;iRow++)
spare[iRow] = x[iRow]*rowScale[iRow];
if (!(flags_&2)) {
CoinBigIndex start=columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
CoinBigIndex next=columnStart[iColumn+1];
double value=0.0;
// scaled
for (j=start;j<next;j++) {
int jRow=row[j];
value += spare[jRow]*elementByColumn[j];
}
start=next;
y[iColumn] += value*scalar*columnScale[iColumn];
}
} else {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
double value=0.0;
// scaled
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int jRow=row[j];
value += spare[jRow]*elementByColumn[j];
}
y[iColumn] += value*scalar*columnScale[iColumn];
}
}
// no need to zero out
//for (iRow=0;iRow<numberRows;iRow++)
//spare[iRow] = 0.0;
}
} else {
transposeTimes(scalar,x,y);
}
}
/* Return <code>x * A + y</code> in <code>z</code>.
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix::transposeTimes(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();
int numberRows = model->numberRows();
#ifndef NO_RTTI
ClpPackedMatrix* rowCopy =
dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
#else
ClpPackedMatrix* rowCopy =
static_cast< ClpPackedMatrix*>(model->rowCopy());
#endif
bool packed = rowArray->packedMode();
double factor = 0.30;
// 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 (numberActiveColumns_*sizeof(double)>1000000) {
if (numberRows*10<numberActiveColumns_)
factor *= 0.333333333;
else if (numberRows*4<numberActiveColumns_)
factor *= 0.5;
else if (numberRows*2<numberActiveColumns_)
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;
assert (!y->getNumElements());
double multiplierX=0.8;
double factor2 = factor*multiplierX;
if (packed&&rowCopy_&&numberInRowArray>2&&numberInRowArray>factor2*numberRows&&
numberInRowArray<0.9*numberRows&&0) {
rowCopy_->transposeTimes(model,rowCopy->getPackedMatrix(),rowArray,y,columnArray);
return;
}
if (numberInRowArray>factor*numberRows||!rowCopy) {
// do by column
// If no gaps - can do a bit faster
if (!(flags_&2)||columnCopy_) {
transposeTimesByColumn( model, scalar,
rowArray, y, columnArray);
return;
}
int iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
const double * rowScale = model->rowScale();
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;
if (!rowScale) {
// modify pi so can collapse to one loop
if (scalar==-1.0) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=-piOld[i];
}
} else {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=scalar*piOld[i];
}
}
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
} else {
// scaled
// modify pi so can collapse to one loop
if (scalar==-1.0) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=-piOld[i]*rowScale[iRow];
}
} else {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=scalar*piOld[i]*rowScale[iRow];
}
}
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
const double * columnScale = model->columnScale();
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
value *= columnScale[iColumn];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
}
// zero out
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=0.0;
}
} else {
if (!rowScale) {
if (scalar==-1.0) {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
if (fabs(value)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn]=-value;
}
}
} else {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
value *= scalar;
if (fabs(value)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn]=value;
}
}
}
} else {
// scaled
if (scalar==-1.0) {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
const double * columnScale = model->columnScale();
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
value *= columnScale[iColumn];
if (fabs(value)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn]=-value;
}
}
} else {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
double value = 0.0;
CoinBigIndex j;
const double * columnScale = model->columnScale();
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
value *= scalar*columnScale[iColumn];
if (fabs(value)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn]=value;
}
}
}
}
}
columnArray->setNumElements(numberNonZero);
y->setNumElements(0);
} else {
// do by row
rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
}
if (packed)
columnArray->setPackedMode(true);
if (0) {
columnArray->checkClean();
int numberNonZero=columnArray->getNumElements();;
int * index = columnArray->getIndices();
double * array = columnArray->denseVector();
int i;
for (i=0;i<numberNonZero;i++) {
int j=index[i];
double value;
if (packed)
value=array[i];
else
value=array[j];
printf("Ti %d %d %g\n",i,j,value);
}
}
}
/* Return <code>x * scalar * A + y</code> in <code>z</code>.
Note - If x packed mode - then z packed mode
This does by column and knows no gaps
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
const CoinIndexedVector * rowArray,
CoinIndexedVector * y,
CoinIndexedVector * columnArray) const
{
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();
bool packed = rowArray->packedMode();
// do by column
int iColumn;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const double * elementByColumn = matrix_->getElements();
const double * rowScale = model->rowScale();
assert (!y->getNumElements());
assert (numberActiveColumns_>0);
if (packed) {
// need to expand pi into y
assert(y->capacity()>=model->numberRows());
double * piOld = pi;
pi = y->denseVector();
const int * whichRow = rowArray->getIndices();
int i;
if (!rowScale) {
// modify pi so can collapse to one loop
if (scalar==-1.0) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=-piOld[i];
}
} else {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=scalar*piOld[i];
}
}
#if 0
double value = 0.0;
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
#else
if (!columnCopy_)
gutsOfTransposeTimesUnscaled(pi,columnArray,zeroTolerance);
else
columnCopy_->transposeTimes(model,pi,columnArray);
numberNonZero = columnArray->getNumElements();
#endif
} else {
// scaled
// modify pi so can collapse to one loop
if (scalar==-1.0) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=-piOld[i]*rowScale[iRow];
}
} else {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=scalar*piOld[i]*rowScale[iRow];
}
}
const double * columnScale = model->columnScale();
#if 0
double value = 0.0;
double scale=columnScale[0];
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
value *= scale;
CoinBigIndex start = end;
scale = columnScale[iColumn+1];
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
if (model->getColumnStatus(iColumn+1)!=ClpSimplex::basic) {
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
}
value *= scale;
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
#else
if (!columnCopy_)
gutsOfTransposeTimesScaled(pi,columnScale,columnArray,zeroTolerance);
else
columnCopy_->transposeTimes(model,pi,columnArray);
numberNonZero = columnArray->getNumElements();
#endif
}
// zero out
int numberRows = model->numberRows();
if (numberInRowArray*4<numberRows) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=0.0;
}
} else {
CoinZeroN(pi,numberRows);
}
} else {
if (!rowScale) {
if (scalar==-1.0) {
double value = 0.0;
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[iColumn]=-value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
if (fabs(value)>zeroTolerance) {
array[iColumn]=-value;
index[numberNonZero++]=iColumn;
}
} else {
double value = 0.0;
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
value *= scalar;
CoinBigIndex start = end;
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[iColumn]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
value *= scalar;
if (fabs(value)>zeroTolerance) {
array[iColumn]=value;
index[numberNonZero++]=iColumn;
}
}
} else {
// scaled
if (scalar==-1.0) {
const double * columnScale = model->columnScale();
double value = 0.0;
double scale=columnScale[0];
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
value *= scale;
CoinBigIndex start = end;
end = columnStart[iColumn+2];
scale = columnScale[iColumn+1];
if (fabs(value)>zeroTolerance) {
array[iColumn]=-value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
}
value *= scale;
if (fabs(value)>zeroTolerance) {
array[iColumn]=-value;
index[numberNonZero++]=iColumn;
}
} else {
const double * columnScale = model->columnScale();
double value = 0.0;
double scale=columnScale[0]*scalar;
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
value *= scale;
CoinBigIndex start = end;
end = columnStart[iColumn+2];
scale = columnScale[iColumn+1]*scalar;
if (fabs(value)>zeroTolerance) {
array[iColumn]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
}
value *= scale;
if (fabs(value)>zeroTolerance) {
array[iColumn]=value;
index[numberNonZero++]=iColumn;
}
}
}
}
columnArray->setNumElements(numberNonZero);
y->setNumElements(0);
if (packed)
columnArray->setPackedMode(true);
}
/* Return <code>x * A + y</code> in <code>z</code>.
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix::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 = getIndices();
const CoinBigIndex * rowStart = getVectorStarts();
const double * element = getElements();
const int * whichRow = rowArray->getIndices();
bool packed = rowArray->packedMode();
if (numberInRowArray>2) {
// do by rows
// ** Row copy is already scaled
int iRow;
int i;
int numberOriginal = 0;
if (packed) {
#if 0
numberNonZero=0;
// and set up mark as char array
char * marked = (char *) (index+columnArray->capacity());
double * array2 = y->denseVector();
#ifdef CLP_DEBUG
for (i=0;i<numberActiveColumns_;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=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array2[iColumn] += value*element[j];
}
}
// 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
gutsOfTransposeTimesByRowGE3(rowArray,columnArray,y,zeroTolerance,scalar);
numberNonZero = columnArray->getNumElements();
#endif
} else {
double * markVector = y->denseVector();
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=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
index[numberNonZero++]=iColumn;
}
array[iColumn] += value*element[j];
}
}
// 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
int numberOriginal;
int i;
CoinBigIndex j;
numberNonZero=0;
double value;
if (packed) {
#if 0
int iRow0 = whichRow[0];
int iRow1 = whichRow[1];
double pi0 = pi[0];
double pi1 = pi[1];
if (rowStart[iRow0+1]-rowStart[iRow0]>
rowStart[iRow1+1]-rowStart[iRow1]) {
// do one with fewer first
iRow0=iRow1;
iRow1=whichRow[0];
pi0=pi1;
pi1=pi[0];
}
// and set up mark as char array
char * marked = (char *) (index+columnArray->capacity());
int * lookup = y->getIndices();
value = pi0*scalar;
for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
array[numberNonZero] = value2;
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
index[numberNonZero++]=iColumn;
}
numberOriginal = numberNonZero;
value = pi1*scalar;
for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
// I am assuming no zeros in matrix
if (marked[iColumn]) {
int iLookup = lookup[iColumn];
array[iLookup] += value2;
} else {
if (fabs(value2)>zeroTolerance) {
array[numberNonZero] = value2;
index[numberNonZero++]=iColumn;
}
}
}
// get rid of tiny values and zero out marked
int nDelete=0;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
marked[iColumn]=0;
if (fabs(array[i])<=zeroTolerance)
nDelete++;
}
if (nDelete) {
numberOriginal=numberNonZero;
numberNonZero=0;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
double value = array[i];
array[i]=0.0;
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
}
}
#else
gutsOfTransposeTimesByRowEQ2(rowArray,columnArray,y,zeroTolerance,scalar);
numberNonZero = columnArray->getNumElements();
#endif
} else {
int iRow = whichRow[0];
value = pi[iRow]*scalar;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
index[numberNonZero++]=iColumn;
array[iColumn] = value2;
}
iRow = whichRow[1];
value = pi[iRow]*scalar;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
// I am assuming no zeros in matrix
if (array[iColumn])
value2 += array[iColumn];
else
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;
CoinBigIndex j;
if (packed) {
#if 0
double value = pi[0]*scalar;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
if (fabs(value2)>zeroTolerance) {
array[numberNonZero] = value2;
index[numberNonZero++]=iColumn;
}
}
#else
gutsOfTransposeTimesByRowEQ1(rowArray,columnArray,zeroTolerance,scalar);
numberNonZero = columnArray->getNumElements();
#endif
} else {
double value = pi[iRow]*scalar;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double value2 = value*element[j];
if (fabs(value2)>zeroTolerance) {
index[numberNonZero++]=iColumn;
array[iColumn] = value2;
}
}
}
}
columnArray->setNumElements(numberNonZero);
y->setNumElements(0);
}
// Meat of transposeTimes by column when not scaled
void
ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * pi,CoinIndexedVector * output, const double zeroTolerance) const
{
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const double * elementByColumn = matrix_->getElements();
double value = 0.0;
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
int iColumn;
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
output->setNumElements(numberNonZero);
}
// Meat of transposeTimes by column when scaled
void
ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * pi,const double * columnScale,
CoinIndexedVector * output, const double zeroTolerance) const
{
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const double * elementByColumn = matrix_->getElements();
double value = 0.0;
double scale=columnScale[0];
CoinBigIndex j;
CoinBigIndex end = columnStart[1];
for (j=columnStart[0]; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
int iColumn;
for (iColumn=0;iColumn<numberActiveColumns_-1;iColumn++) {
value *= scale;
CoinBigIndex start = end;
scale = columnScale[iColumn+1];
end = columnStart[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
value *= scale;
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=iColumn;
}
output->setNumElements(numberNonZero);
}
// Meat of transposeTimes by row n > 2 if packed
void
ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * piVector, CoinIndexedVector * output,
CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
{
double * pi = piVector->denseVector();
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
int numberInRowArray = piVector->getNumElements();
const int * column = getIndices();
const CoinBigIndex * rowStart = getVectorStarts();
const double * element = getElements();
const int * whichRow = piVector->getIndices();
// ** Row copy is already scaled
int iRow;
int i;
// and set up mark as char array
char * marked = (char *) (index+output->capacity());
int * lookup = spareVector->getIndices();
#if 1
for (i=0;i<numberInRowArray;i++) {
iRow = whichRow[i];
double value = pi[i]*scalar;
CoinBigIndex j;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double elValue = element[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
array[numberNonZero] = value*elValue;
index[numberNonZero++]=iColumn;
} else {
int k = lookup[iColumn];
array[k] += value*elValue;
}
}
}
#else
int nextRow = whichRow[0];
CoinBigIndex nextStart = rowStart[nextRow];
CoinBigIndex nextEnd = rowStart[nextRow+1];
whichRow[numberInRowArray]=0; // for electricfence etc
for (i=0;i<numberInRowArray;i++) {
iRow = nextRow;
nextRow = whichRow[i+1];
CoinBigIndex start = nextStart;
CoinBigIndex end = nextEnd;
double value = pi[i]*scalar;
CoinBigIndex j;
nextRow = whichRow[i+1];
nextStart = rowStart[nextRow];
nextEnd = rowStart[nextRow+1];
for (j=start;j<end;j++) {
int iColumn = column[j];
double elValue = element[j];
if (!marked[iColumn]) {
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
array[numberNonZero] = value*elValue;
index[numberNonZero++]=iColumn;
} else {
int k = lookup[iColumn];
array[k] += value*elValue;
}
}
}
#endif
// get rid of tiny values and zero out marked
for (i=0;i<numberNonZero;i++) {
int iColumn = index[i];
marked[iColumn]=0;
double value = array[i];
while (fabs(value)<=tolerance) {
numberNonZero--;
value = array[numberNonZero];
iColumn = index[numberNonZero];
marked[iColumn]=0;
if (i<numberNonZero) {
array[numberNonZero]=0.0;
array[i] = value;
index[i] = iColumn;
} else {
array[i]=0.0;
value =1.0; // to force end of while
}
}
}
output->setNumElements(numberNonZero);
spareVector->setNumElements(0);
}
// Meat of transposeTimes by row n == 2 if packed
void
ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
{
double * pi = piVector->denseVector();
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
const int * column = getIndices();
const CoinBigIndex * rowStart = getVectorStarts();
const double * element = getElements();
const int * whichRow = piVector->getIndices();
int iRow0 = whichRow[0];
int iRow1 = whichRow[1];
double pi0 = pi[0];
double pi1 = pi[1];
if (rowStart[iRow0+1]-rowStart[iRow0]>
rowStart[iRow1+1]-rowStart[iRow1]) {
// do one with fewer first
iRow0=iRow1;
iRow1=whichRow[0];
pi0=pi1;
pi1=pi[0];
}
// and set up mark as char array
char * marked = (char *) (index+output->capacity());
int * lookup = spareVector->getIndices();
double value = pi0*scalar;
CoinBigIndex j;
for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
int iColumn = column[j];
double elValue = element[j];
double value2 = value*elValue;
array[numberNonZero] = value2;
marked[iColumn]=1;
lookup[iColumn]=numberNonZero;
index[numberNonZero++]=iColumn;
}
int numberOriginal = numberNonZero;
value = pi1*scalar;
for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
int iColumn = column[j];
double elValue = element[j];
double value2 = value*elValue;
// I am assuming no zeros in matrix
if (marked[iColumn]) {
int iLookup = lookup[iColumn];
array[iLookup] += value2;
} else {
if (fabs(value2)>tolerance) {
array[numberNonZero] = value2;
index[numberNonZero++]=iColumn;
}
}
}
// get rid of tiny values and zero out marked
int i;
int iFirst=numberNonZero;
for (i=0;i<numberOriginal;i++) {
int iColumn = index[i];
marked[iColumn]=0;
if (fabs(array[i])<=tolerance) {
if (numberNonZero>numberOriginal) {
numberNonZero--;
double value = array[numberNonZero];
array[numberNonZero]=0.0;
array[i]=value;
index[i]=index[numberNonZero];
} else {
iFirst = i;
}
}
}
if (iFirst<numberNonZero) {
int n=iFirst;
for (i=n;i<numberOriginal;i++) {
int iColumn = index[i];
double value = array[i];
array[i]=0.0;
if (fabs(value)>tolerance) {
array[n]=value;
index[n++]=iColumn;
}
}
for (;i<numberNonZero;i++) {
int iColumn = index[i];
double value = array[i];
array[i]=0.0;
array[n]=value;
index[n++]=iColumn;
}
numberNonZero=n;
}
output->setNumElements(numberNonZero);
spareVector->setNumElements(0);
}
// Meat of transposeTimes by row n == 1 if packed
void
ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
const double tolerance, const double scalar) const
{
double * pi = piVector->denseVector();
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
const int * column = getIndices();
const CoinBigIndex * rowStart = getVectorStarts();
const double * element = getElements();
int iRow=piVector->getIndices()[0];
numberNonZero=0;
CoinBigIndex j;
double value = pi[0]*scalar;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
double elValue = element[j];
double value2 = value*elValue;
if (fabs(value2)>tolerance) {
array[numberNonZero] = value2;
index[numberNonZero++]=iColumn;
}
}
output->setNumElements(numberNonZero);
}
/* Return <code>x *A in <code>z</code> but
just for indices in y.
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix::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;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
const double * rowScale = model->rowScale();
int numberToDo = y->getNumElements();
const int * which = y->getIndices();
assert (!rowArray->packedMode());
columnArray->setPacked();
if (!(flags_&2)&&numberToDo>5) {
// no gaps and a reasonable number
if (!rowScale) {
int iColumn = which[0];
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn+1];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
int iColumn = which[jColumn+1];
CoinBigIndex start = columnStart[iColumn];
CoinBigIndex end = columnStart[iColumn+1];
array[jColumn]=value;
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
}
array[jColumn]=value;
} else {
// scaled
const double * columnScale = model->columnScale();
int iColumn = which[0];
double value = 0.0;
double scale=columnScale[iColumn];
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn+1];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
for (jColumn=0;jColumn<numberToDo-1;jColumn++) {
int iColumn = which[jColumn+1];
value *= scale;
scale = columnScale[iColumn];
CoinBigIndex start = columnStart[iColumn];
CoinBigIndex end = columnStart[iColumn+1];
array[jColumn]=value;
value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
}
value *= scale;
array[jColumn]=value;
}
} else {
// gaps
if (!rowScale) {
for (jColumn=0;jColumn<numberToDo;jColumn++) {
int iColumn = which[jColumn];
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j];
}
array[jColumn]=value;
}
} else {
// scaled
const double * columnScale = model->columnScale();
for (jColumn=0;jColumn<numberToDo;jColumn++) {
int iColumn = which[jColumn];
double value = 0.0;
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
value *= columnScale[iColumn];
array[jColumn]=value;
}
}
}
}
/* Returns true if can combine transposeTimes and subsetTransposeTimes
and if it would be faster */
bool
ClpPackedMatrix::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.30;
// 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 (numberActiveColumns_*sizeof(double)>1000000) {
if (numberRows*10<numberActiveColumns_)
factor *= 0.333333333;
else if (numberRows*4<numberActiveColumns_)
factor *= 0.5;
else if (numberRows*2<numberActiveColumns_)
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())&&!(flags_&2));
}
#ifndef CLP_ALL_ONE_FILE
// These have to match ClpPrimalColumnSteepest version
#define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
#endif
// Updates two arrays for steepest
void
ClpPackedMatrix::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;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const double * elementByColumn = matrix_->getElements();
const double * rowScale = model->rowScale();
assert (!spare->getNumElements());
assert (numberActiveColumns_>0);
double * piWeight = pi2->denseVector();
assert (!pi2->packedMode());
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=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;
if (!rowScale) {
// modify pi so can collapse to one loop
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=piOld[i];
}
if (!columnCopy_) {
CoinBigIndex j;
CoinBigIndex end = columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+1];
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value -= pi[iRow]*elementByColumn[j];
}
if (fabs(value)>zeroTolerance) {
// and do other array
double modification = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j];
}
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;
}
}
}
} else {
// use special column copy
// reset back
if (killDjs)
scaleFactor=0.0;
columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
reference,weights,scaleFactor);
numberNonZero = dj1->getNumElements();
}
} else {
// scaled
// modify pi so can collapse to one loop
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=piOld[i]*rowScale[iRow];
}
// can also scale piWeight as not used again
int numberWeight = pi2->getNumElements();
const int * indexWeight = pi2->getIndices();
for (i=0;i<numberWeight;i++) {
int iRow = indexWeight[i];
piWeight[iRow] *= rowScale[iRow];
}
if (!columnCopy_) {
const double * columnScale = model->columnScale();
CoinBigIndex j;
CoinBigIndex end = columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+1];
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double scale=columnScale[iColumn];
double value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value -= pi[iRow]*elementByColumn[j];
}
value *= scale;
if (fabs(value)>zeroTolerance) {
double modification = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j];
}
modification *= scale;
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;
}
}
}
} else {
// use special column copy
// reset back
if (killDjs)
scaleFactor=0.0;
columnCopy_->transposeTimes2(model,pi,dj1,piWeight,referenceIn,devex,
reference,weights,scaleFactor);
numberNonZero = dj1->getNumElements();
}
}
// zero out
int numberRows = model->numberRows();
if (numberInRowArray*4<numberRows) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
pi[iRow]=0.0;
}
} else {
CoinZeroN(pi,numberRows);
}
} else {
if (!rowScale) {
CoinBigIndex j;
CoinBigIndex end = columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+1];
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value -= pi[iRow]*elementByColumn[j];
}
if (fabs(value)>zeroTolerance) {
// and do other array
double modification = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j];
}
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;
}
}
}
} else {
// scaled
// can also scale piWeight as not used again
int numberWeight = pi2->getNumElements();
const int * indexWeight = pi2->getIndices();
for (int i=0;i<numberWeight;i++) {
int iRow = indexWeight[i];
piWeight[iRow] *= rowScale[iRow];
}
const double * columnScale = model->columnScale();
CoinBigIndex j;
CoinBigIndex end = columnStart[0];
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex start = end;
end = columnStart[iColumn+1];
ClpSimplex::Status status = model->getStatus(iColumn);
if (status==ClpSimplex::basic||status==ClpSimplex::isFixed) continue;
double scale=columnScale[iColumn];
double value = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
value -= pi[iRow]*elementByColumn[j]*rowScale[iRow];
}
value *= scale;
if (fabs(value)>zeroTolerance) {
double modification = 0.0;
for (j=start; j<end;j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j];
}
modification *= scale;
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
ClpPackedMatrix::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());
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
const double * rowScale = model->rowScale();
double * piWeight = pi2->denseVector();
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=1.0;
if (!rowScale) {
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;
for (CoinBigIndex j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j];
}
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;
}
} else {
// scaled
const double * columnScale = model->columnScale();
for (int k=0;k<number;k++) {
int iColumn = index[k];
double pivot = array[k]*scaleFactor;
double scale=columnScale[iColumn];
if (killDjs)
array[k]=0.0;
// and do other array
double modification = 0.0;
for (CoinBigIndex j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
modification += piWeight[iRow]*elementByColumn[j]*rowScale[iRow];
}
double thisWeight = weights[iColumn];
modification *= scale;
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;
}
}
}
/// returns number of elements in column part of basis,
CoinBigIndex
ClpPackedMatrix::countBasis(ClpSimplex * model,
const int * whichColumn,
int numberBasic,
int & numberColumnBasic)
{
const int * columnLength = matrix_->getVectorLengths();
int i;
CoinBigIndex numberElements=0;
// just count - can be over so ignore zero problem
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
numberElements += columnLength[iColumn];
}
return numberElements;
}
void
ClpPackedMatrix::fillBasis(ClpSimplex * model,
const int * whichColumn,
int & numberColumnBasic,
int * indexRowU, int * start,
int * rowCount, int * columnCount,
double * elementU)
{
const int * columnLength = matrix_->getVectorLengths();
int i;
CoinBigIndex numberElements=start[0];
// fill
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const double * rowScale = model->rowScale();
const int * row = matrix_->getIndices();
const double * elementByColumn = matrix_->getElements();
if ((flags_&1)==0) {
if (!rowScale) {
// no scaling
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow=row[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=elementByColumn[j];
}
start[i+1]=numberElements;
columnCount[i]=columnLength[iColumn];
}
} else {
// scaling
const double * columnScale = model->columnScale();
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
CoinBigIndex j;
double scale = columnScale[iColumn];
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=
elementByColumn[j]*scale*rowScale[iRow];
}
start[i+1]=numberElements;
columnCount[i]=columnLength[iColumn];
}
}
} else {
// there are zero elements so need to look more closely
if (!rowScale) {
// no scaling
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
CoinBigIndex j;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
double value = elementByColumn[j];
if (value) {
int iRow = row[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=value;
}
}
start[i+1]=numberElements;
columnCount[i]=numberElements-start[i];
}
} else {
// scaling
const double * columnScale = model->columnScale();
for (i=0;i<numberColumnBasic;i++) {
int iColumn = whichColumn[i];
CoinBigIndex j;
double scale = columnScale[iColumn];
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[i];j++) {
double value = elementByColumn[j];
if (value) {
int iRow = row[j];
indexRowU[numberElements]=iRow;
rowCount[iRow]++;
elementU[numberElements++]=value*scale*rowScale[iRow];
}
}
start[i+1]=numberElements;
columnCount[i]=numberElements-start[i];
}
}
}
}
// Creates scales for column copy (rowCopy in model may be modified)
int
ClpPackedMatrix::scale(ClpModel * model) const
{
#ifndef NDEBUG
//checkFlags();
#endif
int numberRows = model->numberRows();
int numberColumns = matrix_->getNumCols();
// If empty - return as sanityCheck will trap
if (!numberRows||!numberColumns)
return 1;
ClpMatrixBase * rowCopyBase=model->rowCopy();
if (!rowCopyBase) {
// temporary copy
rowCopyBase = reverseOrderedCopy();
}
#ifndef NO_RTTI
ClpPackedMatrix* rowCopy =
dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
// Make sure it is really a ClpPackedMatrix
assert (rowCopy!=NULL);
#else
ClpPackedMatrix* rowCopy =
static_cast< ClpPackedMatrix*>(rowCopyBase);
#endif
const int * column = rowCopy->getIndices();
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
const double * element = rowCopy->getElements();
int scaleLength = ((model->specialOptions()&131072)==0) ? 1 : 2;
double * rowScale = new double [numberRows*scaleLength];
double * columnScale = new double [numberColumns*scaleLength];
// we are going to mark bits we are interested in
char * usefulRow = new char [numberRows];
char * usefulColumn = new char [numberColumns];
double * rowLower = model->rowLower();
double * rowUpper = model->rowUpper();
double * columnLower = model->columnLower();
double * columnUpper = model->columnUpper();
int iColumn, iRow;
//#define LEAVE_FIXED
// mark free rows
for (iRow=0;iRow<numberRows;iRow++) {
#ifndef LEAVE_FIXED
usefulRow[iRow]=0;
if (rowUpper[iRow]<1.0e20||
rowLower[iRow]>-1.0e20)
#endif
usefulRow[iRow]=1;
}
// mark empty and fixed columns
// also see if worth scaling
assert (model->scalingFlag()<4); // dynamic not implemented
double largest=0.0;
double smallest=1.0e50;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
char useful=0;
#ifndef LEAVE_FIXED
if (columnUpper[iColumn]>
columnLower[iColumn]+1.0e-12) {
#endif
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
if(elementByColumn[j]&&usefulRow[iRow]) {
useful=1;
largest = CoinMax(largest,fabs(elementByColumn[j]));
smallest = CoinMin(smallest,fabs(elementByColumn[j]));
}
}
#ifndef LEAVE_FIXED
}
#endif
usefulColumn[iColumn]=useful;
}
model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL,*model->messagesPointer())
<<smallest<<largest
<<CoinMessageEol;
if (smallest>=0.5&&largest<=2.0) {
// don't bother scaling
model->messageHandler()->message(CLP_PACKEDSCALE_FORGET,*model->messagesPointer())
<<CoinMessageEol;
delete [] rowScale;
delete [] usefulRow;
delete [] columnScale;
delete [] usefulColumn;
if (!model->rowCopy())
delete rowCopyBase; // get rid of temporary row copy
return 1;
} else {
// need to scale
int scalingMethod = model->scalingFlag();
// and see if there any empty rows
for (iRow=0;iRow<numberRows;iRow++) {
if (usefulRow[iRow]) {
CoinBigIndex j;
int useful=0;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
if (usefulColumn[iColumn]) {
useful=1;
break;
}
}
usefulRow[iRow]=useful;
}
}
double savedOverallRatio=0.0;
double tolerance = 5.0*model->primalTolerance();
double overallLargest=-1.0e-20;
double overallSmallest=1.0e20;
bool finished=false;
// if scalingMethod 3 then may change
while (!finished) {
ClpFillN ( rowScale, numberRows,1.0);
ClpFillN ( columnScale, numberColumns,1.0);
overallLargest=-1.0e-20;
overallSmallest=1.0e20;
if (scalingMethod==1||scalingMethod==3) {
// Maximum in each row
for (iRow=0;iRow<numberRows;iRow++) {
if (usefulRow[iRow]) {
CoinBigIndex j;
largest=1.0e-10;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
if (usefulColumn[iColumn]) {
double value = fabs(element[j]);
largest = CoinMax(largest,value);
assert (largest<1.0e40);
}
}
rowScale[iRow]=1.0/largest;
overallLargest = CoinMax(overallLargest,largest);
overallSmallest = CoinMin(overallSmallest,largest);
}
}
} else {
int numberPass=3;
#ifdef USE_OBJECTIVE
// This will be used to help get scale factors
double * objective = new double[numberColumns];
memcpy(objective,model->costRegion(1),numberColumns*sizeof(double));
double objScale=1.0;
#endif
while (numberPass) {
overallLargest=0.0;
overallSmallest=1.0e50;
numberPass--;
// Geometric mean on row scales
for (iRow=0;iRow<numberRows;iRow++) {
if (usefulRow[iRow]) {
CoinBigIndex j;
largest=1.0e-20;
smallest=1.0e50;
for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
int iColumn = column[j];
if (usefulColumn[iColumn]) {
double value = fabs(element[j]);
// Don't bother with tiny elements
if (value>1.0e-20) {
value *= columnScale[iColumn];
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
}
}
rowScale[iRow]=1.0/sqrt(smallest*largest);
rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
overallLargest = CoinMax(largest*rowScale[iRow],overallLargest);
overallSmallest = CoinMin(smallest*rowScale[iRow],overallSmallest);
}
}
#ifdef USE_OBJECTIVE
largest=1.0e-20;
smallest=1.0e50;
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (usefulColumn[iColumn]) {
double value = fabs(objective[iColumn]);
// Don't bother with tiny elements
if (value>1.0e-20) {
value *= columnScale[iColumn];
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
}
}
objScale=1.0/sqrt(smallest*largest);
#endif
model->messageHandler()->message(CLP_PACKEDSCALE_WHILE,*model->messagesPointer())
<<overallSmallest
<<overallLargest
<<CoinMessageEol;
// skip last column round
if (numberPass==1)
break;
// Geometric mean on column scales
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (usefulColumn[iColumn]) {
CoinBigIndex j;
largest=1.0e-20;
smallest=1.0e50;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
double value = fabs(elementByColumn[j]);
// Don't bother with tiny elements
if (value>1.0e-20&&usefulRow[iRow]) {
value *= rowScale[iRow];
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
}
#ifdef USE_OBJECTIVE
if (fabs(objective[iColumn])>1.0e-20) {
double value = fabs(objective[iColumn])*objScale;
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
#endif
columnScale[iColumn]=1.0/sqrt(smallest*largest);
columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
}
}
}
#ifdef USE_OBJECTIVE
delete [] objective;
printf("obj scale %g - use it if you want\n",objScale);
#endif
}
// If ranges will make horrid then scale
for (iRow=0;iRow<numberRows;iRow++) {
if (usefulRow[iRow]) {
double difference = rowUpper[iRow]-rowLower[iRow];
double scaledDifference = difference*rowScale[iRow];
if (scaledDifference>tolerance&&scaledDifference<1.0e-4) {
// make gap larger
rowScale[iRow] *= 1.0e-4/scaledDifference;
rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
//printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
// scaledDifference,difference*rowScale[iRow]);
}
}
}
// final pass to scale columns so largest is reasonable
// See what smallest will be if largest is 1.0
overallSmallest=1.0e50;
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (usefulColumn[iColumn]) {
CoinBigIndex j;
largest=1.0e-20;
smallest=1.0e50;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
if(elementByColumn[j]&&usefulRow[iRow]) {
double value = fabs(elementByColumn[j]*rowScale[iRow]);
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
}
if (overallSmallest*largest>smallest)
overallSmallest = smallest/largest;
}
}
if (scalingMethod==1||scalingMethod==2) {
finished=true;
} else if (savedOverallRatio==0.0&&scalingMethod!=4) {
savedOverallRatio=overallSmallest;
scalingMethod=4;
} else {
assert (scalingMethod==4);
if (overallSmallest>2.0*savedOverallRatio)
finished=true; // geometric was better
else
scalingMethod=1; // redo equilibrium
#if 0
if (model->logLevel()>2) {
if (finished)
printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
savedOverallRatio,overallSmallest);
else
printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
savedOverallRatio,overallSmallest);
}
#endif
}
}
//#define RANDOMIZE
#ifdef RANDOMIZE
// randomize by up to 10%
for (iRow=0;iRow<numberRows;iRow++) {
double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
rowScale[iRow] *= (1.0+0.1*value);
}
#endif
overallLargest=1.0;
if (overallSmallest<1.0e-1)
overallLargest = 1.0/sqrt(overallSmallest);
overallLargest = CoinMin(100.0,overallLargest);
overallSmallest=1.0e50;
//printf("scaling %d\n",model->scalingFlag());
for (iColumn=0;iColumn<numberColumns;iColumn++) {
if (columnUpper[iColumn]>
columnLower[iColumn]+1.0e-12) {
//if (usefulColumn[iColumn]) {
CoinBigIndex j;
largest=1.0e-20;
smallest=1.0e50;
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
iRow=row[j];
if(elementByColumn[j]&&usefulRow[iRow]) {
double value = fabs(elementByColumn[j]*rowScale[iRow]);
largest = CoinMax(largest,value);
smallest = CoinMin(smallest,value);
}
}
columnScale[iColumn]=overallLargest/largest;
columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
#ifdef RANDOMIZE
double value = 0.5-CoinDrand48();//between -0.5 to + 0.5
columnScale[iColumn] *= (1.0+0.1*value);
#endif
double difference = columnUpper[iColumn]-columnLower[iColumn];
if (difference<1.0e-5*columnScale[iColumn]) {
// make gap larger
columnScale[iColumn] = difference/1.0e-5;
//printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
// scaledDifference,difference*columnScale[iColumn]);
}
overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
}
}
model->messageHandler()->message(CLP_PACKEDSCALE_FINAL,*model->messagesPointer())
<<overallSmallest
<<overallLargest
<<CoinMessageEol;
delete [] usefulRow;
delete [] usefulColumn;
#ifndef SLIM_CLP
// If quadratic then make symmetric
ClpObjective * obj = model->objectiveAsObject();
#ifndef NO_RTTI
ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
#else
ClpQuadraticObjective * quadraticObj = NULL;
if (obj->type()==2)
quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
#endif
if (quadraticObj) {
CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
int numberXColumns = quadratic->getNumCols();
if (numberXColumns<numberColumns) {
// we assume symmetric
int numberQuadraticColumns=0;
int i;
//const int * columnQuadratic = quadratic->getIndices();
//const int * columnQuadraticStart = quadratic->getVectorStarts();
const int * columnQuadraticLength = quadratic->getVectorLengths();
for (i=0;i<numberXColumns;i++) {
int length=columnQuadraticLength[i];
#ifndef CORRECT_COLUMN_COUNTS
length=1;
#endif
if (length)
numberQuadraticColumns++;
}
int numberXRows = numberRows-numberQuadraticColumns;
numberQuadraticColumns=0;
for (i=0;i<numberXColumns;i++) {
int length=columnQuadraticLength[i];
#ifndef CORRECT_COLUMN_COUNTS
length=1;
#endif
if (length) {
rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
numberQuadraticColumns++;
}
}
int numberQuadraticRows=0;
for (i=0;i<numberXRows;i++) {
// See if any in row quadratic
CoinBigIndex j;
int numberQ=0;
for (j=rowStart[i];j<rowStart[i+1];j++) {
int iColumn = column[j];
if (columnQuadraticLength[iColumn])
numberQ++;
}
#ifndef CORRECT_ROW_COUNTS
numberQ=1;
#endif
if (numberQ) {
columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
numberQuadraticRows++;
}
}
// and make sure Sj okay
for (iColumn=numberQuadraticRows+numberXColumns;iColumn<numberColumns;iColumn++) {
CoinBigIndex j=columnStart[iColumn];
assert(columnLength[iColumn]==1);
int iRow=row[j];
double value = fabs(elementByColumn[j]*rowScale[iRow]);
columnScale[iColumn]=1.0/value;
}
}
}
#endif
if (scaleLength>1) {
// make copy
double * inverseScale = rowScale + numberRows;
for (iRow=0;iRow<numberRows;iRow++)
inverseScale[iRow] = 1.0/rowScale[iRow] ;
inverseScale = columnScale + numberColumns;
for (iColumn=0;iColumn<numberColumns;iColumn++)
inverseScale[iColumn] = 1.0/columnScale[iColumn] ;
}
model->setRowScale(rowScale);
model->setColumnScale(columnScale);
if (model->rowCopy()) {
// need to replace row by row
double * newElement = new double[numberColumns];
// scale row copy
for (iRow=0;iRow<numberRows;iRow++) {
CoinBigIndex j;
double scale = rowScale[iRow];
const double * elementsInThisRow = element + rowStart[iRow];
const int * columnsInThisRow = column + rowStart[iRow];
int number = rowStart[iRow+1]-rowStart[iRow];
assert (number<=numberColumns);
for (j=0;j<number;j++) {
int iColumn = columnsInThisRow[j];
newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
}
rowCopy->replaceVector(iRow,number,newElement);
}
delete [] newElement;
} else {
// no row copy
delete rowCopyBase;
}
return 0;
}
}
/* Unpacks a column into an CoinIndexedvector
*/
void
ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
int iColumn) const
{
const double * rowScale = model->rowScale();
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
CoinBigIndex i;
if (!rowScale) {
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
rowArray->add(row[i],elementByColumn[i]);
}
} else {
// apply scaling
double scale = model->columnScale()[iColumn];
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
}
}
}
/* Unpacks a column into a CoinIndexedVector
** in packed format
Note that model is NOT const. Bounds and objective could
be modified if doing column generation (just for this variable) */
void
ClpPackedMatrix::unpackPacked(ClpSimplex * model,
CoinIndexedVector * rowArray,
int iColumn) const
{
const double * rowScale = model->rowScale();
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
CoinBigIndex i;
int * index = rowArray->getIndices();
double * array = rowArray->denseVector();
int number = 0;
if (!rowScale) {
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
double value = elementByColumn[i];
if (value) {
array[number]=value;
index[number++]=iRow;
}
}
rowArray->setNumElements(number);
rowArray->setPackedMode(true);
} else {
// apply scaling
double scale = model->columnScale()[iColumn];
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
double value = elementByColumn[i]*scale*rowScale[iRow];
if (value) {
array[number]=value;
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
ClpPackedMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
int iColumn, double multiplier) const
{
const double * rowScale = model->rowScale();
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
CoinBigIndex i;
if (!rowScale) {
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
rowArray->quickAdd(iRow,multiplier*elementByColumn[i]);
}
} else {
// apply scaling
double scale = model->columnScale()[iColumn]*multiplier;
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
rowArray->quickAdd(iRow,elementByColumn[i]*scale*rowScale[iRow]);
}
}
}
/* Adds multiple of a column into an array */
void
ClpPackedMatrix::add(const ClpSimplex * model,double * array,
int iColumn, double multiplier) const
{
const double * rowScale = model->rowScale();
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
CoinBigIndex i;
if (!rowScale) {
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
array[iRow] += multiplier*elementByColumn[i];
}
} else {
// apply scaling
double scale = model->columnScale()[iColumn]*multiplier;
for (i=columnStart[iColumn];
i<columnStart[iColumn]+columnLength[iColumn];i++) {
int iRow = row[i];
array[iRow] += elementByColumn[i]*scale*rowScale[iRow];
}
}
}
/* Checks if all elements are in valid range. Can just
return true if you are not paranoid. For Clp I will
probably expect no zeros. Code can modify matrix to get rid of
small elements.
check bits (can be turned off to save time) :
1 - check if matrix has gaps
2 - check if zero elements
4 - check and compress duplicates
8 - report on large and small
*/
bool
ClpPackedMatrix::allElementsInRange(ClpModel * model,
double smallest, double largest,
int check)
{
int iColumn;
// make sure matrix correct size
matrix_->setDimensions(model->numberRows(),model->numberColumns());
CoinBigIndex numberLarge=0;;
CoinBigIndex numberSmall=0;;
CoinBigIndex numberDuplicate=0;;
int firstBadColumn=-1;
int firstBadRow=-1;
double firstBadElement=0.0;
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
int numberRows = matrix_->getNumRows();
int numberColumns = matrix_->getNumCols();
// Say no gaps
flags_ &= ~2;
if (check==14||check==10) {
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex start = columnStart[iColumn];
CoinBigIndex end = start + columnLength[iColumn];
if (end!=columnStart[iColumn+1]) {
flags_ |= 2;
break;
}
}
return true;
}
assert (check==15||check==11);
if (check==15) {
int * mark = new int [numberRows];
int i;
for (i=0;i<numberRows;i++)
mark[i]=-1;
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
CoinBigIndex start = columnStart[iColumn];
CoinBigIndex end = start + columnLength[iColumn];
if (end!=columnStart[iColumn+1])
flags_ |= 2;
for (j=start;j<end;j++) {
double value = fabs(elementByColumn[j]);
int iRow = row[j];
if (iRow<0||iRow>=numberRows) {
#ifndef COIN_BIG_INDEX
printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#elif COIN_BIG_INDEX==1
printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#else
printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#endif
return false;
}
if (mark[iRow]==-1) {
mark[iRow]=j;
} else {
// duplicate
numberDuplicate++;
}
//printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
if (!value)
flags_ |= 1; // there are zero elements
if (value<smallest) {
numberSmall++;
} else if (!(value<=largest)) {
numberLarge++;
if (firstBadColumn<0) {
firstBadColumn=iColumn;
firstBadRow=row[j];
firstBadElement=elementByColumn[j];
}
}
}
//clear mark
for (j=columnStart[iColumn];
j<columnStart[iColumn]+columnLength[iColumn];j++) {
int iRow = row[j];
mark[iRow]=-1;
}
}
delete [] mark;
} else {
// just check for out of range - not for duplicates
for (iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
CoinBigIndex start = columnStart[iColumn];
CoinBigIndex end = start + columnLength[iColumn];
if (end!=columnStart[iColumn+1])
flags_ |= 2;
for (j=start;j<end;j++) {
double value = fabs(elementByColumn[j]);
int iRow = row[j];
if (iRow<0||iRow>=numberRows) {
#ifndef COIN_BIG_INDEX
printf("Out of range %d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#elif COIN_BIG_INDEX==1
printf("Out of range %d %ld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#else
printf("Out of range %d %lld %d %g\n",iColumn,j,row[j],elementByColumn[j]);
#endif
return false;
}
if (!value)
flags_ |= 1; // there are zero elements
if (value<smallest) {
numberSmall++;
} else if (!(value<=largest)) {
numberLarge++;
if (firstBadColumn<0) {
firstBadColumn=iColumn;
firstBadRow=iRow;
firstBadElement=value;
}
}
}
}
}
if (numberLarge) {
model->messageHandler()->message(CLP_BAD_MATRIX,model->messages())
<<numberLarge
<<firstBadColumn<<firstBadRow<<firstBadElement
<<CoinMessageEol;
return false;
}
if (numberSmall)
model->messageHandler()->message(CLP_SMALLELEMENTS,model->messages())
<<numberSmall
<<CoinMessageEol;
if (numberDuplicate)
model->messageHandler()->message(CLP_DUPLICATEELEMENTS,model->messages())
<<numberDuplicate
<<CoinMessageEol;
if (numberDuplicate)
matrix_->eliminateDuplicates(smallest);
else if (numberSmall)
matrix_->compress(smallest);
// If smallest >0.0 then there can't be zero elements
if (smallest>0.0)
flags_ &= ~1;;
return true;
}
/* Given positive integer weights for each row fills in sum of weights
for each column (and slack).
Returns weights vector
*/
CoinBigIndex *
ClpPackedMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
{
int numberRows = model->numberRows();
int numberColumns = matrix_->getNumCols();
int number = numberRows+numberColumns;
CoinBigIndex * weights = new CoinBigIndex[number];
// get matrix data pointers
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
int i;
for (i=0;i<numberColumns;i++) {
CoinBigIndex j;
CoinBigIndex count=0;
for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
int iRow=row[j];
count += inputWeights[iRow];
}
weights[i]=count;
}
for (i=0;i<numberRows;i++) {
weights[i+numberColumns]=inputWeights[i];
}
return weights;
}
/* Returns largest and smallest elements of both signs.
Largest refers to largest absolute value.
*/
void
ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
double & smallestPositive, double & largestPositive)
{
smallestNegative=-COIN_DBL_MAX;
largestNegative=0.0;
smallestPositive=COIN_DBL_MAX;
largestPositive=0.0;
// get matrix data pointers
const double * elementByColumn = matrix_->getElements();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
int numberColumns = matrix_->getNumCols();
int i;
for (i=0;i<numberColumns;i++) {
CoinBigIndex j;
for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
double value = elementByColumn[j];
if (value>0.0) {
smallestPositive = CoinMin(smallestPositive,value);
largestPositive = CoinMax(largestPositive,value);
} else if (value<0.0) {
smallestNegative = CoinMax(smallestNegative,value);
largestNegative = CoinMin(largestNegative,value);
}
}
}
}
// Says whether it can do partial pricing
bool
ClpPackedMatrix::canDoPartialPricing() const
{
return true;
}
// Partial pricing
void
ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
int & bestSequence, int & numberWanted)
{
numberWanted=currentWanted_;
int start = (int) (startFraction*numberActiveColumns_);
int end = CoinMin((int) (endFraction*numberActiveColumns_+1),numberActiveColumns_);
const double * element =matrix_->getElements();
const int * row = matrix_->getIndices();
const CoinBigIndex * startColumn = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
const double * rowScale = model->rowScale();
const double * columnScale = model->columnScale();
int iSequence;
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(model->clpMatrix()->reducedCost(model,bestSequence));
else
bestDj=tolerance;
int sequenceOut = model->sequenceOut();
int saveSequence = bestSequence;
int lastScan = minimumObjectsScan_<0 ? end : start+minimumObjectsScan_;
int minNeg = minimumGoodReducedCosts_==-1 ? numberWanted : minimumGoodReducedCosts_;
if (rowScale) {
// scaled
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=0.0;
// scaled
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j]*rowScale[jRow];
}
value = fabs(cost[iSequence] +value*columnScale[iSequence]);
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=0.0;
// scaled
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j]*rowScale[jRow];
}
value = cost[iSequence] +value*columnScale[iSequence];
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=0.0;
// scaled
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j]*rowScale[jRow];
}
value = -(cost[iSequence] +value*columnScale[iSequence]);
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+minNeg<originalWanted_&&iSequence>lastScan) {
// give up
break;
}
if (!numberWanted)
break;
}
if (bestSequence!=saveSequence) {
// recompute dj
double value=0.0;
// scaled
for (j=startColumn[bestSequence];
j<startColumn[bestSequence]+length[bestSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j]*rowScale[jRow];
}
reducedCost[bestSequence] = cost[bestSequence] +value*columnScale[bestSequence];
savedBestSequence_ = bestSequence;
savedBestDj_ = reducedCost[savedBestSequence_];
}
} else {
// not scaled
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];
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j];
}
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];
// scaled
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j];
}
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];
for (j=startColumn[iSequence];
j<startColumn[iSequence]+length[iSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j];
}
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+minNeg<originalWanted_&&iSequence>lastScan) {
// give up
break;
}
if (!numberWanted)
break;
}
if (bestSequence!=saveSequence) {
// recompute dj
double value=cost[bestSequence];
for (j=startColumn[bestSequence];
j<startColumn[bestSequence]+length[bestSequence];j++) {
int jRow=row[j];
value -= duals[jRow]*element[j];
}
reducedCost[bestSequence] = value;
savedBestSequence_ = bestSequence;
savedBestDj_ = reducedCost[savedBestSequence_];
}
}
currentWanted_=numberWanted;
}
// Sets up an effective RHS
void
ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
{
delete [] rhsOffset_;
int numberRows = model->numberRows();
rhsOffset_ = new double[numberRows];
rhsOffset(model,true);
}
// Gets rid of special copies
void
ClpPackedMatrix::clearCopies()
{
delete rowCopy_;
delete columnCopy_;
rowCopy_=NULL;
columnCopy_=NULL;
flags_ &= ~(4+8);
}
// makes sure active columns correct
int
ClpPackedMatrix::refresh(ClpSimplex * model)
{
numberActiveColumns_ = matrix_->getNumCols();
#if 0
ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
ClpPackedMatrix* rowCopy =
dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
// Make sure it is really a ClpPackedMatrix
assert (rowCopy!=NULL);
const int * column = rowCopy->getIndices();
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
const int * rowLength = rowCopy->getVectorLengths();
const double * element = rowCopy->getElements();
for (int i=0;i<rowCopy->getNumRows();i++) {
if (!rowLength[i])
printf("zero row %d\n",i);
}
delete rowCopy;
#endif
return 0;
}
/* Scales rowCopy if column copy scaled
Only called if scales already exist */
void
ClpPackedMatrix::scaleRowCopy(ClpModel * model) const
{
if (model->rowCopy()) {
// need to replace row by row
int numberRows = model->numberRows();
int numberColumns = matrix_->getNumCols();
double * newElement = new double[numberColumns];
ClpMatrixBase * rowCopyBase=model->rowCopy();
#ifndef NO_RTTI
ClpPackedMatrix* rowCopy =
dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
// Make sure it is really a ClpPackedMatrix
assert (rowCopy!=NULL);
#else
ClpPackedMatrix* rowCopy =
static_cast< ClpPackedMatrix*>(rowCopyBase);
#endif
const int * column = rowCopy->getIndices();
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
const double * element = rowCopy->getElements();
const double * rowScale = model->rowScale();
const double * columnScale = model->columnScale();
// scale row copy
for (int iRow=0;iRow<numberRows;iRow++) {
CoinBigIndex j;
double scale = rowScale[iRow];
const double * elementsInThisRow = element + rowStart[iRow];
const int * columnsInThisRow = column + rowStart[iRow];
int number = rowStart[iRow+1]-rowStart[iRow];
assert (number<=numberColumns);
for (j=0;j<number;j++) {
int iColumn = columnsInThisRow[j];
newElement[j] = elementsInThisRow[j]*scale*columnScale[iColumn];
}
rowCopy->replaceVector(iRow,number,newElement);
}
delete [] newElement;
}
}
/* Realy really scales column copy
Only called if scales already exist.
Up to user ro delete */
ClpMatrixBase *
ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const
{
// need to replace column by column
int numberRows = model->numberRows();
int numberColumns = matrix_->getNumCols();
double * newElement = new double[numberRows];
ClpPackedMatrix * copy = new ClpPackedMatrix(*this);
const int * row = copy->getIndices();
const CoinBigIndex * columnStart = copy->getVectorStarts();
const int * length = copy->getVectorLengths();
const double * element = copy->getElements();
const double * rowScale = model->rowScale();
const double * columnScale = model->columnScale();
// scale column copy
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
double scale = columnScale[iColumn];
const double * elementsInThisColumn = element + columnStart[iColumn];
const int * rowsInThisColumn = row + columnStart[iColumn];
int number = length[iColumn];
assert (number<=numberRows);
for (j=0;j<number;j++) {
int iRow = rowsInThisColumn[j];
newElement[j] = elementsInThisColumn[j]*scale*rowScale[iRow];
}
copy->replaceVector(iColumn,number,newElement);
}
delete [] newElement;
return copy;
}
// Really scale matrix
void
ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
{
clearCopies();
int numberColumns = matrix_->getNumCols();
const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * length = matrix_->getVectorLengths();
double * element = matrix_->getMutableElements();
// scale
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
CoinBigIndex j;
double scale = columnScale[iColumn];
for (j=columnStart[iColumn];j<columnStart[iColumn]+length[iColumn];j++) {
int iRow = row[j];
element[j] *= scale*rowScale[iRow];
}
}
}
/* Delete the columns whose indices are listed in <code>indDel</code>. */
void
ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
{
if (matrix_->getNumCols())
matrix_->deleteCols(numDel,indDel);
clearCopies();
numberActiveColumns_ = matrix_->getNumCols();
// may now have gaps
flags_ |= 2;
matrix_->setExtraGap(1.0e-50);
}
/* Delete the rows whose indices are listed in <code>indDel</code>. */
void
ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
{
if (matrix_->getNumRows())
matrix_->deleteRows(numDel,indDel);
clearCopies();
numberActiveColumns_ = matrix_->getNumCols();
// may now have gaps
flags_ |= 2;
matrix_->setExtraGap(1.0e-50);
}
#ifndef CLP_NO_VECTOR
// Append Columns
void
ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
{
matrix_->appendCols(number,columns);
numberActiveColumns_ = matrix_->getNumCols();
clearCopies();
}
// Append Rows
void
ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
{
matrix_->appendRows(number,rows);
numberActiveColumns_ = matrix_->getNumCols();
// may now have gaps
flags_ |= 2;
clearCopies();
}
#endif
/* 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
ClpPackedMatrix::setDimensions(int numrows, int numcols)
{
matrix_->setDimensions(numrows,numcols);
}
/* 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
ClpPackedMatrix::appendMatrix(int number, int type,
const CoinBigIndex * starts, const int * index,
const double * element, int numberOther)
{
int numberErrors=0;
// make sure other dimension is big enough
if (type==0) {
// rows
if (matrix_->isColOrdered()&&numberOther>matrix_->getNumCols())
matrix_->setDimensions(-1,numberOther);
numberErrors=matrix_->appendRows(number,starts,index,element,numberOther);
} else {
// columns
if (!matrix_->isColOrdered()&&numberOther>matrix_->getNumRows())
matrix_->setDimensions(numberOther,-1);
numberErrors=matrix_->appendCols(number,starts,index,element,numberOther);
}
clearCopies();
numberActiveColumns_ = matrix_->getNumCols();
return numberErrors;
}
void
ClpPackedMatrix::specialRowCopy(ClpSimplex * model,const ClpMatrixBase * rowCopy)
{
delete rowCopy_;
rowCopy_ = new ClpPackedMatrix2(model,rowCopy->getPackedMatrix());
// See if anything in it
if (!rowCopy_->usefulInfo()) {
delete rowCopy_;
rowCopy_=NULL;
flags_ &= ~4;
} else {
flags_ |= 4;
}
}
void
ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
{
delete columnCopy_;
if ((flags_&8)!=0)
columnCopy_ = new ClpPackedMatrix3(model,matrix_);
else
columnCopy_=NULL;
}
// Say we don't want special column copy
void
ClpPackedMatrix::releaseSpecialColumnCopy()
{
flags_ &= ~8;
delete columnCopy_;
columnCopy_=NULL;
}
// Correct sequence in and out to give true value
void
ClpPackedMatrix::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
{
if (columnCopy_) {
if (sequenceIn!=-999) {
if (sequenceIn!=sequenceOut) {
if (sequenceIn<numberActiveColumns_)
columnCopy_->swapOne(model,this,sequenceIn);
if (sequenceOut<numberActiveColumns_)
columnCopy_->swapOne(model,this,sequenceOut);
}
} else {
// do all
columnCopy_->sortBlocks(model);
}
}
}
// Check validity
void
ClpPackedMatrix::checkFlags() const
{
int iColumn;
// get matrix data pointers
//const int * row = matrix_->getIndices();
const CoinBigIndex * columnStart = matrix_->getVectorStarts();
const int * columnLength = matrix_->getVectorLengths();
const double * elementByColumn = matrix_->getElements();
if (!zeros()) {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
CoinBigIndex j;
for (j=columnStart[iColumn];j<columnStart[iColumn]+columnLength[iColumn];j++) {
if (!elementByColumn[j])
abort();
}
}
}
if ((flags_&2)==0) {
for (iColumn=0;iColumn<numberActiveColumns_;iColumn++) {
if (columnStart[iColumn+1]!=columnStart[iColumn]+columnLength[iColumn]) {
abort();
}
}
}
}
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpPackedMatrix2::ClpPackedMatrix2 ()
: numberBlocks_(0),
numberRows_(0),
offset_(NULL),
count_(NULL),
rowStart_(NULL),
column_(NULL),
work_(NULL)
{
#ifdef THREAD
threadId_ = NULL;
info_ = NULL;
#endif
}
//-------------------------------------------------------------------
// Useful Constructor
//-------------------------------------------------------------------
ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * model,const CoinPackedMatrix * rowCopy)
: numberBlocks_(0),
numberRows_(0),
offset_(NULL),
count_(NULL),
rowStart_(NULL),
column_(NULL),
work_(NULL)
{
#ifdef THREAD
threadId_ = NULL;
info_ = NULL;
#endif
numberRows_ = rowCopy->getNumRows();
if (!numberRows_)
return;
int numberColumns = rowCopy->getNumCols();
const int * column = rowCopy->getIndices();
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
const int * length = rowCopy->getVectorLengths();
const double * element = rowCopy->getElements();
int chunk = 32768; // tune
//chunk=100;
// tune
#if 0
int chunkY[7]={1024,2048,4096,8192,16384,32768,65535};
int its = model->maximumIterations();
if (its>=1000000&&its<1000999) {
its -= 1000000;
its = its/10;
if (its>=7) abort();
chunk=chunkY[its];
printf("chunk size %d\n",chunk);
#define cpuid(func,ax,bx,cx,dx)\
__asm__ __volatile__ ("cpuid":\
"=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
unsigned int a,b,c,d;
int func=0;
cpuid(func,a,b,c,d);
{
int i;
unsigned int value;
value=b;
for (i=0;i<4;i++) {
printf("%c",(value&0xff));
value = value >>8;
}
value=d;
for (i=0;i<4;i++) {
printf("%c",(value&0xff));
value = value >>8;
}
value=c;
for (i=0;i<4;i++) {
printf("%c",(value&0xff));
value = value >>8;
}
printf("\n");
int maxfunc = a;
if (maxfunc>10) {
printf("not intel?\n");
abort();
}
for (func=1;func<=maxfunc;func++) {
cpuid(func,a,b,c,d);
printf("func %d, %x %x %x %x\n",func,a,b,c,d);
}
}
#else
if (numberColumns>10000||chunk==100) {
#endif
} else {
//printf("no chunk\n");
return;
}
// Could also analyze matrix to get natural breaks
numberBlocks_ = (numberColumns+chunk-1)/chunk;
#ifdef THREAD
// Get work areas
threadId_ = new pthread_t [numberBlocks_];
info_ = new dualColumn0Struct[numberBlocks_];
#endif
// Even out
chunk = (numberColumns+numberBlocks_-1)/numberBlocks_;
offset_ = new int[numberBlocks_+1];
offset_[numberBlocks_]=numberColumns;
int nRow = numberBlocks_*numberRows_;
count_ = new unsigned short[nRow];
memset(count_,0,nRow*sizeof(unsigned short));
rowStart_ = new CoinBigIndex[nRow+numberRows_+1];
CoinBigIndex nElement = rowStart[numberRows_];
rowStart_[nRow+numberRows_]=nElement;
column_ = new unsigned short [nElement];
// assumes int <= double
int sizeWork = 6*numberBlocks_;
work_ = new double[sizeWork];;
int iBlock;
int nZero=0;
for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
int start = iBlock*chunk;
offset_[iBlock]=start;
int end = start+chunk;
for (int iRow=0;iRow<numberRows_;iRow++) {
if (rowStart[iRow+1]!=rowStart[iRow]+length[iRow]) {
printf("not packed correctly - gaps\n");
abort();
}
bool lastFound=false;
int nFound=0;
for (CoinBigIndex j = rowStart[iRow];
j<rowStart[iRow]+length[iRow];j++) {
int iColumn = column[j];
if (iColumn>=start) {
if (iColumn<end) {
if (!element[j]) {
printf("not packed correctly - zero element\n");
abort();
}
column_[j]=iColumn-start;
nFound++;
if (lastFound) {
printf("not packed correctly - out of order\n");
abort();
}
} else {
//can't find any more
lastFound=true;
}
}
}
count_[iRow*numberBlocks_+iBlock]=nFound;
if (!nFound)
nZero++;
}
}
//double fraction = ((double) nZero)/((double) (numberBlocks_*numberRows_));
//printf("%d empty blocks, %g%%\n",nZero,100.0*fraction);
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
: numberBlocks_(rhs.numberBlocks_),
numberRows_(rhs.numberRows_)
{
if (numberBlocks_) {
offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
int nRow = numberBlocks_*numberRows_;
count_ = CoinCopyOfArray(rhs.count_,nRow);
rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
CoinBigIndex nElement = rowStart_[nRow+numberRows_];
column_ = CoinCopyOfArray(rhs.column_,nElement);
int sizeWork = 6*numberBlocks_;
work_ = CoinCopyOfArray(rhs.work_,sizeWork);
#ifdef THREAD
threadId_ = new pthread_t [numberBlocks_];
info_ = new dualColumn0Struct[numberBlocks_];
#endif
} else {
offset_ = NULL;
count_ = NULL;
rowStart_ = NULL;
column_ = NULL;
work_ = NULL;
#ifdef THREAD
threadId_ = NULL;
info_ = NULL;
#endif
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpPackedMatrix2::~ClpPackedMatrix2 ()
{
delete [] offset_;
delete [] count_;
delete [] rowStart_;
delete [] column_;
delete [] work_;
#ifdef THREAD
delete [] threadId_;
delete [] info_;
#endif
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpPackedMatrix2 &
ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
{
if (this != &rhs) {
numberBlocks_ = rhs.numberBlocks_;
numberRows_ = rhs.numberRows_;
delete [] offset_;
delete [] count_;
delete [] rowStart_;
delete [] column_;
delete [] work_;
#ifdef THREAD
delete [] threadId_;
delete [] info_;
#endif
if (numberBlocks_) {
offset_ = CoinCopyOfArray(rhs.offset_,numberBlocks_+1);
int nRow = numberBlocks_*numberRows_;
count_ = CoinCopyOfArray(rhs.count_,nRow);
rowStart_ = CoinCopyOfArray(rhs.rowStart_,nRow+numberRows_+1);
CoinBigIndex nElement = rowStart_[nRow+numberRows_];
column_ = CoinCopyOfArray(rhs.column_,nElement);
int sizeWork = 6*numberBlocks_;
work_ = CoinCopyOfArray(rhs.work_,sizeWork);
#ifdef THREAD
threadId_ = new pthread_t [numberBlocks_];
info_ = new dualColumn0Struct[numberBlocks_];
#endif
} else {
offset_ = NULL;
count_ = NULL;
rowStart_ = NULL;
column_ = NULL;
work_ = NULL;
#ifdef THREAD
threadId_ = NULL;
info_ = NULL;
#endif
}
}
return *this;
}
static int dualColumn0(const ClpSimplex * model,double * spare,
int * spareIndex,const double * arrayTemp,
const int * indexTemp,int numberIn,
int offset, double acceptablePivot,double * bestPossiblePtr,
double * upperThetaPtr,int * posFreePtr, double * freePivotPtr)
{
// do dualColumn0
int i;
int numberRemaining=0;
double bestPossible=0.0;
double upperTheta=1.0e31;
double freePivot = acceptablePivot;
int posFree=-1;
const double * reducedCost=model->djRegion(1);
double dualTolerance = model->dualTolerance();
// We can also see if infeasible or pivoting on free
double tentativeTheta = 1.0e25;
for (i=0;i<numberIn;i++) {
double alpha = arrayTemp[i];
int iSequence = indexTemp[i]+offset;
double oldValue;
double value;
bool keep;
switch(model->getStatus(iSequence)) {
case ClpSimplex::basic:
case ClpSimplex::isFixed:
break;
case ClpSimplex::isFree:
case ClpSimplex::superBasic:
bestPossible = CoinMax(bestPossible,fabs(alpha));
oldValue = reducedCost[iSequence];
// If free has to be very large - should come in via dualRow
if (model->getStatus(iSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
break;
if (oldValue>dualTolerance) {
keep = true;
} else if (oldValue<-dualTolerance) {
keep = true;
} else {
if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5))
keep = true;
else
keep=false;
}
if (keep) {
// free - choose largest
if (fabs(alpha)>freePivot) {
freePivot=fabs(alpha);
posFree = i;
}
}
break;
case ClpSimplex::atUpperBound:
oldValue = reducedCost[iSequence];
value = oldValue-tentativeTheta*alpha;
//assert (oldValue<=dualTolerance*1.0001);
if (value>dualTolerance) {
bestPossible = CoinMax(bestPossible,-alpha);
value = oldValue-upperTheta*alpha;
if (value>dualTolerance && -alpha>=acceptablePivot)
upperTheta = (oldValue-dualTolerance)/alpha;
// add to list
spare[numberRemaining]=alpha;
spareIndex[numberRemaining++]=iSequence;
}
break;
case ClpSimplex::atLowerBound:
oldValue = reducedCost[iSequence];
value = oldValue-tentativeTheta*alpha;
//assert (oldValue>=-dualTolerance*1.0001);
if (value<-dualTolerance) {
bestPossible = CoinMax(bestPossible,alpha);
value = oldValue-upperTheta*alpha;
if (value<-dualTolerance && alpha>=acceptablePivot)
upperTheta = (oldValue+dualTolerance)/alpha;
// add to list
spare[numberRemaining]=alpha;
spareIndex[numberRemaining++]=iSequence;
}
break;
}
}
*bestPossiblePtr = bestPossible;
*upperThetaPtr = upperTheta;
*freePivotPtr = freePivot;
*posFreePtr = posFree;
return numberRemaining;
}
static int doOneBlock(double * array, int * index,
const double * pi,const CoinBigIndex * rowStart, const double * element,
const unsigned short * column, int numberInRowArray, int numberLook)
{
int iWhich=0;
int nextN=0;
CoinBigIndex nextStart=0;
double nextPi=0.0;
for (;iWhich<numberInRowArray;iWhich++) {
nextStart = rowStart[0];
nextN = rowStart[numberInRowArray]-nextStart;
rowStart++;
if (nextN) {
nextPi = pi[iWhich];
break;
}
}
int i;
while (iWhich<numberInRowArray) {
double value = nextPi;
CoinBigIndex j= nextStart;
int n = nextN;
// get next
iWhich++;
for (;iWhich<numberInRowArray;iWhich++) {
nextStart = rowStart[0];
nextN = rowStart[numberInRowArray]-nextStart;
rowStart++;
if (nextN) {
coin_prefetch(element+nextStart);
nextPi = pi[iWhich];
break;
}
}
CoinBigIndex end =j+n;
//coin_prefetch(element+rowStart_[i+1]);
//coin_prefetch(column_+rowStart_[i+1]);
if (n<100) {
if ((n&1)!=0) {
unsigned int jColumn = column[j];
array[jColumn] -= value*element[j];
j++;
}
coin_prefetch(column+nextStart);
for (;j<end;j+=2) {
unsigned int jColumn0 = column[j];
double value0 = value*element[j];
unsigned int jColumn1 = column[j+1];
double value1 = value*element[j+1];
array[jColumn0] -= value0;
array[jColumn1] -= value1;
}
} else {
if ((n&1)!=0) {
unsigned int jColumn = column[j];
array[jColumn] -= value*element[j];
j++;
}
if ((n&2)!=0) {
unsigned int jColumn0 = column[j];
double value0 = value*element[j];
unsigned int jColumn1 = column[j+1];
double value1 = value*element[j+1];
array[jColumn0] -= value0;
array[jColumn1] -= value1;
j+=2;
}
if ((n&4)!=0) {
unsigned int jColumn0 = column[j];
double value0 = value*element[j];
unsigned int jColumn1 = column[j+1];
double value1 = value*element[j+1];
unsigned int jColumn2 = column[j+2];
double value2 = value*element[j+2];
unsigned int jColumn3 = column[j+3];
double value3 = value*element[j+3];
array[jColumn0] -= value0;
array[jColumn1] -= value1;
array[jColumn2] -= value2;
array[jColumn3] -= value3;
j+=4;
}
//coin_prefetch(column+nextStart);
for (;j<end;j+=8) {
coin_prefetch(element+j+16);
unsigned int jColumn0 = column[j];
double value0 = value*element[j];
unsigned int jColumn1 = column[j+1];
double value1 = value*element[j+1];
unsigned int jColumn2 = column[j+2];
double value2 = value*element[j+2];
unsigned int jColumn3 = column[j+3];
double value3 = value*element[j+3];
array[jColumn0] -= value0;
array[jColumn1] -= value1;
array[jColumn2] -= value2;
array[jColumn3] -= value3;
coin_prefetch(column+j+16);
jColumn0 = column[j+4];
value0 = value*element[j+4];
jColumn1 = column[j+5];
value1 = value*element[j+5];
jColumn2 = column[j+6];
value2 = value*element[j+6];
jColumn3 = column[j+7];
value3 = value*element[j+7];
array[jColumn0] -= value0;
array[jColumn1] -= value1;
array[jColumn2] -= value2;
array[jColumn3] -= value3;
}
}
}
// get rid of tiny values
int nSmall = numberLook;
int numberNonZero=0;
for (i=0;i<nSmall;i++) {
double value = array[i];
array[i]=0.0;
if (fabs(value)>1.0e-12) {
array[numberNonZero]=value;
index[numberNonZero++]=i;
}
}
for (;i<numberLook;i+=4) {
double value0 = array[i+0];
double value1 = array[i+1];
double value2 = array[i+2];
double value3 = array[i+3];
array[i+0]=0.0;
array[i+1]=0.0;
array[i+2]=0.0;
array[i+3]=0.0;
if (fabs(value0)>1.0e-12) {
array[numberNonZero]=value0;
index[numberNonZero++]=i+0;
}
if (fabs(value1)>1.0e-12) {
array[numberNonZero]=value1;
index[numberNonZero++]=i+1;
}
if (fabs(value2)>1.0e-12) {
array[numberNonZero]=value2;
index[numberNonZero++]=i+2;
}
if (fabs(value3)>1.0e-12) {
array[numberNonZero]=value3;
index[numberNonZero++]=i+3;
}
}
return numberNonZero;
}
#ifdef THREAD
static void * doOneBlockThread(void * voidInfo)
{
dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
*(info->numberInPtr) = doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
info->rowStart,info->element,info->column,
info->numberInRowArray,info->numberLook);
return NULL;
}
static void * doOneBlockAnd0Thread(void * voidInfo)
{
dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
*(info->numberInPtr) = doOneBlock(info->arrayTemp,info->indexTemp,info->pi,
info->rowStart,info->element,info->column,
info->numberInRowArray,info->numberLook);
*(info->numberOutPtr) = dualColumn0(info->model,info->spare,
info->spareIndex,(const double *)info->arrayTemp,
(const int *) info->indexTemp,*(info->numberInPtr),
info->offset, info->acceptablePivot,info->bestPossiblePtr,
info->upperThetaPtr,info->posFreePtr, info->freePivotPtr);
return NULL;
}
#endif
/* Return <code>x * scalar * A in <code>z</code>.
Note - x packed and z will be packed mode
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix2::transposeTimes(const ClpSimplex * model,
const CoinPackedMatrix * rowCopy,
const CoinIndexedVector * rowArray,
CoinIndexedVector * spareArray,
CoinIndexedVector * columnArray) const
{
// See if dualColumn0 coding wanted
bool dualColumn = model->spareIntArray_[0]==1;
double acceptablePivot = model->spareDoubleArray_[0];
double bestPossible=0.0;
double upperTheta=1.0e31;
double freePivot = acceptablePivot;
int posFree=-1;
int numberRemaining=0;
//if (model->numberIterations()>=200000) {
//printf("time %g\n",CoinCpuTime()-startTime);
//exit(77);
//}
double * pi = rowArray->denseVector();
int numberNonZero=0;
int * index = columnArray->getIndices();
double * array = columnArray->denseVector();
int numberInRowArray = rowArray->getNumElements();
const int * whichRow = rowArray->getIndices();
double * element = const_cast<double *>(rowCopy->getElements());
const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
int i;
CoinBigIndex * rowStart2 = rowStart_;
if (!dualColumn) {
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
CoinBigIndex start = rowStart[iRow];
*rowStart2 = start;
unsigned short * count1 = count_+iRow*numberBlocks_;
int put=0;
for (int j=0;j<numberBlocks_;j++) {
put += numberInRowArray;
start += count1[j];
rowStart2[put]=start;
}
rowStart2 ++;
}
} else {
// also do dualColumn stuff
double * spare = spareArray->denseVector();
int * spareIndex = spareArray->getIndices();
const double * reducedCost=model->djRegion(0);
double dualTolerance = model->dualTolerance();
// We can also see if infeasible or pivoting on free
double tentativeTheta = 1.0e25;
int addSequence=model->numberColumns();
for (i=0;i<numberInRowArray;i++) {
int iRow = whichRow[i];
double alpha=pi[i];
double oldValue;
double value;
bool keep;
switch(model->getStatus(iRow+addSequence)) {
case ClpSimplex::basic:
case ClpSimplex::isFixed:
break;
case ClpSimplex::isFree:
case ClpSimplex::superBasic:
bestPossible = CoinMax(bestPossible,fabs(alpha));
oldValue = reducedCost[iRow];
// If free has to be very large - should come in via dualRow
if (model->getStatus(iRow+addSequence)==ClpSimplex::isFree&&fabs(alpha)<1.0e-3)
break;
if (oldValue>dualTolerance) {
keep = true;
} else if (oldValue<-dualTolerance) {
keep = true;
} else {
if (fabs(alpha)>CoinMax(10.0*acceptablePivot,1.0e-5))
keep = true;
else
keep=false;
}
if (keep) {
// free - choose largest
if (fabs(alpha)>freePivot) {
freePivot=fabs(alpha);
posFree = i + addSequence;
}
}
break;
case ClpSimplex::atUpperBound:
oldValue = reducedCost[iRow];
value = oldValue-tentativeTheta*alpha;
//assert (oldValue<=dualTolerance*1.0001);
if (value>dualTolerance) {
bestPossible = CoinMax(bestPossible,-alpha);
value = oldValue-upperTheta*alpha;
if (value>dualTolerance && -alpha>=acceptablePivot)
upperTheta = (oldValue-dualTolerance)/alpha;
// add to list
spare[numberRemaining]=alpha;
spareIndex[numberRemaining++]=iRow+addSequence;
}
break;
case ClpSimplex::atLowerBound:
oldValue = reducedCost[iRow];
value = oldValue-tentativeTheta*alpha;
//assert (oldValue>=-dualTolerance*1.0001);
if (value<-dualTolerance) {
bestPossible = CoinMax(bestPossible,alpha);
value = oldValue-upperTheta*alpha;
if (value<-dualTolerance && alpha>=acceptablePivot)
upperTheta = (oldValue+dualTolerance)/alpha;
// add to list
spare[numberRemaining]=alpha;
spareIndex[numberRemaining++]=iRow+addSequence;
}
break;
}
CoinBigIndex start = rowStart[iRow];
*rowStart2 = start;
unsigned short * count1 = count_+iRow*numberBlocks_;
int put=0;
for (int j=0;j<numberBlocks_;j++) {
put += numberInRowArray;
start += count1[j];
rowStart2[put]=start;
}
rowStart2 ++;
}
}
double * spare = spareArray->denseVector();
int * spareIndex = spareArray->getIndices();
int saveNumberRemaining=numberRemaining;
int iBlock;
for (iBlock=0;iBlock<numberBlocks_;iBlock++) {
double * dwork = work_+6*iBlock;
int * iwork = (int *) (dwork+3);
if (!dualColumn) {
#ifndef THREAD
int offset = offset_[iBlock];
int offset3 = offset;
offset = numberNonZero;
double * arrayTemp = array+offset;
int * indexTemp = index+offset;
iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
element,column_,numberInRowArray,offset_[iBlock+1]-offset);
int number = iwork[0];
for (i=0;i<number;i++) {
//double value = arrayTemp[i];
//arrayTemp[i]=0.0;
//array[numberNonZero]=value;
index[numberNonZero++]=indexTemp[i]+offset3;
}
#else
int offset = offset_[iBlock];
double * arrayTemp = array+offset;
int * indexTemp = index+offset;
dualColumn0Struct * infoPtr = info_+iBlock;
infoPtr->arrayTemp=arrayTemp;
infoPtr->indexTemp=indexTemp;
infoPtr->numberInPtr=&iwork[0];
infoPtr->pi=pi;
infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
infoPtr->element=element;
infoPtr->column=column_;
infoPtr->numberInRowArray=numberInRowArray;
infoPtr->numberLook=offset_[iBlock+1]-offset;
pthread_create(&threadId_[iBlock],NULL,doOneBlockThread,infoPtr);
#endif
} else {
#ifndef THREAD
int offset = offset_[iBlock];
// allow for already saved
int offset2 = offset + saveNumberRemaining;
int offset3 = offset;
offset = numberNonZero;
offset2 = numberRemaining;
double * arrayTemp = array+offset;
int * indexTemp = index+offset;
iwork[0] = doOneBlock(arrayTemp,indexTemp,pi,rowStart_+numberInRowArray*iBlock,
element,column_,numberInRowArray,offset_[iBlock+1]-offset);
iwork[1] = dualColumn0(model,spare+offset2,
spareIndex+offset2,
arrayTemp,indexTemp,
iwork[0],offset3,acceptablePivot,
&dwork[0],&dwork[1],&iwork[2],
&dwork[2]);
int number = iwork[0];
int numberLook = iwork[1];
#if 1
numberRemaining += numberLook;
#else
double * spareTemp = spare+offset2;
const int * spareIndexTemp = spareIndex+offset2;
for (i=0;i<numberLook;i++) {
double value = spareTemp[i];
spareTemp[i]=0.0;
spare[numberRemaining] = value;
spareIndex[numberRemaining++] = spareIndexTemp[i];
}
#endif
if (dwork[2]>freePivot) {
freePivot=dwork[2];
posFree = iwork[2]+numberNonZero;
}
upperTheta = CoinMin(dwork[1],upperTheta);
bestPossible = CoinMax(dwork[0],bestPossible);
for (i=0;i<number;i++) {
// double value = arrayTemp[i];
//arrayTemp[i]=0.0;
//array[numberNonZero]=value;
index[numberNonZero++]=indexTemp[i]+offset3;
}
#else
int offset = offset_[iBlock];
// allow for already saved
int offset2 = offset + saveNumberRemaining;
double * arrayTemp = array+offset;
int * indexTemp = index+offset;
dualColumn0Struct * infoPtr = info_+iBlock;
infoPtr->model=model;
infoPtr->spare=spare+offset2;
infoPtr->spareIndex=spareIndex+offset2;
infoPtr->arrayTemp=arrayTemp;
infoPtr->indexTemp=indexTemp;
infoPtr->numberInPtr=&iwork[0];
infoPtr->offset=offset;
infoPtr->acceptablePivot=acceptablePivot;
infoPtr->bestPossiblePtr=&dwork[0];
infoPtr->upperThetaPtr=&dwork[1];
infoPtr->posFreePtr=&iwork[2];
infoPtr->freePivotPtr=&dwork[2];
infoPtr->numberOutPtr=&iwork[1];
infoPtr->pi=pi;
infoPtr->rowStart=rowStart_+numberInRowArray*iBlock;
infoPtr->element=element;
infoPtr->column=column_;
infoPtr->numberInRowArray=numberInRowArray;
infoPtr->numberLook=offset_[iBlock+1]-offset;
if (iBlock>=2)
pthread_join(threadId_[iBlock-2],NULL);
pthread_create(threadId_+iBlock,NULL,doOneBlockAnd0Thread,infoPtr);
//pthread_join(threadId_[iBlock],NULL);
#endif
}
}
for ( iBlock=CoinMax(0,numberBlocks_-2);iBlock<numberBlocks_;iBlock++) {
#ifdef THREAD
pthread_join(threadId_[iBlock],NULL);
#endif
}
#ifdef THREAD
for ( iBlock=0;iBlock<numberBlocks_;iBlock++) {
//pthread_join(threadId_[iBlock],NULL);
int offset = offset_[iBlock];
double * dwork = work_+6*iBlock;
int * iwork = (int *) (dwork+3);
int number = iwork[0];
if (dualColumn) {
// allow for already saved
int offset2 = offset + saveNumberRemaining;
int numberLook = iwork[1];
double * spareTemp = spare+offset2;
const int * spareIndexTemp = spareIndex+offset2;
for (i=0;i<numberLook;i++) {
double value = spareTemp[i];
spareTemp[i]=0.0;
spare[numberRemaining] = value;
spareIndex[numberRemaining++] = spareIndexTemp[i];
}
if (dwork[2]>freePivot) {
freePivot=dwork[2];
posFree = iwork[2]+numberNonZero;
}
upperTheta = CoinMin(dwork[1],upperTheta);
bestPossible = CoinMax(dwork[0],bestPossible);
}
double * arrayTemp = array+offset;
const int * indexTemp = index+offset;
for (i=0;i<number;i++) {
double value = arrayTemp[i];
arrayTemp[i]=0.0;
array[numberNonZero]=value;
index[numberNonZero++]=indexTemp[i]+offset;
}
}
#endif
columnArray->setNumElements(numberNonZero);
columnArray->setPackedMode(true);
if (dualColumn) {
model->spareDoubleArray_[0]=upperTheta;
model->spareDoubleArray_[1]=bestPossible;
// and theta and alpha and sequence
if (posFree<0) {
model->spareIntArray_[1]=-1;
} else {
const double * reducedCost=model->djRegion(0);
double alpha;
int numberColumns=model->numberColumns();
if (posFree<numberColumns) {
alpha = columnArray->denseVector()[posFree];
posFree = columnArray->getIndices()[posFree];
} else {
alpha = rowArray->denseVector()[posFree-numberColumns];
posFree = rowArray->getIndices()[posFree-numberColumns]+numberColumns;
}
model->spareDoubleArray_[2]=fabs(reducedCost[posFree]/alpha);;
model->spareDoubleArray_[3]=alpha;
model->spareIntArray_[1]=posFree;
}
spareArray->setNumElements(numberRemaining);
// signal done
model->spareIntArray_[0]=-1;
}
}
/* Default constructor. */
ClpPackedMatrix3::ClpPackedMatrix3()
: numberBlocks_(0),
numberColumns_(0),
column_(NULL),
start_(NULL),
row_(NULL),
element_(NULL),
block_(NULL)
{
}
/* Constructor from copy. */
ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model,const CoinPackedMatrix * columnCopy)
: numberBlocks_(0),
numberColumns_(0),
column_(NULL),
start_(NULL),
row_(NULL),
element_(NULL),
block_(NULL)
{
#define MINBLOCK 6
#define MAXBLOCK 100
#define MAXUNROLL 10
numberColumns_ = columnCopy->getNumCols();
int numberRows = columnCopy->getNumRows();
int * counts = new int[numberRows+1];
CoinZeroN(counts,numberRows+1);
CoinBigIndex nels=0;
CoinBigIndex nZeroEl=0;
int iColumn;
// get matrix data pointers
const int * row = columnCopy->getIndices();
const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
const int * columnLength = columnCopy->getVectorLengths();
const double * elementByColumn = columnCopy->getElements();
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
CoinBigIndex start = columnStart[iColumn];
int n = columnLength[iColumn];
CoinBigIndex end = start + n;
int kZero=0;
for (CoinBigIndex j=start;j<end;j++) {
if(!elementByColumn[j])
kZero++;
}
n -= kZero;
nZeroEl += kZero;
nels += n;
counts[n]++;
}
int nZeroColumns = counts[0];
counts[0]=-1;
numberColumns_ -= nZeroColumns;
column_ = new int [2*numberColumns_];
int * lookup = column_ + numberColumns_;
row_ = new int[nels];
element_ = new double[nels];
int nOdd=0;
CoinBigIndex nInOdd=0;
int i;
for (i=1;i<=numberRows;i++) {
int n = counts[i];
if (n) {
if (n<MINBLOCK||i>MAXBLOCK) {
nOdd += n;
counts[i]=-1;
nInOdd += n*i;
} else {
numberBlocks_++;
}
} else {
counts[i]=-1;
}
}
start_ = new CoinBigIndex [nOdd+1];
// even if no blocks do a dummy one
numberBlocks_ = CoinMax(numberBlocks_,1);
block_ = (blockStruct *) new char [numberBlocks_*sizeof(blockStruct)];
memset(block_,0,numberBlocks_*sizeof(blockStruct));
// Fill in what we can
int nTotal=nOdd;
// in case no blocks
block_->startIndices_=nTotal;
nels=nInOdd;
int nBlock=0;
for (i=0;i<=CoinMin(MAXBLOCK,numberRows);i++) {
if (counts[i]>0) {
blockStruct * block = block_ + nBlock;
int n=counts[i];
counts[i]= nBlock; // backward pointer
nBlock++;
block->startIndices_ = nTotal;
block->startElements_ = nels;
block->numberElements_ = i;
// up counts
nTotal += n;
nels += n*i;
}
}
// fill
start_[0]=0;
nOdd=0;
nInOdd=0;
const double * columnScale = model->columnScale();
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
CoinBigIndex start = columnStart[iColumn];
int n = columnLength[iColumn];
CoinBigIndex end = start + n;
int kZero=0;
for (CoinBigIndex j=start;j<end;j++) {
if(!elementByColumn[j])
kZero++;
}
n -= kZero;
if (n) {
int iBlock = counts[n];
if (iBlock>=0) {
blockStruct * block = block_ + iBlock;
int k = block->numberInBlock_;
block->numberInBlock_ ++;
column_[block->startIndices_+k]=iColumn;
lookup[iColumn]=k;
CoinBigIndex put = block->startElements_+k*n;
for (CoinBigIndex j=start;j<end;j++) {
double value = elementByColumn[j];
if(value) {
if (columnScale)
value *= columnScale[iColumn];
element_[put] = value;
row_[put++]=row[j];
}
}
} else {
// odd ones
for (CoinBigIndex j=start;j<end;j++) {
double value = elementByColumn[j];
if(value) {
if (columnScale)
value *= columnScale[iColumn];
element_[nInOdd] = value;
row_[nInOdd++]=row[j];
}
}
column_[nOdd]=iColumn;
lookup[iColumn]=-1;
nOdd++;
start_[nOdd]=nInOdd;
}
}
}
delete [] counts;
}
/* Destructor */
ClpPackedMatrix3::~ClpPackedMatrix3()
{
delete [] column_;
delete [] start_;
delete [] row_;
delete [] element_;
delete [] block_;
}
/* The copy constructor. */
ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
: numberBlocks_(rhs.numberBlocks_),
numberColumns_(rhs.numberColumns_),
column_(NULL),
start_(NULL),
row_(NULL),
element_(NULL),
block_(NULL)
{
if (rhs.numberBlocks_) {
block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
int numberOdd = block_->startIndices_;
start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
blockStruct * lastBlock = block_ + (numberBlocks_-1);
CoinBigIndex numberElements = lastBlock->startElements_+
lastBlock->numberInBlock_*lastBlock->numberElements_;
row_ = CoinCopyOfArray(rhs.row_,numberElements);
element_ = CoinCopyOfArray(rhs.element_,numberElements);
}
}
ClpPackedMatrix3&
ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
{
if (this != &rhs) {
delete [] column_;
delete [] start_;
delete [] row_;
delete [] element_;
delete [] block_;
numberBlocks_ = rhs.numberBlocks_;
numberColumns_ = rhs.numberColumns_;
if (rhs.numberBlocks_) {
block_ = CoinCopyOfArray(rhs.block_,numberBlocks_);
column_ = CoinCopyOfArray(rhs.column_,2*numberColumns_);
int numberOdd = block_->startIndices_;
start_ = CoinCopyOfArray(rhs.start_,numberOdd+1);
blockStruct * lastBlock = block_ + (numberBlocks_-1);
CoinBigIndex numberElements = lastBlock->startElements_+
lastBlock->numberInBlock_*lastBlock->numberElements_;
row_ = CoinCopyOfArray(rhs.row_,numberElements);
element_ = CoinCopyOfArray(rhs.element_,numberElements);
} else {
column_ = NULL;
start_ = NULL;
row_ = NULL;
element_ = NULL;
block_ = NULL;
}
}
return *this;
}
/* Sort blocks */
void
ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
{
int * lookup = column_ + numberColumns_;
for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
blockStruct * block = block_ + iBlock;
int numberInBlock = block->numberInBlock_;
int nel = block->numberElements_;
int * row = row_ + block->startElements_;
double * element = element_ + block->startElements_;
int * column = column_ + block->startIndices_;
int lastPrice=0;
int firstNotPrice=numberInBlock-1;
while (lastPrice<=firstNotPrice) {
// find first basic or fixed
int iColumn=numberInBlock;
for (;lastPrice<=firstNotPrice;lastPrice++) {
iColumn = column[lastPrice];
if (model->getColumnStatus(iColumn)==ClpSimplex::basic||
model->getColumnStatus(iColumn)==ClpSimplex::isFixed)
break;
}
// find last non basic or fixed
int jColumn=-1;
for (;firstNotPrice>lastPrice;firstNotPrice--) {
jColumn = column[firstNotPrice];
if (model->getColumnStatus(jColumn)!=ClpSimplex::basic&&
model->getColumnStatus(jColumn)!=ClpSimplex::isFixed)
break;
}
if (firstNotPrice>lastPrice) {
assert (column[lastPrice]==iColumn);
assert (column[firstNotPrice]==jColumn);
// need to swap
column[firstNotPrice]=iColumn;
lookup[iColumn]=firstNotPrice;
column[lastPrice]=jColumn;
lookup[jColumn]=lastPrice;
double * elementA = element + lastPrice*nel;
int * rowA = row + lastPrice*nel;
double * elementB = element + firstNotPrice*nel;
int * rowB = row + firstNotPrice*nel;
for (int i=0;i<nel;i++) {
int temp = rowA[i];
double tempE = elementA[i];
rowA[i]=rowB[i];
elementA[i]=elementB[i];
rowB[i]=temp;
elementB[i]=tempE;
}
firstNotPrice--;
lastPrice++;
} else if (lastPrice==firstNotPrice) {
// make sure correct side
iColumn = column[lastPrice];
if (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
model->getColumnStatus(iColumn)!=ClpSimplex::isFixed)
lastPrice++;
break;
}
}
block->numberPrice_=lastPrice;
#ifndef NDEBUG
// check
int i;
for (i=0;i<lastPrice;i++) {
int iColumn = column[i];
assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
assert (lookup[iColumn]==i);
}
for (;i<numberInBlock;i++) {
int iColumn = column[i];
assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
assert (lookup[iColumn]==i);
}
#endif
}
}
// Swap one variable
void
ClpPackedMatrix3::swapOne(const ClpSimplex * model,const ClpPackedMatrix * matrix,
int iColumn)
{
int * lookup = column_ + numberColumns_;
// position in block
int kA = lookup[iColumn];
if (kA<0)
return; // odd one
// get matrix data pointers
const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
//const int * row = columnCopy->getIndices();
const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
const int * columnLength = columnCopy->getVectorLengths();
const double * elementByColumn = columnCopy->getElements();
CoinBigIndex start = columnStart[iColumn];
int n = columnLength[iColumn];
if (matrix->zeros()) {
CoinBigIndex end = start + n;
for (CoinBigIndex j=start;j<end;j++) {
if(!elementByColumn[j])
n--;
}
}
// find block - could do binary search
int iBlock=CoinMin(n,numberBlocks_)-1;
while (block_[iBlock].numberElements_!=n)
iBlock--;
blockStruct * block = block_ + iBlock;
int nel = block->numberElements_;
int * row = row_ + block->startElements_;
double * element = element_ + block->startElements_;
int * column = column_ + block->startIndices_;
assert (column[kA]==iColumn);
bool moveUp = (model->getColumnStatus(iColumn)==ClpSimplex::basic||
model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
int lastPrice=block->numberPrice_;
int kB;
if (moveUp) {
kB=lastPrice-1;
block->numberPrice_--;
} else {
kB=lastPrice;
block->numberPrice_++;
}
int jColumn = column[kB];
column[kA]=jColumn;
lookup[jColumn]=kA;
column[kB]=iColumn;
lookup[iColumn]=kB;
double * elementA = element + kB*nel;
int * rowA = row + kB*nel;
double * elementB = element + kA*nel;
int * rowB = row + kA*nel;
int i;
for (i=0;i<nel;i++) {
int temp = rowA[i];
double tempE = elementA[i];
rowA[i]=rowB[i];
elementA[i]=elementB[i];
rowB[i]=temp;
elementB[i]=tempE;
}
#if 1
#ifndef NDEBUG
// check
for (i=0;i<block->numberPrice_;i++) {
int iColumn = column[i];
if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
assert (model->getColumnStatus(iColumn)!=ClpSimplex::basic&&
model->getColumnStatus(iColumn)!=ClpSimplex::isFixed);
assert (lookup[iColumn]==i);
}
int numberInBlock = block->numberInBlock_;
for (;i<numberInBlock;i++) {
int iColumn = column[i];
if (iColumn!=model->sequenceIn()&&iColumn!=model->sequenceOut())
assert (model->getColumnStatus(iColumn)==ClpSimplex::basic||
model->getColumnStatus(iColumn)==ClpSimplex::isFixed);
assert (lookup[iColumn]==i);
}
#endif
#endif
}
/* Return <code>x * -1 * A in <code>z</code>.
Note - x packed and z will be packed mode
Squashes small elements and knows about ClpSimplex */
void
ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
const double * pi,
CoinIndexedVector * output) const
{
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
double zeroTolerance = model->factorization()->zeroTolerance();
double value = 0.0;
CoinBigIndex j;
int numberOdd = block_->startIndices_;
if (numberOdd) {
// A) as probably long may be worth unrolling
CoinBigIndex end = start_[1];
for (j=start_[0]; j<end;j++) {
int iRow = row_[j];
value += pi[iRow]*element_[j];
}
int iColumn;
// int jColumn=column_[0];
for (iColumn=0;iColumn<numberOdd-1;iColumn++) {
CoinBigIndex start = end;
end = start_[iColumn+2];
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=column_[iColumn];
//index[numberNonZero++]=jColumn;
}
// jColumn = column_[iColumn+1];
value = 0.0;
//if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
for (j=start; j<end;j++) {
int iRow = row_[j];
value += pi[iRow]*element_[j];
}
//}
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=column_[iColumn];
//index[numberNonZero++]=jColumn;
}
}
for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
// B) Can sort so just do nonbasic (and nonfixed)
// C) Can do two at a time (if so put odd one into start_)
// D) can use switch
blockStruct * block = block_ + iBlock;
//int numberPrice = block->numberInBlock_;
int numberPrice = block->numberPrice_;
int nel = block->numberElements_;
int * row = row_ + block->startElements_;
double * element = element_ + block->startElements_;
int * column = column_ + block->startIndices_;
#if 0
// two at a time
if ((numberPrice&1)!=0) {
double value = 0.0;
int nel2=nel;
for (; nel2;nel2--) {
int iRow = *row++;
value += pi[iRow]*(*element++);
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=*column;
}
column++;
}
numberPrice = numberPrice>>1;
switch ((nel%2)) {
// 2 k +0
case 0:
for (;numberPrice;numberPrice--) {
double value0 = 0.0;
double value1 = 0.0;
int nel2=nel;
for (; nel2;nel2--) {
int iRow0 = *row;
int iRow1 = *(row+nel);
row++;
double element0 = *element;
double element1 = *(element+nel);
element++;
value0 += pi[iRow0]*element0;
value1 += pi[iRow1]*element1;
}
row += nel;
element += nel;
if (fabs(value0)>zeroTolerance) {
array[numberNonZero]=value0;
index[numberNonZero++]=*column;
}
column++;
if (fabs(value1)>zeroTolerance) {
array[numberNonZero]=value1;
index[numberNonZero++]=*column;
}
column++;
}
break;
// 2 k +1
case 1:
for (;numberPrice;numberPrice--) {
double value0;
double value1;
int nel2=nel-1;
{
int iRow0 = row[0];
int iRow1 = row[nel];
double pi0 = pi[iRow0];
double pi1 = pi[iRow1];
value0 = pi0*element[0];
value1 = pi1*element[nel];
row++;
element++;
}
for (; nel2;nel2--) {
int iRow0 = *row;
int iRow1 = *(row+nel);
row++;
double element0 = *element;
double element1 = *(element+nel);
element++;
value0 += pi[iRow0]*element0;
value1 += pi[iRow1]*element1;
}
row += nel;
element += nel;
if (fabs(value0)>zeroTolerance) {
array[numberNonZero]=value0;
index[numberNonZero++]=*column;
}
column++;
if (fabs(value1)>zeroTolerance) {
array[numberNonZero]=value1;
index[numberNonZero++]=*column;
}
column++;
}
break;
}
#else
for (;numberPrice;numberPrice--) {
double value = 0.0;
int nel2=nel;
for (; nel2;nel2--) {
int iRow = *row++;
value += pi[iRow]*(*element++);
}
if (fabs(value)>zeroTolerance) {
array[numberNonZero]=value;
index[numberNonZero++]=*column;
}
column++;
}
#endif
}
output->setNumElements(numberNonZero);
}
// Updates two arrays for steepest
void
ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
const double * pi, CoinIndexedVector * output,
const double * piWeight,
double referenceIn, double devex,
// Array for exact devex to say what is in reference framework
unsigned int * reference,
double * weights, double scaleFactor)
{
int numberNonZero=0;
int * index = output->getIndices();
double * array = output->denseVector();
double zeroTolerance = model->factorization()->zeroTolerance();
double value = 0.0;
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=1.0;
int numberOdd = block_->startIndices_;
int iColumn;
CoinBigIndex end = start_[0];
for (iColumn=0;iColumn<numberOdd;iColumn++) {
CoinBigIndex start = end;
CoinBigIndex j;
int jColumn = column_[iColumn];
end = start_[iColumn+1];
value = 0.0;
if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
for (j=start; j<end;j++) {
int iRow = row_[j];
value -= pi[iRow]*element_[j];
}
if (fabs(value)>zeroTolerance) {
// and do other array
double modification = 0.0;
for (j=start; j<end;j++) {
int iRow = row_[j];
modification += piWeight[iRow]*element_[j];
}
double thisWeight = weights[jColumn];
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(jColumn))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[jColumn] = thisWeight;
if (!killDjs) {
array[numberNonZero]=value;
index[numberNonZero++]=jColumn;
}
}
}
}
for (int iBlock=0;iBlock<numberBlocks_;iBlock++) {
// B) Can sort so just do nonbasic (and nonfixed)
// C) Can do two at a time (if so put odd one into start_)
// D) can use switch
blockStruct * block = block_ + iBlock;
//int numberPrice = block->numberInBlock_;
int numberPrice = block->numberPrice_;
int nel = block->numberElements_;
int * row = row_ + block->startElements_;
double * element = element_ + block->startElements_;
int * column = column_ + block->startIndices_;
for (;numberPrice;numberPrice--) {
double value = 0.0;
int nel2=nel;
for (; nel2;nel2--) {
int iRow = *row++;
value -= pi[iRow]*(*element++);
}
if (fabs(value)>zeroTolerance) {
int jColumn = *column;
// back to beginning
row -= nel;
element -= nel;
// and do other array
double modification = 0.0;
nel2=nel;
for (; nel2;nel2--) {
int iRow = *row++;
modification += piWeight[iRow]*(*element++);
}
double thisWeight = weights[jColumn];
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(jColumn))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[jColumn] = thisWeight;
if (!killDjs) {
array[numberNonZero]=value;
index[numberNonZero++]=jColumn;
}
}
column++;
}
}
output->setNumElements(numberNonZero);
output->setPackedMode(true);
}
#ifdef CLP_ALL_ONE_FILE
#undef reference
#endif
syntax highlighted by Code2HTML, v. 0.9.1