/* $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;
}