// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include "CoinPragma.hpp"
#include <iostream>
#include "CoinIndexedVector.hpp"
#include "CoinHelperFunctions.hpp"
#include "ClpMatrixBase.hpp"
#include "ClpSimplex.hpp"
//#############################################################################
// Constructors / Destructor / Assignment
//#############################################################################
//-------------------------------------------------------------------
// Default Constructor
//-------------------------------------------------------------------
ClpMatrixBase::ClpMatrixBase () :
rhsOffset_(NULL),
startFraction_(0.0),
endFraction_(1.0),
savedBestDj_(0.0),
originalWanted_(0),
currentWanted_(0),
savedBestSequence_(-1),
type_(-1),
lastRefresh_(-1),
refreshFrequency_(0),
minimumObjectsScan_(-1),
minimumGoodReducedCosts_(-1),
trueSequenceIn_(-1),
trueSequenceOut_(-1),
skipDualCheck_(false)
{
}
//-------------------------------------------------------------------
// Copy constructor
//-------------------------------------------------------------------
ClpMatrixBase::ClpMatrixBase (const ClpMatrixBase & rhs) :
type_(rhs.type_),
skipDualCheck_(rhs.skipDualCheck_)
{
startFraction_ = rhs.startFraction_;
endFraction_ = rhs.endFraction_;
savedBestDj_ = rhs.savedBestDj_;
originalWanted_ = rhs.originalWanted_;
currentWanted_ = rhs.currentWanted_;
savedBestSequence_ = rhs.savedBestSequence_;
lastRefresh_ = rhs.lastRefresh_;
refreshFrequency_ = rhs.refreshFrequency_;
minimumObjectsScan_ = rhs.minimumObjectsScan_;
minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
trueSequenceIn_ = rhs.trueSequenceIn_;
trueSequenceOut_ = rhs.trueSequenceOut_;
skipDualCheck_ = rhs.skipDualCheck_;
int numberRows = rhs.getNumRows();
if (rhs.rhsOffset_&&numberRows) {
rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
} else {
rhsOffset_=NULL;
}
}
//-------------------------------------------------------------------
// Destructor
//-------------------------------------------------------------------
ClpMatrixBase::~ClpMatrixBase ()
{
delete [] rhsOffset_;
}
//----------------------------------------------------------------
// Assignment operator
//-------------------------------------------------------------------
ClpMatrixBase &
ClpMatrixBase::operator=(const ClpMatrixBase& rhs)
{
if (this != &rhs) {
type_ = rhs.type_;
delete [] rhsOffset_;
int numberRows = rhs.getNumRows();
if (rhs.rhsOffset_&&numberRows) {
rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_,numberRows);
} else {
rhsOffset_=NULL;
}
startFraction_ = rhs.startFraction_;
endFraction_ = rhs.endFraction_;
savedBestDj_ = rhs.savedBestDj_;
originalWanted_ = rhs.originalWanted_;
currentWanted_ = rhs.currentWanted_;
savedBestSequence_ = rhs.savedBestSequence_;
lastRefresh_ = rhs.lastRefresh_;
refreshFrequency_ = rhs.refreshFrequency_;
minimumObjectsScan_ = rhs.minimumObjectsScan_;
minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
trueSequenceIn_ = rhs.trueSequenceIn_;
trueSequenceOut_ = rhs.trueSequenceOut_;
skipDualCheck_ = rhs.skipDualCheck_;
}
return *this;
}
// And for scaling - default aborts for when scaling not supported
void
ClpMatrixBase::times(double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale) const
{
if (rowScale) {
std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
abort();
} else {
times(scalar,x,y);
}
}
// And for scaling - default aborts for when scaling not supported
void
ClpMatrixBase::transposeTimes(double scalar,
const double * x, double * y,
const double * rowScale,
const double * columnScale,
double * spare) const
{
if (rowScale) {
std::cerr<<"Scaling not supported - ClpMatrixBase"<<std::endl;
abort();
} else {
transposeTimes(scalar,x,y);
}
}
/* Subset clone (without gaps). Duplicates are allowed
and order is as given.
Derived classes need not provide this as it may not always make
sense */
ClpMatrixBase *
ClpMatrixBase::subsetClone (
int numberRows, const int * whichRows,
int numberColumns, const int * whichColumns) const
{
std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
abort();
return NULL;
}
/* Given positive integer weights for each row fills in sum of weights
for each column (and slack).
Returns weights vector
Default returns vector of ones
*/
CoinBigIndex *
ClpMatrixBase::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
{
int number = model->numberRows()+model->numberColumns();
CoinBigIndex * weights = new CoinBigIndex[number];
int i;
for (i=0;i<number;i++)
weights[i]=1;
return weights;
}
#ifndef CLP_NO_VECTOR
// Append Columns
void
ClpMatrixBase::appendCols(int number, const CoinPackedVectorBase * const * columns)
{
std::cerr<<"appendCols not supported - ClpMatrixBase"<<std::endl;
abort();
}
// Append Rows
void
ClpMatrixBase::appendRows(int number, const CoinPackedVectorBase * const * rows)
{
std::cerr<<"appendRows not supported - ClpMatrixBase"<<std::endl;
abort();
}
#endif
/* Returns largest and smallest elements of both signs.
Largest refers to largest absolute value.
*/
void
ClpMatrixBase::rangeOfElements(double & smallestNegative, double & largestNegative,
double & smallestPositive, double & largestPositive)
{
smallestNegative=0.0;
largestNegative=0.0;
smallestPositive=0.0;
largestPositive=0.0;
}
/* The length of a major-dimension vector. */
int
ClpMatrixBase::getVectorLength(int index) const
{
return getVectorLengths()[index];
}
// Says whether it can do partial pricing
bool
ClpMatrixBase::canDoPartialPricing() const
{
return false; // default is no
}
/* Return <code>x *A</code> in <code>z</code> but
just for number indices in y.
Default cheats with fake CoinIndexedVector and
then calls subsetTransposeTimes */
void
ClpMatrixBase::listTransposeTimes(const ClpSimplex * model,
double * x,
int * y,
int number,
double * z) const
{
CoinIndexedVector pi;
CoinIndexedVector list;
CoinIndexedVector output;
int * saveIndices = list.getIndices();
list.setNumElements(number);
list.setIndexVector(y);
double * savePi = pi.denseVector();
pi.setDenseVector(x);
double * saveOutput = output.denseVector();
output.setDenseVector(z);
output.setPacked();
subsetTransposeTimes(model,&pi,&list,&output);
// restore settings
list.setIndexVector(saveIndices);
pi.setDenseVector(savePi);
output.setDenseVector(saveOutput);
}
// Partial pricing
void
ClpMatrixBase::partialPricing(ClpSimplex * model, double start, double end,
int & bestSequence, int & numberWanted)
{
std::cerr<<"partialPricing not supported - ClpMatrixBase"<<std::endl;
abort();
}
/* expands an updated column to allow for extra rows which the main
solver does not know about and returns number added. If the arrays are NULL
then returns number of extra entries needed.
This will normally be a no-op - it is in for GUB!
*/
int
ClpMatrixBase::extendUpdated(ClpSimplex * model,CoinIndexedVector * update,int mode)
{
return 0;
}
/*
utility primal function for dealing with dynamic constraints
mode=n see ClpGubMatrix.hpp for definition
Remember to update here when settled down
*/
void
ClpMatrixBase::primalExpanded(ClpSimplex * model,int mode)
{
}
/*
utility dual function for dealing with dynamic constraints
mode=n see ClpGubMatrix.hpp for definition
Remember to update here when settled down
*/
void
ClpMatrixBase::dualExpanded(ClpSimplex * model,
CoinIndexedVector * array,
double * other,int mode)
{
}
/*
general utility function for dealing with dynamic constraints
mode=n see ClpGubMatrix.hpp for definition
Remember to update here when settled down
*/
int
ClpMatrixBase::generalExpanded(ClpSimplex * model,int mode, int &number)
{
int returnCode=0;
switch (mode) {
// Fill in pivotVariable but not for key variables
case 0:
{
int i;
int numberBasic=number;
int numberColumns = model->numberColumns();
// Use different array so can build from true pivotVariable_
//int * pivotVariable = model->pivotVariable();
int * pivotVariable = model->rowArray(0)->getIndices();
for (i=0;i<numberColumns;i++) {
if (model->getColumnStatus(i) == ClpSimplex::basic)
pivotVariable[numberBasic++]=i;
}
number = numberBasic;
}
break;
// Do initial extra rows + maximum basic
case 2:
{
number = model->numberRows();
}
break;
// To see if can dual or primal
case 4:
{
returnCode= 3;
}
break;
default:
break;
}
return returnCode;
}
// Sets up an effective RHS
void
ClpMatrixBase::useEffectiveRhs(ClpSimplex * model)
{
std::cerr<<"useEffectiveRhs not supported - ClpMatrixBase"<<std::endl;
abort();
}
/* Returns effective RHS if it is being used. This is used for long problems
or big gub or anywhere where going through full columns is
expensive. This may re-compute */
double *
ClpMatrixBase::rhsOffset(ClpSimplex * model,bool forceRefresh,bool check)
{
if (rhsOffset_) {
#ifdef CLP_DEBUG
if (check) {
// no need - but check anyway
// zero out basic
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
double * solution = new double [numberColumns];
double * rhs = new double[numberRows];
const double * solutionSlack = model->solutionRegion(0);
CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
int iRow;
for (iRow=0;iRow<numberRows;iRow++) {
if (model->getRowStatus(iRow)!=ClpSimplex::basic)
rhs[iRow]=solutionSlack[iRow];
else
rhs[iRow]=0.0;
}
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
solution[iColumn]=0.0;
}
times(-1.0,solution,rhs);
delete [] solution;
for (iRow=0;iRow<numberRows;iRow++) {
if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
}
}
#endif
if (forceRefresh||(refreshFrequency_&&model->numberIterations()>=
lastRefresh_+refreshFrequency_)) {
// zero out basic
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
double * solution = new double [numberColumns];
const double * solutionSlack = model->solutionRegion(0);
CoinMemcpyN(model->solutionRegion(),numberColumns,solution);
for (int iRow=0;iRow<numberRows;iRow++) {
if (model->getRowStatus(iRow)!=ClpSimplex::basic)
rhsOffset_[iRow]=solutionSlack[iRow];
else
rhsOffset_[iRow]=0.0;
}
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
if (model->getColumnStatus(iColumn)==ClpSimplex::basic)
solution[iColumn]=0.0;
}
times(-1.0,solution,rhsOffset_);
delete [] solution;
lastRefresh_ = model->numberIterations();
}
}
return rhsOffset_;
}
/*
update information for a pivot (and effective rhs)
*/
int
ClpMatrixBase::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
{
if (rhsOffset_) {
// update effective rhs
int sequenceIn = model->sequenceIn();
int sequenceOut = model->sequenceOut();
double * solution = model->solutionRegion();
int numberColumns = model->numberColumns();
if (sequenceIn==sequenceOut) {
if (sequenceIn<numberColumns)
add(model,rhsOffset_,sequenceIn,oldInValue-solution[sequenceIn]);
} else {
if (sequenceIn<numberColumns)
add(model,rhsOffset_,sequenceIn,oldInValue);
if (sequenceOut<numberColumns)
add(model,rhsOffset_,sequenceOut,-solution[sequenceOut]);
}
}
return 0;
}
int
ClpMatrixBase::hiddenRows() const
{
return 0;
}
/* Creates a variable. This is called after partial pricing and may modify matrix.
May update bestSequence.
*/
void
ClpMatrixBase::createVariable(ClpSimplex * model, int & bestSequence)
{
}
// Returns reduced cost of a variable
double
ClpMatrixBase::reducedCost(ClpSimplex * model,int sequence) const
{
int numberRows = model->numberRows();
int numberColumns = model->numberColumns();
if (sequence<numberRows+numberColumns)
return model->djRegion()[sequence];
else
return savedBestDj_;
}
/* Just for debug if odd type matrix.
Returns number and sum of primal infeasibilities.
*/
int
ClpMatrixBase::checkFeasible(ClpSimplex * model, double & sum) const
{
int numberRows = model->numberRows();
double * rhs = new double[numberRows];
int numberColumns = model->numberColumns();
int iRow;
CoinZeroN(rhs,numberRows);
times(1.0,model->solutionRegion(),rhs,model->rowScale(),model->columnScale());
int iColumn;
int logLevel = model->messageHandler()->logLevel();
int numberInfeasible=0;
const double * rowLower = model->lowerRegion(0);
const double * rowUpper = model->upperRegion(0);
const double * solution;
solution = model->solutionRegion(0);
double tolerance = model->primalTolerance()*1.01;
sum=0.0;
for (iRow=0;iRow<numberRows;iRow++) {
double value=rhs[iRow];
double value2= solution[iRow];
if (logLevel>3) {
if (fabs(value-value2)>1.0e-8)
printf("Row %d stored %g, computed %g\n",iRow,value2,value);
}
if (value<rowLower[iRow]-tolerance||
value>rowUpper[iRow]+tolerance) {
numberInfeasible++;
sum += CoinMax(rowLower[iRow]-value,value-rowUpper[iRow]);
}
if (value2>rowLower[iRow]+tolerance&&
value2<rowUpper[iRow]-tolerance&&
model->getRowStatus(iRow)!=ClpSimplex::basic) {
assert (model->getRowStatus(iRow)==ClpSimplex::superBasic);
}
}
const double * columnLower = model->lowerRegion(1);
const double * columnUpper = model->upperRegion(1);
solution = model->solutionRegion(1);
for (iColumn=0;iColumn<numberColumns;iColumn++) {
double value= solution[iColumn];
if (value<columnLower[iColumn]-tolerance||
value>columnUpper[iColumn]+tolerance) {
numberInfeasible++;
sum += CoinMax(columnLower[iColumn]-value,value-columnUpper[iColumn]);
}
if (value>columnLower[iColumn]+tolerance&&
value<columnUpper[iColumn]-tolerance&&
model->getColumnStatus(iColumn)!=ClpSimplex::basic) {
assert (model->getColumnStatus(iColumn)==ClpSimplex::superBasic);
}
}
delete [] rhs;
return numberInfeasible;
}
// These have to match ClpPrimalColumnSteepest version
#define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
// Updates second array for steepest and does devex weights (need not be coded)
void
ClpMatrixBase::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)
{
// get subset which have nonzero tableau elements
subsetTransposeTimes(model,pi2,dj1,dj2);
bool killDjs = (scaleFactor==0.0);
if (!scaleFactor)
scaleFactor=1.0;
// columns
int number = dj1->getNumElements();
const int * index = dj1->getIndices();
double * updateBy = dj1->denseVector();
double * updateBy2 = dj2->denseVector();
for (int j=0;j<number;j++) {
double thisWeight;
double pivot;
double pivotSquared;
int iSequence = index[j];
double value2 = updateBy[j];
if (killDjs)
updateBy[j]=0.0;
double modification=updateBy2[j];
updateBy2[j]=0.0;
ClpSimplex::Status status = model->getStatus(iSequence);
if (status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) {
thisWeight = weights[iSequence];
pivot = value2*scaleFactor;
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(iSequence))
thisWeight += 1.0;
thisWeight = CoinMax(thisWeight,DEVEX_TRY_NORM);
}
}
weights[iSequence] = thisWeight;
}
}
dj2->setNumElements(0);
}
// Correct sequence in and out to give true value
void
ClpMatrixBase::correctSequence(const ClpSimplex * model,int & sequenceIn, int & sequenceOut)
{
}
// Really scale matrix
void
ClpMatrixBase::reallyScale(const double * rowScale, const double * columnScale)
{
std::cerr<<"reallyScale not supported - ClpMatrixBase"<<std::endl;
abort();
}
// Updates two arrays for steepest
void
ClpMatrixBase::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)
{
std::cerr<<"transposeTimes2 not supported - ClpMatrixBase"<<std::endl;
abort();
}
/* 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
ClpMatrixBase::setDimensions(int numrows, int numcols){
// If odd matrix assume user knows what they are doing
}
/* 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
ClpMatrixBase::appendMatrix(int number, int type,
const CoinBigIndex * starts, const int * index,
const double * element, int numberOther)
{
std::cerr<<"appendMatrix not supported - ClpMatrixBase"<<std::endl;
abort();
return -1;
}
/* Modify one element of packed matrix. An element may be added.
This works for either ordering If the new element is zero it will be
deleted unless keepZero true */
void
ClpMatrixBase::modifyCoefficient(int row, int column, double newElement,
bool keepZero)
{
std::cerr<<"modifyCoefficient not supported - ClpMatrixBase"<<std::endl;
abort();
}
syntax highlighted by Code2HTML, v. 0.9.1