// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
#include "CoinPragma.hpp"
#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "ClpInterior.hpp"
#include "ClpMatrixBase.hpp"
#ifdef PDCO
#include "ClpLsqr.hpp"
#include "ClpPdcoBase.hpp"
#endif
#include "CoinDenseVector.hpp"
#include "ClpMessage.hpp"
#include "ClpHelperFunctions.hpp"
#include "ClpCholeskyDense.hpp"
#include "ClpLinearObjective.hpp"
#include "ClpQuadraticObjective.hpp"
#include <cfloat>
#include <string>
#include <stdio.h>
#include <iostream>
//#############################################################################
ClpInterior::ClpInterior () :
ClpModel(),
largestPrimalError_(0.0),
largestDualError_(0.0),
sumDualInfeasibilities_(0.0),
sumPrimalInfeasibilities_(0.0),
worstComplementarity_(0.0),
xsize_(0.0),
zsize_(0.0),
lower_(NULL),
rowLowerWork_(NULL),
columnLowerWork_(NULL),
upper_(NULL),
rowUpperWork_(NULL),
columnUpperWork_(NULL),
cost_(NULL),
rhs_(NULL),
x_(NULL),
y_(NULL),
dj_(NULL),
lsqrObject_(NULL),
pdcoStuff_(NULL),
mu_(0.0),
objectiveNorm_(1.0e-12),
rhsNorm_(1.0e-12),
solutionNorm_(1.0e-12),
dualObjective_(0.0),
primalObjective_(0.0),
diagonalNorm_(1.0e-12),
stepLength_(0.995),
linearPerturbation_(1.0e-12),
diagonalPerturbation_(1.0e-15),
gamma_(0.0),
delta_(0),
targetGap_(1.0e-12),
projectionTolerance_(1.0e-7),
maximumRHSError_(0.0),
maximumBoundInfeasibility_(0.0),
maximumDualError_(0.0),
diagonalScaleFactor_(0.0),
scaleFactor_(1.0),
actualPrimalStep_(0.0),
actualDualStep_(0.0),
smallestInfeasibility_(0.0),
complementarityGap_(0.0),
baseObjectiveNorm_(0.0),
worstDirectionAccuracy_(0.0),
maximumRHSChange_(0.0),
errorRegion_(NULL),
rhsFixRegion_(NULL),
upperSlack_(NULL),
lowerSlack_(NULL),
diagonal_(NULL),
solution_(NULL),
workArray_(NULL),
deltaX_(NULL),
deltaY_(NULL),
deltaZ_(NULL),
deltaW_(NULL),
deltaSU_(NULL),
deltaSL_(NULL),
primalR_(NULL),
dualR_(NULL),
rhsB_(NULL),
rhsU_(NULL),
rhsL_(NULL),
rhsZ_(NULL),
rhsW_(NULL),
rhsC_(NULL),
zVec_(NULL),
wVec_(NULL),
cholesky_(NULL),
numberComplementarityPairs_(0),
numberComplementarityItems_(0),
maximumBarrierIterations_(200),
gonePrimalFeasible_(false),
goneDualFeasible_(false),
algorithm_(-1)
{
memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
solveType_=3; // say interior based life form
cholesky_ = new ClpCholeskyDense(); // put in placeholder
}
// Subproblem constructor
ClpInterior::ClpInterior ( const ClpModel * rhs,
int numberRows, const int * whichRow,
int numberColumns, const int * whichColumn,
bool dropNames, bool dropIntegers)
: ClpModel(rhs, numberRows, whichRow,
numberColumns,whichColumn,dropNames,dropIntegers),
largestPrimalError_(0.0),
largestDualError_(0.0),
sumDualInfeasibilities_(0.0),
sumPrimalInfeasibilities_(0.0),
worstComplementarity_(0.0),
xsize_(0.0),
zsize_(0.0),
lower_(NULL),
rowLowerWork_(NULL),
columnLowerWork_(NULL),
upper_(NULL),
rowUpperWork_(NULL),
columnUpperWork_(NULL),
cost_(NULL),
rhs_(NULL),
x_(NULL),
y_(NULL),
dj_(NULL),
lsqrObject_(NULL),
pdcoStuff_(NULL),
mu_(0.0),
objectiveNorm_(1.0e-12),
rhsNorm_(1.0e-12),
solutionNorm_(1.0e-12),
dualObjective_(0.0),
primalObjective_(0.0),
diagonalNorm_(1.0e-12),
stepLength_(0.99995),
linearPerturbation_(1.0e-12),
diagonalPerturbation_(1.0e-15),
gamma_(0.0),
delta_(0),
targetGap_(1.0e-12),
projectionTolerance_(1.0e-7),
maximumRHSError_(0.0),
maximumBoundInfeasibility_(0.0),
maximumDualError_(0.0),
diagonalScaleFactor_(0.0),
scaleFactor_(0.0),
actualPrimalStep_(0.0),
actualDualStep_(0.0),
smallestInfeasibility_(0.0),
complementarityGap_(0.0),
baseObjectiveNorm_(0.0),
worstDirectionAccuracy_(0.0),
maximumRHSChange_(0.0),
errorRegion_(NULL),
rhsFixRegion_(NULL),
upperSlack_(NULL),
lowerSlack_(NULL),
diagonal_(NULL),
solution_(NULL),
workArray_(NULL),
deltaX_(NULL),
deltaY_(NULL),
deltaZ_(NULL),
deltaW_(NULL),
deltaSU_(NULL),
deltaSL_(NULL),
primalR_(NULL),
dualR_(NULL),
rhsB_(NULL),
rhsU_(NULL),
rhsL_(NULL),
rhsZ_(NULL),
rhsW_(NULL),
rhsC_(NULL),
zVec_(NULL),
wVec_(NULL),
cholesky_(NULL),
numberComplementarityPairs_(0),
numberComplementarityItems_(0),
maximumBarrierIterations_(200),
gonePrimalFeasible_(false),
goneDualFeasible_(false),
algorithm_(-1)
{
memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
solveType_=3; // say interior based life form
cholesky_= new ClpCholeskyDense();
}
//-----------------------------------------------------------------------------
ClpInterior::~ClpInterior ()
{
gutsOfDelete();
}
//#############################################################################
/*
This does housekeeping
*/
int
ClpInterior::housekeeping()
{
numberIterations_++;
return 0;
}
// Copy constructor.
ClpInterior::ClpInterior(const ClpInterior &rhs) :
ClpModel(rhs),
largestPrimalError_(0.0),
largestDualError_(0.0),
sumDualInfeasibilities_(0.0),
sumPrimalInfeasibilities_(0.0),
worstComplementarity_(0.0),
xsize_(0.0),
zsize_(0.0),
lower_(NULL),
rowLowerWork_(NULL),
columnLowerWork_(NULL),
upper_(NULL),
rowUpperWork_(NULL),
columnUpperWork_(NULL),
cost_(NULL),
rhs_(NULL),
x_(NULL),
y_(NULL),
dj_(NULL),
lsqrObject_(NULL),
pdcoStuff_(NULL),
errorRegion_(NULL),
rhsFixRegion_(NULL),
upperSlack_(NULL),
lowerSlack_(NULL),
diagonal_(NULL),
solution_(NULL),
workArray_(NULL),
deltaX_(NULL),
deltaY_(NULL),
deltaZ_(NULL),
deltaW_(NULL),
deltaSU_(NULL),
deltaSL_(NULL),
primalR_(NULL),
dualR_(NULL),
rhsB_(NULL),
rhsU_(NULL),
rhsL_(NULL),
rhsZ_(NULL),
rhsW_(NULL),
rhsC_(NULL),
zVec_(NULL),
wVec_(NULL),
cholesky_(NULL)
{
gutsOfDelete();
gutsOfCopy(rhs);
solveType_=3; // say interior based life form
}
// Copy constructor from model
ClpInterior::ClpInterior(const ClpModel &rhs) :
ClpModel(rhs),
largestPrimalError_(0.0),
largestDualError_(0.0),
sumDualInfeasibilities_(0.0),
sumPrimalInfeasibilities_(0.0),
worstComplementarity_(0.0),
xsize_(0.0),
zsize_(0.0),
lower_(NULL),
rowLowerWork_(NULL),
columnLowerWork_(NULL),
upper_(NULL),
rowUpperWork_(NULL),
columnUpperWork_(NULL),
cost_(NULL),
rhs_(NULL),
x_(NULL),
y_(NULL),
dj_(NULL),
lsqrObject_(NULL),
pdcoStuff_(NULL),
mu_(0.0),
objectiveNorm_(1.0e-12),
rhsNorm_(1.0e-12),
solutionNorm_(1.0e-12),
dualObjective_(0.0),
primalObjective_(0.0),
diagonalNorm_(1.0e-12),
stepLength_(0.99995),
linearPerturbation_(1.0e-12),
diagonalPerturbation_(1.0e-15),
gamma_(0.0),
delta_(0),
targetGap_(1.0e-12),
projectionTolerance_(1.0e-7),
maximumRHSError_(0.0),
maximumBoundInfeasibility_(0.0),
maximumDualError_(0.0),
diagonalScaleFactor_(0.0),
scaleFactor_(0.0),
actualPrimalStep_(0.0),
actualDualStep_(0.0),
smallestInfeasibility_(0.0),
complementarityGap_(0.0),
baseObjectiveNorm_(0.0),
worstDirectionAccuracy_(0.0),
maximumRHSChange_(0.0),
errorRegion_(NULL),
rhsFixRegion_(NULL),
upperSlack_(NULL),
lowerSlack_(NULL),
diagonal_(NULL),
solution_(NULL),
workArray_(NULL),
deltaX_(NULL),
deltaY_(NULL),
deltaZ_(NULL),
deltaW_(NULL),
deltaSU_(NULL),
deltaSL_(NULL),
primalR_(NULL),
dualR_(NULL),
rhsB_(NULL),
rhsU_(NULL),
rhsL_(NULL),
rhsZ_(NULL),
rhsW_(NULL),
rhsC_(NULL),
zVec_(NULL),
wVec_(NULL),
cholesky_(NULL),
numberComplementarityPairs_(0),
numberComplementarityItems_(0),
maximumBarrierIterations_(200),
gonePrimalFeasible_(false),
goneDualFeasible_(false),
algorithm_(-1)
{
memset(historyInfeasibility_,0,LENGTH_HISTORY*sizeof(double));
solveType_=3; // say interior based life form
cholesky_= new ClpCholeskyDense();
}
// Assignment operator. This copies the data
ClpInterior &
ClpInterior::operator=(const ClpInterior & rhs)
{
if (this != &rhs) {
gutsOfDelete();
ClpModel::operator=(rhs);
gutsOfCopy(rhs);
}
return *this;
}
void
ClpInterior::gutsOfCopy(const ClpInterior & rhs)
{
lower_ = ClpCopyOfArray(rhs.lower_,numberColumns_+numberRows_);
rowLowerWork_ = lower_+numberColumns_;
columnLowerWork_ = lower_;
upper_ = ClpCopyOfArray(rhs.upper_,numberColumns_+numberRows_);
rowUpperWork_ = upper_+numberColumns_;
columnUpperWork_ = upper_;
//cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
cost_ = ClpCopyOfArray(rhs.cost_,numberColumns_);
rhs_ = ClpCopyOfArray(rhs.rhs_,numberRows_);
x_ = ClpCopyOfArray(rhs.x_,numberColumns_);
y_ = ClpCopyOfArray(rhs.y_,numberRows_);
dj_ = ClpCopyOfArray(rhs.dj_,numberColumns_+numberRows_);
#ifdef PDCO
lsqrObject_= new ClpLsqr(*rhs.lsqrObject_);
pdcoStuff_ = rhs.pdcoStuff_->clone();
#endif
largestPrimalError_ = rhs.largestPrimalError_;
largestDualError_ = rhs.largestDualError_;
sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
worstComplementarity_ = rhs.worstComplementarity_;
xsize_ = rhs.xsize_;
zsize_ = rhs.zsize_;
solveType_=rhs.solveType_;
mu_ = rhs.mu_;
objectiveNorm_ = rhs.objectiveNorm_;
rhsNorm_ = rhs.rhsNorm_;
solutionNorm_ = rhs.solutionNorm_;
dualObjective_ = rhs.dualObjective_;
primalObjective_ = rhs.primalObjective_;
diagonalNorm_ = rhs.diagonalNorm_;
stepLength_ = rhs.stepLength_;
linearPerturbation_ = rhs.linearPerturbation_;
diagonalPerturbation_ = rhs.diagonalPerturbation_;
gamma_=rhs.gamma_;
delta_=rhs.delta_;
targetGap_ = rhs.targetGap_;
projectionTolerance_ = rhs.projectionTolerance_;
maximumRHSError_ = rhs.maximumRHSError_;
maximumBoundInfeasibility_ = rhs.maximumBoundInfeasibility_;
maximumDualError_ = rhs.maximumDualError_;
diagonalScaleFactor_ = rhs.diagonalScaleFactor_;
scaleFactor_ = rhs.scaleFactor_;
actualPrimalStep_ = rhs.actualPrimalStep_;
actualDualStep_ = rhs.actualDualStep_;
smallestInfeasibility_ = rhs.smallestInfeasibility_;
complementarityGap_ = rhs.complementarityGap_;
baseObjectiveNorm_ = rhs.baseObjectiveNorm_;
worstDirectionAccuracy_ = rhs.worstDirectionAccuracy_;
maximumRHSChange_ = rhs.maximumRHSChange_;
errorRegion_ = ClpCopyOfArray(rhs.errorRegion_,numberRows_);
rhsFixRegion_ = ClpCopyOfArray(rhs.rhsFixRegion_,numberRows_);
deltaY_ = ClpCopyOfArray(rhs.deltaY_,numberRows_);
upperSlack_ = ClpCopyOfArray(rhs.upperSlack_,numberRows_+numberColumns_);
lowerSlack_ = ClpCopyOfArray(rhs.lowerSlack_,numberRows_+numberColumns_);
diagonal_ = ClpCopyOfArray(rhs.diagonal_,numberRows_+numberColumns_);
deltaX_ = ClpCopyOfArray(rhs.deltaX_,numberRows_+numberColumns_);
deltaZ_ = ClpCopyOfArray(rhs.deltaZ_,numberRows_+numberColumns_);
deltaW_ = ClpCopyOfArray(rhs.deltaW_,numberRows_+numberColumns_);
deltaSU_ = ClpCopyOfArray(rhs.deltaSU_,numberRows_+numberColumns_);
deltaSL_ = ClpCopyOfArray(rhs.deltaSL_,numberRows_+numberColumns_);
primalR_ = ClpCopyOfArray(rhs.primalR_,numberRows_+numberColumns_);
dualR_ = ClpCopyOfArray(rhs.dualR_,numberRows_+numberColumns_);
rhsB_ = ClpCopyOfArray(rhs.rhsB_,numberRows_);
rhsU_ = ClpCopyOfArray(rhs.rhsU_,numberRows_+numberColumns_);
rhsL_ = ClpCopyOfArray(rhs.rhsL_,numberRows_+numberColumns_);
rhsZ_ = ClpCopyOfArray(rhs.rhsZ_,numberRows_+numberColumns_);
rhsW_ = ClpCopyOfArray(rhs.rhsW_,numberRows_+numberColumns_);
rhsC_ = ClpCopyOfArray(rhs.rhsC_,numberRows_+numberColumns_);
solution_ = ClpCopyOfArray(rhs.solution_,numberRows_+numberColumns_);
workArray_ = ClpCopyOfArray(rhs.workArray_,numberRows_+numberColumns_);
zVec_ = ClpCopyOfArray(rhs.zVec_,numberRows_+numberColumns_);
wVec_ = ClpCopyOfArray(rhs.wVec_,numberRows_+numberColumns_);
cholesky_ = rhs.cholesky_->clone();
numberComplementarityPairs_ = rhs.numberComplementarityPairs_;
numberComplementarityItems_ = rhs.numberComplementarityItems_;
maximumBarrierIterations_ = rhs.maximumBarrierIterations_;
gonePrimalFeasible_ = rhs.gonePrimalFeasible_;
goneDualFeasible_ = rhs.goneDualFeasible_;
algorithm_ = rhs.algorithm_;
}
void
ClpInterior::gutsOfDelete()
{
delete [] lower_;
lower_=NULL;
rowLowerWork_=NULL;
columnLowerWork_=NULL;
delete [] upper_;
upper_=NULL;
rowUpperWork_=NULL;
columnUpperWork_=NULL;
delete [] cost_;
cost_=NULL;
delete [] rhs_;
rhs_ = NULL;
delete [] x_;
x_ = NULL;
delete [] y_;
y_ = NULL;
delete [] dj_;
dj_ = NULL;
#ifdef PDCO
delete lsqrObject_;
lsqrObject_ = NULL;
//delete pdcoStuff_;
pdcoStuff_=NULL;
#endif
delete [] errorRegion_;
errorRegion_ = NULL;
delete [] rhsFixRegion_;
rhsFixRegion_ = NULL;
delete [] deltaY_;
deltaY_ = NULL;
delete [] upperSlack_;
upperSlack_ = NULL;
delete [] lowerSlack_;
lowerSlack_ = NULL;
delete [] diagonal_;
diagonal_ = NULL;
delete [] deltaX_;
deltaX_ = NULL;
delete [] deltaZ_;
deltaZ_ = NULL;
delete [] deltaW_;
deltaW_ = NULL;
delete [] deltaSU_;
deltaSU_ = NULL;
delete [] deltaSL_;
deltaSL_ = NULL;
delete [] primalR_;
primalR_=NULL;
delete [] dualR_;
dualR_=NULL;
delete [] rhsB_;
rhsB_ = NULL;
delete [] rhsU_;
rhsU_ = NULL;
delete [] rhsL_;
rhsL_ = NULL;
delete [] rhsZ_;
rhsZ_ = NULL;
delete [] rhsW_;
rhsW_ = NULL;
delete [] rhsC_;
rhsC_ = NULL;
delete [] solution_;
solution_ = NULL;
delete [] workArray_;
workArray_ = NULL;
delete [] zVec_;
zVec_ = NULL;
delete [] wVec_;
wVec_ = NULL;
delete cholesky_;
}
bool
ClpInterior::createWorkingData()
{
bool goodMatrix=true;
//check matrix
if (!matrix_->allElementsInRange(this,1.0e-12,1.0e20)) {
problemStatus_=4;
goodMatrix= false;
}
int nTotal = numberRows_+numberColumns_;
delete [] solution_;
solution_ = new double[nTotal];
memcpy(solution_,columnActivity_,
numberColumns_*sizeof(double));
memcpy(solution_+numberColumns_,rowActivity_,
numberRows_*sizeof(double));
delete [] cost_;
cost_ = new double[nTotal];
int i;
double direction = optimizationDirection_*objectiveScale_;
// direction is actually scale out not scale in
if (direction)
direction = 1.0/direction;
const double * obj = objective();
for (i=0;i<numberColumns_;i++)
cost_[i] = direction*obj[i];
memset(cost_+numberColumns_,0,numberRows_*sizeof(double));
// do scaling if needed
if (scalingFlag_>0&&!rowScale_) {
if (matrix_->scale(this))
scalingFlag_=-scalingFlag_; // not scaled after all
}
delete [] lower_;
delete [] upper_;
lower_ = new double[nTotal];
upper_ = new double[nTotal];
rowLowerWork_ = lower_+numberColumns_;
columnLowerWork_ = lower_;
rowUpperWork_ = upper_+numberColumns_;
columnUpperWork_ = upper_;
memcpy(rowLowerWork_,rowLower_,numberRows_*sizeof(double));
memcpy(rowUpperWork_,rowUpper_,numberRows_*sizeof(double));
memcpy(columnLowerWork_,columnLower_,numberColumns_*sizeof(double));
memcpy(columnUpperWork_,columnUpper_,numberColumns_*sizeof(double));
// clean up any mismatches on infinity
for (i=0;i<numberColumns_;i++) {
if (columnLowerWork_[i]<-1.0e30)
columnLowerWork_[i] = -COIN_DBL_MAX;
if (columnUpperWork_[i]>1.0e30)
columnUpperWork_[i] = COIN_DBL_MAX;
}
// clean up any mismatches on infinity
for (i=0;i<numberRows_;i++) {
if (rowLowerWork_[i]<-1.0e30)
rowLowerWork_[i] = -COIN_DBL_MAX;
if (rowUpperWork_[i]>1.0e30)
rowUpperWork_[i] = COIN_DBL_MAX;
}
// check rim of problem okay
if (!sanityCheck())
goodMatrix=false;
if(rowScale_) {
for (i=0;i<numberColumns_;i++) {
double multiplier = rhsScale_/columnScale_[i];
cost_[i] *= columnScale_[i];
if (columnLowerWork_[i]>-1.0e50)
columnLowerWork_[i] *= multiplier;
if (columnUpperWork_[i]<1.0e50)
columnUpperWork_[i] *= multiplier;
}
for (i=0;i<numberRows_;i++) {
double multiplier = rhsScale_*rowScale_[i];
if (rowLowerWork_[i]>-1.0e50)
rowLowerWork_[i] *= multiplier;
if (rowUpperWork_[i]<1.0e50)
rowUpperWork_[i] *= multiplier;
}
} else if (rhsScale_!=1.0) {
for (i=0;i<numberColumns_+numberRows_;i++) {
if (lower_[i]>-1.0e50)
lower_[i] *= rhsScale_;
if (upper_[i]<1.0e50)
upper_[i] *= rhsScale_;
}
}
assert (!errorRegion_);
errorRegion_ = new double [numberRows_];
assert (!rhsFixRegion_);
rhsFixRegion_ = new double [numberRows_];
assert (!deltaY_);
deltaY_ = new double [numberRows_];
CoinZeroN(deltaY_,numberRows_);
assert (!upperSlack_);
upperSlack_ = new double [nTotal];
assert (!lowerSlack_);
lowerSlack_ = new double [nTotal];
assert (!diagonal_);
diagonal_ = new double [nTotal];
assert (!deltaX_);
deltaX_ = new double [nTotal];
CoinZeroN(deltaX_,nTotal);
assert (!deltaZ_);
deltaZ_ = new double [nTotal];
CoinZeroN(deltaZ_,nTotal);
assert (!deltaW_);
deltaW_ = new double [nTotal];
CoinZeroN(deltaW_,nTotal);
assert (!deltaSU_);
deltaSU_ = new double [nTotal];
CoinZeroN(deltaSU_,nTotal);
assert (!deltaSL_);
deltaSL_ = new double [nTotal];
CoinZeroN(deltaSL_,nTotal);
assert (!primalR_);
assert (!dualR_);
// create arrays if we are doing KKT
if (cholesky_->type()>=20) {
primalR_ = new double [nTotal];
CoinZeroN(primalR_,nTotal);
dualR_ = new double [numberRows_];
CoinZeroN(dualR_,numberRows_);
}
assert (!rhsB_);
rhsB_ = new double [numberRows_];
CoinZeroN(rhsB_,numberRows_);
assert (!rhsU_);
rhsU_ = new double [nTotal];
CoinZeroN(rhsU_,nTotal);
assert (!rhsL_);
rhsL_ = new double [nTotal];
CoinZeroN(rhsL_,nTotal);
assert (!rhsZ_);
rhsZ_ = new double [nTotal];
CoinZeroN(rhsZ_,nTotal);
assert (!rhsW_);
rhsW_ = new double [nTotal];
CoinZeroN(rhsW_,nTotal);
assert (!rhsC_);
rhsC_ = new double [nTotal];
CoinZeroN(rhsC_,nTotal);
assert (!workArray_);
workArray_ = new double [nTotal];
CoinZeroN(workArray_,nTotal);
assert (!zVec_);
zVec_ = new double [nTotal];
CoinZeroN(zVec_,nTotal);
assert (!wVec_);
wVec_ = new double [nTotal];
CoinZeroN(wVec_,nTotal);
assert (!dj_);
dj_ = new double [nTotal];
if (!status_)
status_ = new unsigned char [numberRows_+numberColumns_];
memset(status_,0,numberRows_+numberColumns_);
return goodMatrix;
}
void
ClpInterior::deleteWorkingData()
{
int i;
if (optimizationDirection_!=1.0||objectiveScale_!=1.0) {
double scaleC = optimizationDirection_/objectiveScale_;
// and modify all dual signs
for (i=0;i<numberColumns_;i++)
reducedCost_[i] = scaleC*dj_[i];
for (i=0;i<numberRows_;i++)
dual_[i] *= scaleC;
}
if (rowScale_) {
double scaleR = 1.0/rhsScale_;
for (i=0;i<numberColumns_;i++) {
double scaleFactor = columnScale_[i];
double valueScaled = columnActivity_[i];
columnActivity_[i] = valueScaled*scaleFactor*scaleR;
double valueScaledDual = reducedCost_[i];
reducedCost_[i] = valueScaledDual/scaleFactor;
}
for (i=0;i<numberRows_;i++) {
double scaleFactor = rowScale_[i];
double valueScaled = rowActivity_[i];
rowActivity_[i] = (valueScaled*scaleR)/scaleFactor;
double valueScaledDual = dual_[i];
dual_[i] = valueScaledDual*scaleFactor;
}
} else if (rhsScale_!=1.0) {
double scaleR = 1.0/rhsScale_;
for (i=0;i<numberColumns_;i++) {
double valueScaled = columnActivity_[i];
columnActivity_[i] = valueScaled*scaleR;
}
for (i=0;i<numberRows_;i++) {
double valueScaled = rowActivity_[i];
rowActivity_[i] = valueScaled*scaleR;
}
}
delete [] cost_;
cost_ = NULL;
delete [] solution_;
solution_ = NULL;
delete [] lower_;
lower_ = NULL;
delete [] upper_;
upper_ = NULL;
delete [] errorRegion_;
errorRegion_ = NULL;
delete [] rhsFixRegion_;
rhsFixRegion_ = NULL;
delete [] deltaY_;
deltaY_ = NULL;
delete [] upperSlack_;
upperSlack_ = NULL;
delete [] lowerSlack_;
lowerSlack_ = NULL;
delete [] diagonal_;
diagonal_ = NULL;
delete [] deltaX_;
deltaX_ = NULL;
delete [] workArray_;
workArray_ = NULL;
delete [] zVec_;
zVec_ = NULL;
delete [] wVec_;
wVec_ = NULL;
delete [] dj_;
dj_ = NULL;
}
// Sanity check on input data - returns true if okay
bool
ClpInterior::sanityCheck()
{
// bad if empty
if (!numberColumns_||((!numberRows_||!matrix_->getNumElements())&&objective_->type()<2)) {
problemStatus_=emptyProblem();
return false;
}
int numberBad ;
double largestBound, smallestBound, minimumGap;
double smallestObj, largestObj;
int firstBad;
int modifiedBounds=0;
int i;
numberBad=0;
firstBad=-1;
minimumGap=1.0e100;
smallestBound=1.0e100;
largestBound=0.0;
smallestObj=1.0e100;
largestObj=0.0;
// If bounds are too close - fix
double fixTolerance = 1.1*primalTolerance();
for (i=numberColumns_;i<numberColumns_+numberRows_;i++) {
double value;
value = fabs(cost_[i]);
if (value>1.0e50) {
numberBad++;
if (firstBad<0)
firstBad=i;
} else if (value) {
if (value>largestObj)
largestObj=value;
if (value<smallestObj)
smallestObj=value;
}
value=upper_[i]-lower_[i];
if (value<-primalTolerance()) {
numberBad++;
if (firstBad<0)
firstBad=i;
} else if (value<=fixTolerance) {
if (value) {
// modify
upper_[i] = lower_[i];
modifiedBounds++;
}
} else {
if (value<minimumGap)
minimumGap=value;
}
if (lower_[i]>-1.0e100&&lower_[i]) {
value = fabs(lower_[i]);
if (value>largestBound)
largestBound=value;
if (value<smallestBound)
smallestBound=value;
}
if (upper_[i]<1.0e100&&upper_[i]) {
value = fabs(upper_[i]);
if (value>largestBound)
largestBound=value;
if (value<smallestBound)
smallestBound=value;
}
}
if (largestBound)
handler_->message(CLP_RIMSTATISTICS3,messages_)
<<smallestBound
<<largestBound
<<minimumGap
<<CoinMessageEol;
minimumGap=1.0e100;
smallestBound=1.0e100;
largestBound=0.0;
for (i=0;i<numberColumns_;i++) {
double value;
value = fabs(cost_[i]);
if (value>1.0e50) {
numberBad++;
if (firstBad<0)
firstBad=i;
} else if (value) {
if (value>largestObj)
largestObj=value;
if (value<smallestObj)
smallestObj=value;
}
value=upper_[i]-lower_[i];
if (value<-primalTolerance()) {
numberBad++;
if (firstBad<0)
firstBad=i;
} else if (value<=fixTolerance) {
if (value) {
// modify
upper_[i] = lower_[i];
modifiedBounds++;
}
} else {
if (value<minimumGap)
minimumGap=value;
}
if (lower_[i]>-1.0e100&&lower_[i]) {
value = fabs(lower_[i]);
if (value>largestBound)
largestBound=value;
if (value<smallestBound)
smallestBound=value;
}
if (upper_[i]<1.0e100&&upper_[i]) {
value = fabs(upper_[i]);
if (value>largestBound)
largestBound=value;
if (value<smallestBound)
smallestBound=value;
}
}
char rowcol[]={'R','C'};
if (numberBad) {
handler_->message(CLP_BAD_BOUNDS,messages_)
<<numberBad
<<rowcol[isColumn(firstBad)]<<sequenceWithin(firstBad)
<<CoinMessageEol;
problemStatus_=4;
return false;
}
if (modifiedBounds)
handler_->message(CLP_MODIFIEDBOUNDS,messages_)
<<modifiedBounds
<<CoinMessageEol;
handler_->message(CLP_RIMSTATISTICS1,messages_)
<<smallestObj
<<largestObj
<<CoinMessageEol; if (largestBound)
handler_->message(CLP_RIMSTATISTICS2,messages_)
<<smallestBound
<<largestBound
<<minimumGap
<<CoinMessageEol;
return true;
}
/* Loads a problem (the constraints on the
rows are given by lower and upper bounds). If a pointer is 0 then the
following values are the default:
<ul>
<li> <code>colub</code>: all columns have upper bound infinity
<li> <code>collb</code>: all columns have lower bound 0
<li> <code>rowub</code>: all rows have upper bound infinity
<li> <code>rowlb</code>: all rows have lower bound -infinity
<li> <code>obj</code>: all variables have 0 objective coefficient
</ul>
*/
void
ClpInterior::loadProblem ( const ClpMatrixBase& matrix,
const double* collb, const double* colub,
const double* obj,
const double* rowlb, const double* rowub,
const double * rowObjective)
{
ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
rowObjective);
}
void
ClpInterior::loadProblem ( const CoinPackedMatrix& matrix,
const double* collb, const double* colub,
const double* obj,
const double* rowlb, const double* rowub,
const double * rowObjective)
{
ClpModel::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
rowObjective);
}
/* Just like the other loadProblem() method except that the matrix is
given in a standard column major ordered format (without gaps). */
void
ClpInterior::loadProblem ( const int numcols, const int numrows,
const CoinBigIndex* start, const int* index,
const double* value,
const double* collb, const double* colub,
const double* obj,
const double* rowlb, const double* rowub,
const double * rowObjective)
{
ClpModel::loadProblem(numcols, numrows, start, index, value,
collb, colub, obj, rowlb, rowub,
rowObjective);
}
void
ClpInterior::loadProblem ( const int numcols, const int numrows,
const CoinBigIndex* start, const int* index,
const double* value,const int * length,
const double* collb, const double* colub,
const double* obj,
const double* rowlb, const double* rowub,
const double * rowObjective)
{
ClpModel::loadProblem(numcols, numrows, start, index, value, length,
collb, colub, obj, rowlb, rowub,
rowObjective);
}
// Read an mps file from the given filename
int
ClpInterior::readMps(const char *filename,
bool keepNames,
bool ignoreErrors)
{
int status = ClpModel::readMps(filename,keepNames,ignoreErrors);
return status;
}
#ifdef PDCO
#include "ClpPdco.hpp"
/* Pdco algorithm - see ClpPdco.hpp for method */
int
ClpInterior::pdco()
{
return ((ClpPdco *) this)->pdco();
}
// ** Temporary version
int
ClpInterior::pdco( ClpPdcoBase * stuff, Options &options, Info &info, Outfo &outfo)
{
return ((ClpPdco *) this)->pdco(stuff,options,info,outfo);
}
#endif
#include "ClpPredictorCorrector.hpp"
// Primal-Dual Predictor-Corrector barrier
int
ClpInterior::primalDual()
{
return ((ClpPredictorCorrector *) this)->solve();
}
void
ClpInterior::checkSolution()
{
int iRow,iColumn;
CoinMemcpyN(cost_,numberColumns_,reducedCost_);
matrix_->transposeTimes(-1.0,dual_,reducedCost_);
// Now modify reduced costs for quadratic
double quadraticOffset=quadraticDjs(reducedCost_,
solution_,scaleFactor_);
objectiveValue_ = 0.0;
// now look at solution
sumPrimalInfeasibilities_=0.0;
sumDualInfeasibilities_=0.0;
double dualTolerance = 10.0*dblParam_[ClpDualTolerance];
double primalTolerance = dblParam_[ClpPrimalTolerance];
double primalTolerance2 = 10.0*dblParam_[ClpPrimalTolerance];
worstComplementarity_=0.0;
complementarityGap_=0.0;
// Done scaled - use permanent regions for output
// but internal for bounds
const double * lower = lower_+numberColumns_;
const double * upper = upper_+numberColumns_;
for (iRow=0;iRow<numberRows_;iRow++) {
double infeasibility=0.0;
double distanceUp = CoinMin(upper[iRow]-
rowActivity_[iRow],1.0e10);
double distanceDown = CoinMin(rowActivity_[iRow] -
lower[iRow],1.0e10);
if (distanceUp>primalTolerance2) {
double value = dual_[iRow];
// should not be negative
if (value<-dualTolerance) {
sumDualInfeasibilities_ += -dualTolerance-value;
value = - value*distanceUp;
if (value>worstComplementarity_)
worstComplementarity_=value;
complementarityGap_ += value;
}
}
if (distanceDown>primalTolerance2) {
double value = dual_[iRow];
// should not be positive
if (value>dualTolerance) {
sumDualInfeasibilities_ += value-dualTolerance;
value = value*distanceDown;
if (value>worstComplementarity_)
worstComplementarity_=value;
complementarityGap_ += value;
}
}
if (rowActivity_[iRow]>upper[iRow]) {
infeasibility=rowActivity_[iRow]-upper[iRow];
} else if (rowActivity_[iRow]<lower[iRow]) {
infeasibility=lower[iRow]-rowActivity_[iRow];
}
if (infeasibility>primalTolerance) {
sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
}
}
lower = lower_;
upper = upper_;
for (iColumn=0;iColumn<numberColumns_;iColumn++) {
double infeasibility=0.0;
objectiveValue_ += cost_[iColumn]*columnActivity_[iColumn];
double distanceUp = CoinMin(upper[iColumn]-
columnActivity_[iColumn],1.0e10);
double distanceDown = CoinMin(columnActivity_[iColumn] -
lower[iColumn],1.0e10);
if (distanceUp>primalTolerance2) {
double value = reducedCost_[iColumn];
// should not be negative
if (value<-dualTolerance) {
sumDualInfeasibilities_ += -dualTolerance-value;
value = - value*distanceUp;
if (value>worstComplementarity_)
worstComplementarity_=value;
complementarityGap_ += value;
}
}
if (distanceDown>primalTolerance2) {
double value = reducedCost_[iColumn];
// should not be positive
if (value>dualTolerance) {
sumDualInfeasibilities_ += value-dualTolerance;
value = value*distanceDown;
if (value>worstComplementarity_)
worstComplementarity_=value;
complementarityGap_ += value;
}
}
if (columnActivity_[iColumn]>upper[iColumn]) {
infeasibility=columnActivity_[iColumn]-upper[iColumn];
} else if (columnActivity_[iColumn]<lower[iColumn]) {
infeasibility=lower[iColumn]-columnActivity_[iColumn];
}
if (infeasibility>primalTolerance) {
sumPrimalInfeasibilities_ += infeasibility-primalTolerance;
}
}
// add in offset
objectiveValue_ += 0.5*quadraticOffset;
}
// Set cholesky (and delete present one)
void
ClpInterior::setCholesky(ClpCholeskyBase * cholesky)
{
delete cholesky_;
cholesky_= cholesky;
}
/* Borrow model. This is so we dont have to copy large amounts
of data around. It assumes a derived class wants to overwrite
an empty model with a real one - while it does an algorithm.
This is same as ClpModel one. */
void
ClpInterior::borrowModel(ClpModel & otherModel)
{
ClpModel::borrowModel(otherModel);
}
/* Return model - updates any scalars */
void
ClpInterior::returnModel(ClpModel & otherModel)
{
ClpModel::returnModel(otherModel);
}
// Return number fixed to see if worth presolving
int
ClpInterior::numberFixed() const
{
int i;
int nFixed=0;
for (i=0;i<numberColumns_;i++) {
if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) {
if (columnUpper_[i]>columnLower_[i]) {
if (fixedOrFree(i))
nFixed++;
}
}
}
for (i=0;i<numberRows_;i++) {
if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) {
if (rowUpper_[i]>rowLower_[i]) {
if (fixedOrFree(i+numberColumns_))
nFixed++;
}
}
}
return nFixed;
}
// fix variables interior says should be
void
ClpInterior::fixFixed(bool reallyFix)
{
// Arrays for change in columns and rhs
double * columnChange = new double[numberColumns_];
double * rowChange = new double[numberRows_];
CoinZeroN(columnChange,numberColumns_);
CoinZeroN(rowChange,numberRows_);
matrix_->times(1.0,columnChange,rowChange);
int i;
double tolerance = primalTolerance();
for (i=0;i<numberColumns_;i++) {
if (columnUpper_[i]<1.0e20||columnLower_[i]>-1.0e20) {
if (columnUpper_[i]>columnLower_[i]) {
if (fixedOrFree(i)) {
if (columnActivity_[i]-columnLower_[i]<columnUpper_[i]-columnActivity_[i]) {
double change = columnLower_[i]-columnActivity_[i];
if (fabs(change)<tolerance) {
if (reallyFix)
columnUpper_[i]=columnLower_[i];
columnChange[i] = change;
columnActivity_[i]=columnLower_[i];
}
} else {
double change = columnUpper_[i]-columnActivity_[i];
if (fabs(change)<tolerance) {
if (reallyFix)
columnLower_[i]=columnUpper_[i];
columnChange[i] = change;
columnActivity_[i]=columnUpper_[i];
}
}
}
}
}
}
CoinZeroN(rowChange,numberRows_);
matrix_->times(1.0,columnChange,rowChange);
// If makes mess of things then don't do
double newSum=0.0;
for (i=0;i<numberRows_;i++) {
double value=rowActivity_[i]+rowChange[i];
if (value>rowUpper_[i]+tolerance)
newSum += value-rowUpper_[i]-tolerance;
else if (value<rowLower_[i]-tolerance)
newSum -= value-rowLower_[i]+tolerance;
}
if (newSum>1.0e-5+1.5*sumPrimalInfeasibilities_) {
// put back and skip changes
for (i=0;i<numberColumns_;i++)
columnActivity_[i] -= columnChange[i];
} else {
CoinZeroN(rowActivity_,numberRows_);
matrix_->times(1.0,columnActivity_,rowActivity_);
if (reallyFix) {
for (i=0;i<numberRows_;i++) {
if (rowUpper_[i]<1.0e20||rowLower_[i]>-1.0e20) {
if (rowUpper_[i]>rowLower_[i]) {
if (fixedOrFree(i+numberColumns_)) {
if (rowActivity_[i]-rowLower_[i]<rowUpper_[i]-rowActivity_[i]) {
double change = rowLower_[i]-rowActivity_[i];
if (fabs(change)<tolerance) {
if (reallyFix)
rowUpper_[i]=rowLower_[i];
rowActivity_[i]=rowLower_[i];
}
} else {
double change = rowLower_[i]-rowActivity_[i];
if (fabs(change)<tolerance) {
if (reallyFix)
rowLower_[i]=rowUpper_[i];
rowActivity_[i]=rowUpper_[i];
}
}
}
}
}
}
}
}
delete [] rowChange;
delete [] columnChange;
}
/* Modifies djs to allow for quadratic.
returns quadratic offset */
double
ClpInterior::quadraticDjs(double * djRegion, const double * solution,
double scaleFactor)
{
double quadraticOffset=0.0;
#ifndef NO_RTTI
ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
#else
ClpQuadraticObjective * quadraticObj = NULL;
if (objective_->type()==2)
quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
#endif
if (quadraticObj) {
CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
const int * columnQuadratic = quadratic->getIndices();
const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
const int * columnQuadraticLength = quadratic->getVectorLengths();
double * quadraticElement = quadratic->getMutableElements();
int numberColumns = quadratic->getNumCols();
for (int iColumn=0;iColumn<numberColumns;iColumn++) {
double value = 0.0;
for (CoinBigIndex j=columnQuadraticStart[iColumn];
j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
int jColumn = columnQuadratic[j];
double valueJ = solution[jColumn];
double elementValue = quadraticElement[j];
//value += valueI*valueJ*elementValue;
value += valueJ*elementValue;
quadraticOffset += solution[iColumn]*valueJ*elementValue;
}
djRegion[iColumn] += scaleFactor*value;
}
}
return quadraticOffset;
}
syntax highlighted by Code2HTML, v. 0.9.1