limp-cbc-0.3.2.0: cbits/coin/ClpInterior.cpp
/* $Id: ClpInterior.cpp 1941 2013-04-10 16:52:27Z stefan $ */
// Copyright (C) 2002, International Business Machines
// Corporation and others. All Rights Reserved.
// This code is licensed under the terms of the Eclipse Public License (EPL).
#include "CoinPragma.hpp"
#include <math.h>
#include "CoinHelperFunctions.hpp"
#include "ClpInterior.hpp"
#include "ClpMatrixBase.hpp"
#include "ClpLsqr.hpp"
#include "ClpPdcoBase.hpp"
#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(CoinWorkDouble));
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(CoinWorkDouble));
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(CoinWorkDouble));
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_);
lsqrObject_ = rhs.lsqrObject_ != NULL ? new ClpLsqr(*rhs.lsqrObject_) : NULL;
pdcoStuff_ = rhs.pdcoStuff_ != NULL ? rhs.pdcoStuff_->clone() : NULL;
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;
delete lsqrObject_;
lsqrObject_ = NULL;
//delete pdcoStuff_; // FIXME
pdcoStuff_ = 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 [] 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 CoinWorkDouble[nTotal];
CoinMemcpyN(columnActivity_, numberColumns_, solution_);
CoinMemcpyN(rowActivity_, numberRows_, solution_ + numberColumns_);
delete [] cost_;
cost_ = new CoinWorkDouble[nTotal];
int i;
CoinWorkDouble 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(CoinWorkDouble));
// do scaling if needed
if (scalingFlag_ > 0 && !rowScale_) {
if (matrix_->scale(this))
scalingFlag_ = -scalingFlag_; // not scaled after all
}
delete [] lower_;
delete [] upper_;
lower_ = new CoinWorkDouble[nTotal];
upper_ = new CoinWorkDouble[nTotal];
rowLowerWork_ = lower_ + numberColumns_;
columnLowerWork_ = lower_;
rowUpperWork_ = upper_ + numberColumns_;
columnUpperWork_ = upper_;
CoinMemcpyN(rowLower_, numberRows_, rowLowerWork_);
CoinMemcpyN(rowUpper_, numberRows_, rowUpperWork_);
CoinMemcpyN(columnLower_, numberColumns_, columnLowerWork_);
CoinMemcpyN(columnUpper_, numberColumns_, columnUpperWork_);
// 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++) {
CoinWorkDouble 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++) {
CoinWorkDouble 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 CoinWorkDouble [numberRows_];
assert (!rhsFixRegion_);
rhsFixRegion_ = new CoinWorkDouble [numberRows_];
assert (!deltaY_);
deltaY_ = new CoinWorkDouble [numberRows_];
CoinZeroN(deltaY_, numberRows_);
assert (!upperSlack_);
upperSlack_ = new CoinWorkDouble [nTotal];
assert (!lowerSlack_);
lowerSlack_ = new CoinWorkDouble [nTotal];
assert (!diagonal_);
diagonal_ = new CoinWorkDouble [nTotal];
assert (!deltaX_);
deltaX_ = new CoinWorkDouble [nTotal];
CoinZeroN(deltaX_, nTotal);
assert (!deltaZ_);
deltaZ_ = new CoinWorkDouble [nTotal];
CoinZeroN(deltaZ_, nTotal);
assert (!deltaW_);
deltaW_ = new CoinWorkDouble [nTotal];
CoinZeroN(deltaW_, nTotal);
assert (!deltaSU_);
deltaSU_ = new CoinWorkDouble [nTotal];
CoinZeroN(deltaSU_, nTotal);
assert (!deltaSL_);
deltaSL_ = new CoinWorkDouble [nTotal];
CoinZeroN(deltaSL_, nTotal);
assert (!primalR_);
assert (!dualR_);
// create arrays if we are doing KKT
if (cholesky_->type() >= 20) {
primalR_ = new CoinWorkDouble [nTotal];
CoinZeroN(primalR_, nTotal);
dualR_ = new CoinWorkDouble [numberRows_];
CoinZeroN(dualR_, numberRows_);
}
assert (!rhsB_);
rhsB_ = new CoinWorkDouble [numberRows_];
CoinZeroN(rhsB_, numberRows_);
assert (!rhsU_);
rhsU_ = new CoinWorkDouble [nTotal];
CoinZeroN(rhsU_, nTotal);
assert (!rhsL_);
rhsL_ = new CoinWorkDouble [nTotal];
CoinZeroN(rhsL_, nTotal);
assert (!rhsZ_);
rhsZ_ = new CoinWorkDouble [nTotal];
CoinZeroN(rhsZ_, nTotal);
assert (!rhsW_);
rhsW_ = new CoinWorkDouble [nTotal];
CoinZeroN(rhsW_, nTotal);
assert (!rhsC_);
rhsC_ = new CoinWorkDouble [nTotal];
CoinZeroN(rhsC_, nTotal);
assert (!workArray_);
workArray_ = new CoinWorkDouble [nTotal];
CoinZeroN(workArray_, nTotal);
assert (!zVec_);
zVec_ = new CoinWorkDouble [nTotal];
CoinZeroN(zVec_, nTotal);
assert (!wVec_);
wVec_ = new CoinWorkDouble [nTotal];
CoinZeroN(wVec_, nTotal);
assert (!dj_);
dj_ = new CoinWorkDouble [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) {
CoinWorkDouble 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_) {
CoinWorkDouble scaleR = 1.0 / rhsScale_;
for (i = 0; i < numberColumns_; i++) {
CoinWorkDouble scaleFactor = columnScale_[i];
CoinWorkDouble valueScaled = columnActivity_[i];
columnActivity_[i] = valueScaled * scaleFactor * scaleR;
CoinWorkDouble valueScaledDual = reducedCost_[i];
reducedCost_[i] = valueScaledDual / scaleFactor;
}
for (i = 0; i < numberRows_; i++) {
CoinWorkDouble scaleFactor = rowScale_[i];
CoinWorkDouble valueScaled = rowActivity_[i];
rowActivity_[i] = (valueScaled * scaleR) / scaleFactor;
CoinWorkDouble valueScaledDual = dual_[i];
dual_[i] = valueScaledDual * scaleFactor;
}
} else if (rhsScale_ != 1.0) {
CoinWorkDouble scaleR = 1.0 / rhsScale_;
for (i = 0; i < numberColumns_; i++) {
CoinWorkDouble valueScaled = columnActivity_[i];
columnActivity_[i] = valueScaled * scaleR;
}
for (i = 0; i < numberRows_; i++) {
CoinWorkDouble 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 ;
CoinWorkDouble largestBound, smallestBound, minimumGap;
CoinWorkDouble 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
CoinWorkDouble fixTolerance = 1.1 * primalTolerance();
for (i = numberColumns_; i < numberColumns_ + numberRows_; i++) {
CoinWorkDouble value;
value = CoinAbs(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 = CoinAbs(lower_[i]);
if (value > largestBound)
largestBound = value;
if (value < smallestBound)
smallestBound = value;
}
if (upper_[i] < 1.0e100 && upper_[i]) {
value = CoinAbs(upper_[i]);
if (value > largestBound)
largestBound = value;
if (value < smallestBound)
smallestBound = value;
}
}
if (largestBound)
handler_->message(CLP_RIMSTATISTICS3, messages_)
<< static_cast<double>(smallestBound)
<< static_cast<double>(largestBound)
<< static_cast<double>(minimumGap)
<< CoinMessageEol;
minimumGap = 1.0e100;
smallestBound = 1.0e100;
largestBound = 0.0;
for (i = 0; i < numberColumns_; i++) {
CoinWorkDouble value;
value = CoinAbs(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 = CoinAbs(lower_[i]);
if (value > largestBound)
largestBound = value;
if (value < smallestBound)
smallestBound = value;
}
if (upper_[i] < 1.0e100 && upper_[i]) {
value = CoinAbs(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_)
<< static_cast<double>(smallestObj)
<< static_cast<double>(largestObj)
<< CoinMessageEol;
if (largestBound)
handler_->message(CLP_RIMSTATISTICS2, messages_)
<< static_cast<double>(smallestBound)
<< static_cast<double>(largestBound)
<< static_cast<double>(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;
}
#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);
}
#include "ClpPredictorCorrector.hpp"
// Primal-Dual Predictor-Corrector barrier
int
ClpInterior::primalDual()
{
return (static_cast<ClpPredictorCorrector *> (this))->solve();
}
void
ClpInterior::checkSolution()
{
int iRow, iColumn;
CoinWorkDouble * reducedCost = reinterpret_cast<CoinWorkDouble *>(reducedCost_);
CoinWorkDouble * dual = reinterpret_cast<CoinWorkDouble *>(dual_);
CoinMemcpyN(cost_, numberColumns_, reducedCost);
matrix_->transposeTimes(-1.0, dual, reducedCost);
// Now modify reduced costs for quadratic
CoinWorkDouble quadraticOffset = quadraticDjs(reducedCost,
solution_, scaleFactor_);
objectiveValue_ = 0.0;
// now look at solution
sumPrimalInfeasibilities_ = 0.0;
sumDualInfeasibilities_ = 0.0;
CoinWorkDouble dualTolerance = 10.0 * dblParam_[ClpDualTolerance];
CoinWorkDouble primalTolerance = dblParam_[ClpPrimalTolerance];
CoinWorkDouble primalTolerance2 = 10.0 * dblParam_[ClpPrimalTolerance];
worstComplementarity_ = 0.0;
complementarityGap_ = 0.0;
// Done scaled - use permanent regions for output
// but internal for bounds
const CoinWorkDouble * lower = lower_ + numberColumns_;
const CoinWorkDouble * upper = upper_ + numberColumns_;
for (iRow = 0; iRow < numberRows_; iRow++) {
CoinWorkDouble infeasibility = 0.0;
CoinWorkDouble distanceUp = CoinMin(upper[iRow] -
rowActivity_[iRow], static_cast<CoinWorkDouble>(1.0e10));
CoinWorkDouble distanceDown = CoinMin(rowActivity_[iRow] -
lower[iRow], static_cast<CoinWorkDouble>(1.0e10));
if (distanceUp > primalTolerance2) {
CoinWorkDouble 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) {
CoinWorkDouble 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++) {
CoinWorkDouble infeasibility = 0.0;
objectiveValue_ += cost_[iColumn] * columnActivity_[iColumn];
CoinWorkDouble distanceUp = CoinMin(upper[iColumn] -
columnActivity_[iColumn], static_cast<CoinWorkDouble>(1.0e10));
CoinWorkDouble distanceDown = CoinMin(columnActivity_[iColumn] -
lower[iColumn], static_cast<CoinWorkDouble>(1.0e10));
if (distanceUp > primalTolerance2) {
CoinWorkDouble 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) {
CoinWorkDouble 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;
}
}
#if COIN_LONG_WORK
// ok as packs down
CoinMemcpyN(reducedCost, numberColumns_, reducedCost_);
#endif
// 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
CoinWorkDouble * columnChange = new CoinWorkDouble[numberColumns_];
CoinWorkDouble * rowChange = new CoinWorkDouble[numberRows_];
CoinZeroN(columnChange, numberColumns_);
CoinZeroN(rowChange, numberRows_);
matrix_->times(1.0, columnChange, rowChange);
int i;
CoinWorkDouble 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]) {
CoinWorkDouble change = columnLower_[i] - columnActivity_[i];
if (CoinAbs(change) < tolerance) {
if (reallyFix)
columnUpper_[i] = columnLower_[i];
columnChange[i] = change;
columnActivity_[i] = columnLower_[i];
}
} else {
CoinWorkDouble change = columnUpper_[i] - columnActivity_[i];
if (CoinAbs(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
CoinWorkDouble newSum = 0.0;
for (i = 0; i < numberRows_; i++) {
CoinWorkDouble 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]) {
CoinWorkDouble change = rowLower_[i] - rowActivity_[i];
if (CoinAbs(change) < tolerance) {
if (reallyFix)
rowUpper_[i] = rowLower_[i];
rowActivity_[i] = rowLower_[i];
}
} else {
CoinWorkDouble change = rowLower_[i] - rowActivity_[i];
if (CoinAbs(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 */
CoinWorkDouble
ClpInterior::quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
CoinWorkDouble scaleFactor)
{
CoinWorkDouble 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++) {
CoinWorkDouble value = 0.0;
for (CoinBigIndex j = columnQuadraticStart[iColumn];
j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
int jColumn = columnQuadratic[j];
CoinWorkDouble valueJ = solution[jColumn];
CoinWorkDouble elementValue = quadraticElement[j];
//value += valueI*valueJ*elementValue;
value += valueJ * elementValue;
quadraticOffset += solution[iColumn] * valueJ * elementValue;
}
djRegion[iColumn] += scaleFactor * value;
}
}
return quadraticOffset;
}