/* $Id: ClpDynamicMatrix.cpp 1941 2013-04-10 16:52:27Z stefan $ */ // Copyright (C) 2004, International Business Machines // Corporation and others. All Rights Reserved. // This code is licensed under the terms of the Eclipse Public License (EPL). #include #include "CoinPragma.hpp" #include "CoinIndexedVector.hpp" #include "CoinHelperFunctions.hpp" #include "ClpSimplex.hpp" #include "ClpFactorization.hpp" #include "ClpQuadraticObjective.hpp" #include "ClpNonLinearCost.hpp" // at end to get min/max! #include "ClpDynamicMatrix.hpp" #include "ClpMessage.hpp" //#define CLP_DEBUG //#define CLP_DEBUG_PRINT //############################################################################# // Constructors / Destructor / Assignment //############################################################################# //------------------------------------------------------------------- // Default Constructor //------------------------------------------------------------------- ClpDynamicMatrix::ClpDynamicMatrix () : ClpPackedMatrix(), sumDualInfeasibilities_(0.0), sumPrimalInfeasibilities_(0.0), sumOfRelaxedDualInfeasibilities_(0.0), sumOfRelaxedPrimalInfeasibilities_(0.0), savedBestGubDual_(0.0), savedBestSet_(0), backToPivotRow_(NULL), keyVariable_(NULL), toIndex_(NULL), fromIndex_(NULL), numberSets_(0), numberActiveSets_(0), objectiveOffset_(0.0), lowerSet_(NULL), upperSet_(NULL), status_(NULL), model_(NULL), firstAvailable_(0), firstAvailableBefore_(0), firstDynamic_(0), lastDynamic_(0), numberStaticRows_(0), numberElements_(0), numberDualInfeasibilities_(0), numberPrimalInfeasibilities_(0), noCheck_(-1), infeasibilityWeight_(0.0), numberGubColumns_(0), maximumGubColumns_(0), maximumElements_(0), startSet_(NULL), next_(NULL), startColumn_(NULL), row_(NULL), element_(NULL), cost_(NULL), id_(NULL), dynamicStatus_(NULL), columnLower_(NULL), columnUpper_(NULL) { setType(15); } //------------------------------------------------------------------- // Copy constructor //------------------------------------------------------------------- ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs) : ClpPackedMatrix(rhs) { objectiveOffset_ = rhs.objectiveOffset_; numberSets_ = rhs.numberSets_; numberActiveSets_ = rhs.numberActiveSets_; firstAvailable_ = rhs.firstAvailable_; firstAvailableBefore_ = rhs.firstAvailableBefore_; firstDynamic_ = rhs.firstDynamic_; lastDynamic_ = rhs.lastDynamic_; numberStaticRows_ = rhs.numberStaticRows_; numberElements_ = rhs.numberElements_; backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_); keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_); lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); status_ = ClpCopyOfArray(rhs.status_, static_cast(2*numberSets_+4*sizeof(int))); model_ = rhs.model_; sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; savedBestGubDual_ = rhs.savedBestGubDual_; savedBestSet_ = rhs.savedBestSet_; noCheck_ = rhs.noCheck_; infeasibilityWeight_ = rhs.infeasibilityWeight_; // Now secondary data numberGubColumns_ = rhs.numberGubColumns_; maximumGubColumns_ = rhs.maximumGubColumns_; maximumElements_ = rhs.maximumElements_; startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1); next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); row_ = ClpCopyOfArray(rhs.row_, maximumElements_); element_ = ClpCopyOfArray(rhs.element_, maximumElements_); cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_); id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_); columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_); } /* This is the real constructor*/ ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * model, int numberSets, int numberGubColumns, const int * starts, const double * lower, const double * upper, const CoinBigIndex * startColumn, const int * row, const double * element, const double * cost, const double * columnLower, const double * columnUpper, const unsigned char * status, const unsigned char * dynamicStatus) : ClpPackedMatrix() { setType(15); objectiveOffset_ = model->objectiveOffset(); model_ = model; numberSets_ = numberSets; numberGubColumns_ = numberGubColumns; maximumGubColumns_ = numberGubColumns_; if (numberGubColumns_) maximumElements_ = startColumn[numberGubColumns_]; else maximumElements_ = 0; startSet_ = new int [numberSets_+1]; next_ = new int [maximumGubColumns_]; // fill in startSet and next int iSet; if (numberGubColumns_) { for (iSet = 0; iSet < numberSets_; iSet++) { int first = starts[iSet]; int last = starts[iSet+1] - 1; startSet_[iSet] = first; for (int i = first; i < last; i++) next_[i] = i + 1; next_[last] = -iSet - 1; } startSet_[numberSets_] = starts[numberSets_]; } int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); numberStaticRows_ = numberRows; savedBestGubDual_ = 0.0; savedBestSet_ = 0; // Number of columns needed int frequency = model->factorizationFrequency(); int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4; // But we may have two per row + one for incoming (make it two) numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall); // for small problems this could be too big //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_); int numberNeeded = numberGubInSmall + numberColumns; firstAvailable_ = numberColumns; firstAvailableBefore_ = firstAvailable_; firstDynamic_ = numberColumns; lastDynamic_ = numberNeeded; startColumn_ = ClpCopyOfArray(startColumn, numberGubColumns_ + 1); if (!numberGubColumns_) { if (!startColumn_) startColumn_ = new CoinBigIndex [1]; startColumn_[0] = 0; } CoinBigIndex numberElements = startColumn_[numberGubColumns_]; row_ = ClpCopyOfArray(row, numberElements); element_ = new double[numberElements]; CoinBigIndex i; for (i = 0; i < numberElements; i++) element_[i] = element[i]; cost_ = new double[numberGubColumns_]; for (i = 0; i < numberGubColumns_; i++) { cost_[i] = cost[i]; // I don't think I need sorted but ... CoinSort_2(row_ + startColumn_[i], row_ + startColumn_[i+1], element_ + startColumn_[i]); } if (columnLower) { columnLower_ = new double[numberGubColumns_]; for (i = 0; i < numberGubColumns_; i++) columnLower_[i] = columnLower[i]; } else { columnLower_ = NULL; } if (columnUpper) { columnUpper_ = new double[numberGubColumns_]; for (i = 0; i < numberGubColumns_; i++) columnUpper_[i] = columnUpper[i]; } else { columnUpper_ = NULL; } lowerSet_ = new double[numberSets_]; for (i = 0; i < numberSets_; i++) { if (lower[i] > -1.0e20) lowerSet_[i] = lower[i]; else lowerSet_[i] = -1.0e30; } upperSet_ = new double[numberSets_]; for (i = 0; i < numberSets_; i++) { if (upper[i] < 1.0e20) upperSet_[i] = upper[i]; else upperSet_[i] = 1.0e30; } id_ = new int[numberGubInSmall]; for (i = 0; i < numberGubInSmall; i++) id_[i] = -1; ClpPackedMatrix* originalMatrixA = dynamic_cast< ClpPackedMatrix*>(model->clpMatrix()); assert (originalMatrixA); CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix(); originalMatrixA->setMatrixNull(); // so can be deleted safely // guess how much space needed double guess = numberElements; guess /= static_cast (numberColumns); guess *= 2 * numberGubInSmall; numberElements_ = static_cast (guess); numberElements_ = CoinMin(numberElements_, numberElements) + originalMatrix->getNumElements(); matrix_ = originalMatrix; //delete originalMatrixA; flags_ &= ~1; // resize model (matrix stays same) // modify frequency if (frequency>=50) frequency = 50+(frequency-50)/2; int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1; model->resize(newRowSize, numberNeeded); for (i = numberRows; i < newRowSize; i++) model->setRowStatus(i, ClpSimplex::basic); if (columnUpper_) { // set all upper bounds so we have enough space double * columnUpper = model->columnUpper(); for(i = firstDynamic_; i < lastDynamic_; i++) columnUpper[i] = 1.0e10; } // resize matrix // extra 1 is so can keep number of elements handy originalMatrix->reserve(numberNeeded, numberElements_, true); originalMatrix->reserve(numberNeeded + 1, numberElements_, false); originalMatrix->getMutableVectorStarts()[numberColumns] = originalMatrix->getNumElements(); originalMatrix->setDimensions(newRowSize, -1); numberActiveColumns_ = firstDynamic_; // redo number of columns numberColumns = matrix_->getNumCols(); backToPivotRow_ = new int[numberNeeded]; keyVariable_ = new int[numberSets_]; if (status) { status_ = ClpCopyOfArray(status, static_cast(2*numberSets_+4*sizeof(int))); assert (dynamicStatus); dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_); } else { status_ = new unsigned char [2*numberSets_+4*sizeof(int)]; memset(status_, 0, numberSets_); int i; for (i = 0; i < numberSets_; i++) { // make slack key setStatus(i, ClpSimplex::basic); } dynamicStatus_ = new unsigned char [2*numberGubColumns_]; memset(dynamicStatus_, 0, numberGubColumns_); // for clarity for (i = 0; i < numberGubColumns_; i++) setDynamicStatus(i, atLowerBound); } toIndex_ = new int[numberSets_]; for (iSet = 0; iSet < numberSets_; iSet++) toIndex_[iSet] = -1; fromIndex_ = new int [newRowSize-numberStaticRows_+1]; numberActiveSets_ = 0; rhsOffset_ = NULL; if (numberGubColumns_) { if (!status) { gubCrash(); } else { initialProblem(); } } noCheck_ = -1; infeasibilityWeight_ = 0.0; } //------------------------------------------------------------------- // Destructor //------------------------------------------------------------------- ClpDynamicMatrix::~ClpDynamicMatrix () { delete [] backToPivotRow_; delete [] keyVariable_; delete [] toIndex_; delete [] fromIndex_; delete [] lowerSet_; delete [] upperSet_; delete [] status_; delete [] startSet_; delete [] next_; delete [] startColumn_; delete [] row_; delete [] element_; delete [] cost_; delete [] id_; delete [] dynamicStatus_; delete [] columnLower_; delete [] columnUpper_; } //---------------------------------------------------------------- // Assignment operator //------------------------------------------------------------------- ClpDynamicMatrix & ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs) { if (this != &rhs) { ClpPackedMatrix::operator=(rhs); delete [] backToPivotRow_; delete [] keyVariable_; delete [] toIndex_; delete [] fromIndex_; delete [] lowerSet_; delete [] upperSet_; delete [] status_; delete [] startSet_; delete [] next_; delete [] startColumn_; delete [] row_; delete [] element_; delete [] cost_; delete [] id_; delete [] dynamicStatus_; delete [] columnLower_; delete [] columnUpper_; objectiveOffset_ = rhs.objectiveOffset_; numberSets_ = rhs.numberSets_; numberActiveSets_ = rhs.numberActiveSets_; firstAvailable_ = rhs.firstAvailable_; firstAvailableBefore_ = rhs.firstAvailableBefore_; firstDynamic_ = rhs.firstDynamic_; lastDynamic_ = rhs.lastDynamic_; numberStaticRows_ = rhs.numberStaticRows_; numberElements_ = rhs.numberElements_; backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, lastDynamic_); keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + 1 - numberStaticRows_); lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_); upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_); status_ = ClpCopyOfArray(rhs.status_, static_cast(2*numberSets_+4*sizeof(int))); model_ = rhs.model_; sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; savedBestGubDual_ = rhs.savedBestGubDual_; savedBestSet_ = rhs.savedBestSet_; noCheck_ = rhs.noCheck_; infeasibilityWeight_ = rhs.infeasibilityWeight_; // Now secondary data numberGubColumns_ = rhs.numberGubColumns_; maximumGubColumns_ = rhs.maximumGubColumns_; maximumElements_ = rhs.maximumElements_; startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1); next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_); startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1); row_ = ClpCopyOfArray(rhs.row_, maximumElements_); element_ = ClpCopyOfArray(rhs.element_, maximumElements_); cost_ = ClpCopyOfArray(rhs.cost_, maximumGubColumns_); id_ = ClpCopyOfArray(rhs.id_, lastDynamic_ - firstDynamic_); columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_); columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_); dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_); } return *this; } //------------------------------------------------------------------- // Clone //------------------------------------------------------------------- ClpMatrixBase * ClpDynamicMatrix::clone() const { return new ClpDynamicMatrix(*this); } // Partial pricing void ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, int & bestSequence, int & numberWanted) { numberWanted = currentWanted_; assert(!model->rowScale()); if (numberSets_) { // Do packed part before gub // always??? //printf("normal packed price - start %d end %d (passed end %d, first %d)\n", ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); } else { // no gub ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); return; } if (numberWanted > 0) { // and do some proportion of full set int startG2 = static_cast (startFraction * numberSets_); int endG2 = static_cast (endFraction * numberSets_ + 0.1); endG2 = CoinMin(endG2, numberSets_); //printf("gub price - set start %d end %d\n", // startG2,endG2); double tolerance = model->currentDualTolerance(); double * reducedCost = model->djRegion(); const double * duals = model->dualRowSolution(); double bestDj; int numberRows = model->numberRows(); int slackOffset = lastDynamic_ + numberRows; int structuralOffset = slackOffset + numberSets_; // If nothing found yet can go all the way to end int endAll = endG2; if (bestSequence < 0 && !startG2) endAll = numberSets_; if (bestSequence >= 0) { if (bestSequence != savedBestSequence_) bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent else bestDj = savedBestDj_; } else { bestDj = tolerance; } int saveSequence = bestSequence; double djMod = 0.0; double bestDjMod = 0.0; //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(), // startG2,endG2,numberWanted); int bestSet = -1; #if 0 // make sure first available is clean (in case last iteration rejected) cost[firstAvailable_] = 0.0; length[firstAvailable_] = 0; model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0); model->setStatus(firstAvailable_, ClpSimplex::atLowerBound); { for (int i = firstAvailable_; i < lastDynamic_; i++) assert(!cost[i]); } #endif int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; int minNeg = minimumGoodReducedCosts_ < 0 ? 5 : minimumGoodReducedCosts_; for (int iSet = startG2; iSet < endAll; iSet++) { if (numberWanted + minNeg < originalWanted_ && iSet > startG2 + minSet) { // give up numberWanted = 0; break; } else if (iSet == endG2 && bestSequence >= 0) { break; } int gubRow = toIndex_[iSet]; if (gubRow >= 0) { djMod = duals[gubRow+numberStaticRows_]; // have I got sign right? } else { int iBasic = keyVariable_[iSet]; if (iBasic >= maximumGubColumns_) { djMod = 0.0; // set not in } else { // get dj without djMod = 0.0; for (CoinBigIndex j = startColumn_[iBasic]; j < startColumn_[iBasic+1]; j++) { int jRow = row_[j]; djMod -= duals[jRow] * element_[j]; } djMod += cost_[iBasic]; // See if gub slack possible - dj is djMod if (getStatus(iSet) == ClpSimplex::atLowerBound) { double value = -djMod; if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!flagged(iSet)) { bestDj = value; bestSequence = slackOffset + iSet; bestDjMod = djMod; bestSet = iSet; } else { // just to make sure we don't exit before got something numberWanted++; abort(); } } } } else if (getStatus(iSet) == ClpSimplex::atUpperBound) { double value = djMod; if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!flagged(iSet)) { bestDj = value; bestSequence = slackOffset + iSet; bestDjMod = djMod; bestSet = iSet; } else { // just to make sure we don't exit before got something numberWanted++; abort(); } } } } } } int iSequence = startSet_[iSet]; while (iSequence >= 0) { DynamicStatus status = getDynamicStatus(iSequence); if (status == atLowerBound || status == atUpperBound) { double value = cost_[iSequence] - djMod; for (CoinBigIndex j = startColumn_[iSequence]; j < startColumn_[iSequence+1]; j++) { int jRow = row_[j]; value -= duals[jRow] * element_[j]; } // change sign if at lower bound if (status == atLowerBound) value = -value; if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!flagged(iSequence)) { if (false/*status == atLowerBound &&keyValue(iSet)<1.0e-7*/) { // can't come in because // of ones at ub numberWanted++; } else { bestDj = value; bestSequence = structuralOffset + iSequence; bestDjMod = djMod; bestSet = iSet; } } else { // just to make sure we don't exit before got something numberWanted++; } } } } iSequence = next_[iSequence]; //onto next in set } if (numberWanted <= 0) { numberWanted = 0; break; } } if (bestSequence != saveSequence) { savedBestGubDual_ = bestDjMod; savedBestDj_ = bestDj; savedBestSequence_ = bestSequence; savedBestSet_ = bestSet; } // See if may be finished if (!startG2 && bestSequence < 0) infeasibilityWeight_ = model_->infeasibilityCost(); else if (bestSequence >= 0) infeasibilityWeight_ = -1.0; } currentWanted_ = numberWanted; } /* Returns effective RHS if it is being used. This is used for long problems or big gub or anywhere where going through full columns is expensive. This may re-compute */ double * ClpDynamicMatrix::rhsOffset(ClpSimplex * model, bool forceRefresh, bool #ifdef CLP_DEBUG check #endif ) { // forceRefresh=true;printf("take out forceRefresh\n"); if (!model_->numberIterations()) forceRefresh = true; //check=false; #ifdef CLP_DEBUG double * saveE = NULL; if (rhsOffset_ && check) { int numberRows = model->numberRows(); saveE = new double[numberRows]; } #endif if (rhsOffset_) { #ifdef CLP_DEBUG if (check) { // no need - but check anyway int numberRows = model->numberRows(); double * rhs = new double[numberRows]; int iRow; int iSet; CoinZeroN(rhs, numberRows); // do ones at bounds before gub const double * smallSolution = model->solutionRegion(); const double * element = matrix_->getElements(); const int * row = matrix_->getIndices(); const CoinBigIndex * startColumn = matrix_->getVectorStarts(); const int * length = matrix_->getVectorLengths(); int iColumn; double objectiveOffset = 0.0; for (iColumn = 0; iColumn < firstDynamic_; iColumn++) { if (model->getStatus(iColumn) != ClpSimplex::basic) { double value = smallSolution[iColumn]; for (CoinBigIndex j = startColumn[iColumn]; j < startColumn[iColumn] + length[iColumn]; j++) { int jRow = row[j]; rhs[jRow] -= value * element[j]; } } } if (columnLower_ || columnUpper_) { double * solution = new double [numberGubColumns_]; for (iSet = 0; iSet < numberSets_; iSet++) { int j = startSet_[iSet]; while (j >= 0) { double value = 0.0; if (getDynamicStatus(j) != inSmall) { if (getDynamicStatus(j) == atLowerBound) { if (columnLower_) value = columnLower_[j]; } else if (getDynamicStatus(j) == atUpperBound) { value = columnUpper_[j]; } else if (getDynamicStatus(j) == soloKey) { value = keyValue(iSet); } objectiveOffset += value * cost_[j]; } solution[j] = value; j = next_[j]; //onto next in set } } // ones in gub and in small problem for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) { if (model_->getStatus(iColumn) != ClpSimplex::basic) { int jFull = id_[iColumn-firstDynamic_]; solution[jFull] = smallSolution[iColumn]; } } for (iSet = 0; iSet < numberSets_; iSet++) { int kRow = toIndex_[iSet]; if (kRow >= 0) kRow += numberStaticRows_; int j = startSet_[iSet]; while (j >= 0) { double value = solution[j]; if (value) { for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) { int iRow = row_[k]; rhs[iRow] -= element_[k] * value; } if (kRow >= 0) rhs[kRow] -= value; } j = next_[j]; //onto next in set } } delete [] solution; } else { // bounds ClpSimplex::Status iStatus; for (int iSet = 0; iSet < numberSets_; iSet++) { int kRow = toIndex_[iSet]; if (kRow < 0) { int iColumn = keyVariable_[iSet]; if (iColumn < maximumGubColumns_) { // key is not treated as basic double b = 0.0; // key is structural - where is slack iStatus = getStatus(iSet); assert (iStatus != ClpSimplex::basic); if (iStatus == ClpSimplex::atLowerBound) b = lowerSet_[iSet]; else b = upperSet_[iSet]; if (b) { objectiveOffset += b * cost_[iColumn]; for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) { int iRow = row_[j]; rhs[iRow] -= element_[j] * b; } } } } } } if (fabs(model->objectiveOffset() - objectiveOffset_ - objectiveOffset) > 1.0e-1) printf("old offset %g, true %g\n", model->objectiveOffset(), objectiveOffset_ - objectiveOffset); for (iRow = 0; iRow < numberRows; iRow++) { if (fabs(rhs[iRow] - rhsOffset_[iRow]) > 1.0e-3) printf("** bad effective %d - true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]); } CoinMemcpyN(rhs, numberRows, saveE); delete [] rhs; } #endif if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) { int numberRows = model->numberRows(); int iSet; CoinZeroN(rhsOffset_, numberRows); // do ones at bounds before gub const double * smallSolution = model->solutionRegion(); const double * element = matrix_->getElements(); const int * row = matrix_->getIndices(); const CoinBigIndex * startColumn = matrix_->getVectorStarts(); const int * length = matrix_->getVectorLengths(); int iColumn; double objectiveOffset = 0.0; for (iColumn = 0; iColumn < firstDynamic_; iColumn++) { if (model->getStatus(iColumn) != ClpSimplex::basic) { double value = smallSolution[iColumn]; for (CoinBigIndex j = startColumn[iColumn]; j < startColumn[iColumn] + length[iColumn]; j++) { int jRow = row[j]; rhsOffset_[jRow] -= value * element[j]; } } } if (columnLower_ || columnUpper_) { double * solution = new double [numberGubColumns_]; for (iSet = 0; iSet < numberSets_; iSet++) { int j = startSet_[iSet]; while (j >= 0) { double value = 0.0; if (getDynamicStatus(j) != inSmall) { if (getDynamicStatus(j) == atLowerBound) { if (columnLower_) value = columnLower_[j]; } else if (getDynamicStatus(j) == atUpperBound) { value = columnUpper_[j]; assert (value<1.0e30); } else if (getDynamicStatus(j) == soloKey) { value = keyValue(iSet); } objectiveOffset += value * cost_[j]; } solution[j] = value; j = next_[j]; //onto next in set } } // ones in gub and in small problem for (iColumn = firstDynamic_; iColumn < firstAvailable_; iColumn++) { if (model_->getStatus(iColumn) != ClpSimplex::basic) { int jFull = id_[iColumn-firstDynamic_]; solution[jFull] = smallSolution[iColumn]; } } for (iSet = 0; iSet < numberSets_; iSet++) { int kRow = toIndex_[iSet]; if (kRow >= 0) kRow += numberStaticRows_; int j = startSet_[iSet]; while (j >= 0) { //? use DynamicStatus status = getDynamicStatus(j); double value = solution[j]; if (value) { for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) { int iRow = row_[k]; rhsOffset_[iRow] -= element_[k] * value; } if (kRow >= 0) rhsOffset_[kRow] -= value; } j = next_[j]; //onto next in set } } delete [] solution; } else { // bounds ClpSimplex::Status iStatus; for (int iSet = 0; iSet < numberSets_; iSet++) { int kRow = toIndex_[iSet]; if (kRow < 0) { int iColumn = keyVariable_[iSet]; if (iColumn < maximumGubColumns_) { // key is not treated as basic double b = 0.0; // key is structural - where is slack iStatus = getStatus(iSet); assert (iStatus != ClpSimplex::basic); if (iStatus == ClpSimplex::atLowerBound) b = lowerSet_[iSet]; else b = upperSet_[iSet]; if (b) { objectiveOffset += b * cost_[iColumn]; for (CoinBigIndex j = startColumn_[iColumn]; j < startColumn_[iColumn+1]; j++) { int iRow = row_[j]; rhsOffset_[iRow] -= element_[j] * b; } } } } } } model->setObjectiveOffset(objectiveOffset_ - objectiveOffset); #ifdef CLP_DEBUG if (saveE) { int iRow; for (iRow = 0; iRow < numberRows; iRow++) { if (fabs(saveE[iRow] - rhsOffset_[iRow]) > 1.0e-3) printf("** %d - old eff %g new %g\n", iRow, saveE[iRow], rhsOffset_[iRow]); } delete [] saveE; } #endif lastRefresh_ = model->numberIterations(); } } return rhsOffset_; } /* update information for a pivot (and effective rhs) */ int ClpDynamicMatrix::updatePivot(ClpSimplex * model, double oldInValue, double oldOutValue) { // now update working model int sequenceIn = model->sequenceIn(); int sequenceOut = model->sequenceOut(); int numberColumns = model->numberColumns(); if (sequenceIn != sequenceOut && sequenceIn < numberColumns) backToPivotRow_[sequenceIn] = model->pivotRow(); if (sequenceIn >= firstDynamic_ && sequenceIn < numberColumns) { int bigSequence = id_[sequenceIn-firstDynamic_]; if (getDynamicStatus(bigSequence) != inSmall) { firstAvailable_++; setDynamicStatus(bigSequence, inSmall); } } // make sure slack is synchronized if (sequenceIn >= numberColumns + numberStaticRows_) { int iDynamic = sequenceIn - numberColumns - numberStaticRows_; int iSet = fromIndex_[iDynamic]; setStatus(iSet, model->getStatus(sequenceIn)); } if (sequenceOut >= numberColumns + numberStaticRows_) { int iDynamic = sequenceOut - numberColumns - numberStaticRows_; int iSet = fromIndex_[iDynamic]; // out may have gone through barrier - so check double valueOut = model->lowerRegion()[sequenceOut]; if (fabs(valueOut - static_cast (lowerSet_[iSet])) < fabs(valueOut - static_cast (upperSet_[iSet]))) setStatus(iSet, ClpSimplex::atLowerBound); else setStatus(iSet, ClpSimplex::atUpperBound); if (lowerSet_[iSet] == upperSet_[iSet]) setStatus(iSet, ClpSimplex::isFixed); #if 0 if (getStatus(iSet) != model->getStatus(sequenceOut)) printf("** set %d status %d, var status %d\n", iSet, getStatus(iSet), model->getStatus(sequenceOut)); #endif } ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue); #ifdef CLP_DEBUG char * inSmall = new char [numberGubColumns_]; memset(inSmall, 0, numberGubColumns_); const double * solution = model->solutionRegion(); for (int i = 0; i < numberGubColumns_; i++) if (getDynamicStatus(i) == ClpDynamicMatrix::inSmall) inSmall[i] = 1; for (int i = firstDynamic_; i < firstAvailable_; i++) { int k = id_[i-firstDynamic_]; inSmall[k] = 0; //if (k>=23289&&k<23357&&solution[i]) //printf("var %d (in small %d) has value %g\n",k,i,solution[i]); } for (int i = 0; i < numberGubColumns_; i++) assert (!inSmall[i]); delete [] inSmall; #ifndef NDEBUG for (int i = 0; i < numberActiveSets_; i++) { int iSet = fromIndex_[i]; assert (toIndex_[iSet] == i); //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_)) //printf("*** set %d status %d, var status %d\n",iSet, // getStatus(iSet),model->getRowStatus(i+numberStaticRows_)); //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet)); //if (iSet==1035) { //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_], // model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]); //} } #endif #endif if (numberStaticRows_+numberActiveSets_numberRows()) return 0; else return 1; } /* utility dual function for dealing with dynamic constraints mode=n see ClpGubMatrix.hpp for definition Remember to update here when settled down */ void ClpDynamicMatrix::dualExpanded(ClpSimplex * model, CoinIndexedVector * /*array*/, double * /*other*/, int mode) { switch (mode) { // modify costs before transposeUpdate case 0: break; // create duals for key variables (without check on dual infeasible) case 1: break; // as 1 but check slacks and compute djs case 2: { // do pivots int * pivotVariable = model->pivotVariable(); int numberRows = numberStaticRows_ + numberActiveSets_; int numberColumns = model->numberColumns(); for (int iRow = 0; iRow < numberRows; iRow++) { int iPivot = pivotVariable[iRow]; if (iPivot < numberColumns) backToPivotRow_[iPivot] = iRow; } if (noCheck_ >= 0) { if (infeasibilityWeight_ != model_->infeasibilityCost()) { // don't bother checking sumDualInfeasibilities_ = 100.0; numberDualInfeasibilities_ = 1; sumOfRelaxedDualInfeasibilities_ = 100.0; return; } } // In theory we should subtract out ones we have done but .... // If key slack then dual 0.0 // If not then slack could be dual infeasible // dj for key is zero so that defines dual on set int i; double * dual = model->dualRowSolution(); double dualTolerance = model->dualTolerance(); double relaxedTolerance = dualTolerance; // we can't really trust infeasibilities if there is dual error double error = CoinMin(1.0e-2, model->largestDualError()); // allow tolerance at least slightly bigger than standard relaxedTolerance = relaxedTolerance + error; // but we will be using difference relaxedTolerance -= dualTolerance; sumDualInfeasibilities_ = 0.0; numberDualInfeasibilities_ = 0; sumOfRelaxedDualInfeasibilities_ = 0.0; for (i = 0; i < numberSets_; i++) { double value = 0.0; int gubRow = toIndex_[i]; if (gubRow < 0) { int kColumn = keyVariable_[i]; if (kColumn < maximumGubColumns_) { // dj without set value = cost_[kColumn]; for (CoinBigIndex j = startColumn_[kColumn]; j < startColumn_[kColumn+1]; j++) { int iRow = row_[j]; value -= dual[iRow] * element_[j]; } double infeasibility = 0.0; if (getStatus(i) == ClpSimplex::atLowerBound) { if (-value > dualTolerance) infeasibility = -value - dualTolerance; } else if (getStatus(i) == ClpSimplex::atUpperBound) { if (value > dualTolerance) infeasibility = value - dualTolerance; } if (infeasibility > 0.0) { sumDualInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += infeasibility; numberDualInfeasibilities_ ++; } } } else { value = dual[gubRow+numberStaticRows_]; } // Now subtract out from all int k = startSet_[i]; while (k >= 0) { if (getDynamicStatus(k) != inSmall) { double djValue = cost_[k] - value; for (CoinBigIndex j = startColumn_[k]; j < startColumn_[k+1]; j++) { int iRow = row_[j]; djValue -= dual[iRow] * element_[j]; } double infeasibility = 0.0; if (getDynamicStatus(k) == atLowerBound) { if (djValue < -dualTolerance) infeasibility = -djValue - dualTolerance; } else if (getDynamicStatus(k) == atUpperBound) { // at upper bound if (djValue > dualTolerance) infeasibility = djValue - dualTolerance; } if (infeasibility > 0.0) { sumDualInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += infeasibility; numberDualInfeasibilities_ ++; } } k = next_[k]; //onto next in set } } } infeasibilityWeight_ = -1.0; break; // Report on infeasibilities of key variables case 3: { model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_); model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_); model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_); } break; // modify costs before transposeUpdate for partial pricing case 4: break; } } /* general utility function for dealing with dynamic constraints mode=n see ClpGubMatrix.hpp for definition Remember to update here when settled down */ int ClpDynamicMatrix::generalExpanded(ClpSimplex * model, int mode, int &number) { int returnCode = 0; #if 0 //ndef NDEBUG { int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); int * pivotVariable = model->pivotVariable(); if (pivotVariable&&model->numberIterations()) { for (int i=numberStaticRows_+numberActiveSets_;inumberRows()]; rhsOffset(model, true); } int i; int numberBasic = number; int numberColumns = model->numberColumns(); // Use different array so can build from true pivotVariable_ //int * pivotVariable = model->pivotVariable(); int * pivotVariable = model->rowArray(0)->getIndices(); for (i = 0; i < numberColumns; i++) { if (model->getColumnStatus(i) == ClpSimplex::basic) pivotVariable[numberBasic++] = i; } number = numberBasic; } break; // Do initial extra rows + maximum basic case 2: { number = model->numberRows(); } break; // Before normal replaceColumn case 3: { if (numberActiveSets_ + numberStaticRows_ == model_->numberRows()) { // no space - re-factorize returnCode = 4; number = -1; // say no need for normal replaceColumn } } break; // To see if can dual or primal case 4: { returnCode = 1; } break; // save status case 5: { memcpy(status_+numberSets_,status_,numberSets_); memcpy(status_+2*numberSets_,&numberActiveSets_,sizeof(int)); memcpy(dynamicStatus_+maximumGubColumns_, dynamicStatus_,maximumGubColumns_); } break; // restore status case 6: { memcpy(status_,status_+numberSets_,numberSets_); memcpy(&numberActiveSets_,status_+2*numberSets_,sizeof(int)); memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_, maximumGubColumns_); initialProblem(); } break; // unflag all variables case 8: { for (int i = 0; i < numberGubColumns_; i++) { if (flagged(i)) { unsetFlagged(i); returnCode++; } } } break; // redo costs in primal case 9: { double * cost = model->costRegion(); double * solution = model->solutionRegion(); double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); int i; bool doCosts = (number & 4) != 0; bool doBounds = (number & 1) != 0; for ( i = firstDynamic_; i < firstAvailable_; i++) { int jColumn = id_[i-firstDynamic_]; if (doBounds) { if (!columnLower_ && !columnUpper_) { columnLower[i] = 0.0; columnUpper[i] = COIN_DBL_MAX; } else { if (columnLower_) columnLower[i] = columnLower_[jColumn]; else columnLower[i] = 0.0; if (columnUpper_) columnUpper[i] = columnUpper_[jColumn]; else columnUpper[i] = COIN_DBL_MAX; } } if (doCosts) { cost[i] = cost_[jColumn]; // Original bounds if (model->nonLinearCost()) model->nonLinearCost()->setOne(i, solution[i], this->columnLower(jColumn), this->columnUpper(jColumn), cost_[jColumn]); } } // and active sets for (i = 0; i < numberActiveSets_; i++ ) { int iSet = fromIndex_[i]; int iSequence = lastDynamic_ + numberStaticRows_ + i; if (doBounds) { if (lowerSet_[iSet] > -1.0e20) columnLower[iSequence] = lowerSet_[iSet]; else columnLower[iSequence] = -COIN_DBL_MAX; if (upperSet_[iSet] < 1.0e20) columnUpper[iSequence] = upperSet_[iSet]; else columnUpper[iSequence] = COIN_DBL_MAX; } if (doCosts) { if (model->nonLinearCost()) { double trueLower; if (lowerSet_[iSet] > -1.0e20) trueLower = lowerSet_[iSet]; else trueLower = -COIN_DBL_MAX; double trueUpper; if (upperSet_[iSet] < 1.0e20) trueUpper = upperSet_[iSet]; else trueUpper = COIN_DBL_MAX; model->nonLinearCost()->setOne(iSequence, solution[iSequence], trueLower, trueUpper, 0.0); } } } } break; // return 1 if there may be changing bounds on variable (column generation) case 10: { // return 1 as bounds on rhs will change returnCode = 1; } break; // make sure set is clean case 7: { // first flag if (number >= firstDynamic_ && number < lastDynamic_) { int sequence = id_[number-firstDynamic_]; setFlagged(sequence); //model->clearFlagged(number); } else if (number>=model_->numberColumns()+numberStaticRows_) { // slack int iSet = fromIndex_[number-model_->numberColumns()- numberStaticRows_]; setFlaggedSlack(iSet); //model->clearFlagged(number); } } case 11: { //int sequenceIn = model->sequenceIn(); if (number >= firstDynamic_ && number < lastDynamic_) { //assert (number == model->sequenceIn()); // take out variable (but leave key) double * cost = model->costRegion(); double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); double * solution = model->solutionRegion(); int * length = matrix_->getMutableVectorLengths(); #ifndef NDEBUG if (length[number]) { int * row = matrix_->getMutableIndices(); CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); int iRow = row[startColumn[number] + length[number] - 1]; int which = iRow - numberStaticRows_; assert (which >= 0); int iSet = fromIndex_[which]; assert (toIndex_[iSet] == which); } #endif // no need firstAvailable_--; solution[firstAvailable_] = 0.0; cost[firstAvailable_] = 0.0; length[firstAvailable_] = 0; model->nonLinearCost()->setOne(firstAvailable_, 0.0, 0.0, COIN_DBL_MAX, 0.0); model->setStatus(firstAvailable_, ClpSimplex::atLowerBound); columnLower[firstAvailable_] = 0.0; columnUpper[firstAvailable_] = COIN_DBL_MAX; // not really in small problem int iBig = id_[number-firstDynamic_]; if (model->getStatus(number) == ClpSimplex::atLowerBound) { setDynamicStatus(iBig, atLowerBound); if (columnLower_) modifyOffset(number, columnLower_[iBig]); } else { setDynamicStatus(iBig, atUpperBound); modifyOffset(number, columnUpper_[iBig]); } } else if (number>=model_->numberColumns()+numberStaticRows_) { // slack int iSet = fromIndex_[number-model_->numberColumns()- numberStaticRows_]; printf("what now - set %d\n",iSet); } } break; default: break; } return returnCode; } /* Purely for column generation and similar ideas. Allows matrix and any bounds or costs to be updated (sensibly). Returns non-zero if any changes. */ int ClpDynamicMatrix::refresh(ClpSimplex * model) { // If at beginning then create initial problem if (firstDynamic_ == firstAvailable_) { initialProblem(); return 1; } else if (!model->nonLinearCost()) { // will be same as last time return 1; } #ifndef NDEBUG { int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); int * pivotVariable = model->pivotVariable(); for (int i=numberStaticRows_+numberActiveSets_;inumberColumns(); double * element = matrix_->getMutableElements(); int * row = matrix_->getMutableIndices(); CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); int * length = matrix_->getMutableVectorLengths(); double * cost = model->costRegion(); double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); CoinBigIndex numberElements = startColumn[firstDynamic_]; // first just do lookup and basic stuff int currentNumber = firstAvailable_; firstAvailable_ = firstDynamic_; double objectiveChange = 0.0; double * solution = model->solutionRegion(); int currentNumberActiveSets = numberActiveSets_; for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) { int iRow = row[startColumn[iColumn] + length[iColumn] - 1]; int which = iRow - numberStaticRows_; assert (which >= 0); int iSet = fromIndex_[which]; assert (toIndex_[iSet] == which); if (model->getStatus(iColumn) == ClpSimplex::basic) { active[which]++; // may as well make key keyVariable_[iSet] = id_[iColumn-firstDynamic_]; } } int i; numberActiveSets_ = 0; int numberDeleted = 0; for (i = 0; i < currentNumberActiveSets; i++) { int iRow = i + numberStaticRows_; int numberActive = active[i]; int iSet = fromIndex_[i]; if (model->getRowStatus(iRow) == ClpSimplex::basic) { numberActive++; // may as well make key keyVariable_[iSet] = maximumGubColumns_ + iSet; } if (numberActive > 1) { // keep active[i] = numberActiveSets_; numberActiveSets_++; } else { active[i] = -1; } } for (iColumn = firstDynamic_; iColumn < currentNumber; iColumn++) { int iRow = row[startColumn[iColumn] + length[iColumn] - 1]; int which = iRow - numberStaticRows_; int jColumn = id_[iColumn-firstDynamic_]; if (model->getStatus(iColumn) == ClpSimplex::atLowerBound || model->getStatus(iColumn) == ClpSimplex::atUpperBound) { double value = solution[iColumn]; bool toLowerBound = true; assert (jColumn>=0);assert (iColumn>=0); if (columnUpper_) { if (!columnLower_) { if (fabs(value - columnUpper_[jColumn]) < fabs(value)) toLowerBound = false; } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) { toLowerBound = false; } } if (toLowerBound) { // throw out to lower bound if (columnLower_) { setDynamicStatus(jColumn, atLowerBound); // treat solution as if exactly at a bound double value = columnLower[iColumn]; objectiveChange += cost[iColumn] * value; } else { setDynamicStatus(jColumn, atLowerBound); } } else { // throw out to upper bound setDynamicStatus(jColumn, atUpperBound); double value = columnUpper[iColumn]; objectiveChange += cost[iColumn] * value; } numberDeleted++; } else { assert(model->getStatus(iColumn) != ClpSimplex::superBasic); // deal with later int iPut = active[which]; if (iPut >= 0) { // move id_[firstAvailable_-firstDynamic_] = jColumn; int numberThis = startColumn_[jColumn+1] - startColumn_[jColumn] + 1; length[firstAvailable_] = numberThis; cost[firstAvailable_] = cost[iColumn]; columnLower[firstAvailable_] = columnLower[iColumn]; columnUpper[firstAvailable_] = columnUpper[iColumn]; model->nonLinearCost()->setOne(firstAvailable_, solution[iColumn], 0.0, COIN_DBL_MAX, cost_[jColumn]); CoinBigIndex base = startColumn_[jColumn]; for (int j = 0; j < numberThis - 1; j++) { row[numberElements] = row_[base+j]; element[numberElements++] = element_[base+j]; } row[numberElements] = iPut + numberStaticRows_; element[numberElements++] = 1.0; model->setStatus(firstAvailable_, model->getStatus(iColumn)); solution[firstAvailable_] = solution[iColumn]; firstAvailable_++; startColumn[firstAvailable_] = numberElements; } } } // zero solution CoinZeroN(solution + firstAvailable_, currentNumber - firstAvailable_); // zero cost CoinZeroN(cost + firstAvailable_, currentNumber - firstAvailable_); // zero lengths CoinZeroN(length + firstAvailable_, currentNumber - firstAvailable_); for ( iColumn = firstAvailable_; iColumn < currentNumber; iColumn++) { model->nonLinearCost()->setOne(iColumn, 0.0, 0.0, COIN_DBL_MAX, 0.0); model->setStatus(iColumn, ClpSimplex::atLowerBound); columnLower[iColumn] = 0.0; columnUpper[iColumn] = COIN_DBL_MAX; } // move rhs and set rest to infinite numberActiveSets_ = 0; for (i = 0; i < currentNumberActiveSets; i++) { int iSet = fromIndex_[i]; assert (toIndex_[iSet] == i); int iRow = i + numberStaticRows_; int iPut = active[i]; if (iPut >= 0) { int in = iRow + numberColumns; int out = iPut + numberColumns + numberStaticRows_; solution[out] = solution[in]; columnLower[out] = columnLower[in]; columnUpper[out] = columnUpper[in]; cost[out] = cost[in]; double trueLower; if (lowerSet_[iSet] > -1.0e20) trueLower = lowerSet_[iSet]; else trueLower = -COIN_DBL_MAX; double trueUpper; if (upperSet_[iSet] < 1.0e20) trueUpper = upperSet_[iSet]; else trueUpper = COIN_DBL_MAX; model->nonLinearCost()->setOne(out, solution[out], trueLower, trueUpper, 0.0); model->setStatus(out, model->getStatus(in)); toIndex_[iSet] = numberActiveSets_; rhsOffset_[numberActiveSets_+numberStaticRows_] = rhsOffset_[i+numberStaticRows_]; fromIndex_[numberActiveSets_++] = fromIndex_[i]; } else { // adjust offset // put out as key int jColumn = keyVariable_[iSet]; toIndex_[iSet] = -1; if (jColumn < maximumGubColumns_) { setDynamicStatus(jColumn, soloKey); double value = keyValue(iSet); objectiveChange += cost_[jColumn] * value; modifyOffset(jColumn, -value); } } } model->setObjectiveOffset(model->objectiveOffset() - objectiveChange); delete [] active; for (i = numberActiveSets_; i < currentNumberActiveSets; i++) { int iSequence = i + numberStaticRows_ + numberColumns; solution[iSequence] = 0.0; columnLower[iSequence] = -COIN_DBL_MAX; columnUpper[iSequence] = COIN_DBL_MAX; cost[iSequence] = 0.0; model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence], columnUpper[iSequence], 0.0); model->setStatus(iSequence, ClpSimplex::basic); rhsOffset_[i+numberStaticRows_] = 0.0; } if (currentNumberActiveSets != numberActiveSets_ || numberDeleted) { // deal with pivotVariable int * pivotVariable = model->pivotVariable(); int sequence = firstDynamic_; int iRow = 0; int base1 = firstDynamic_; int base2 = lastDynamic_; int base3 = numberColumns + numberStaticRows_; int numberRows = numberStaticRows_ + currentNumberActiveSets; while (sequence < firstAvailable_) { int iPivot = pivotVariable[iRow]; while (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) { iPivot = pivotVariable[++iRow]; } pivotVariable[iRow++] = sequence; sequence++; } // move normal basic ones down int iPut = iRow; for (; iRow < numberRows; iRow++) { int iPivot = pivotVariable[iRow]; if (iPivot < base1 || (iPivot >= base2 && iPivot < base3)) { pivotVariable[iPut++] = iPivot; } } // and basic others for (i = 0; i < numberActiveSets_; i++) { if (model->getRowStatus(i + numberStaticRows_) == ClpSimplex::basic) { pivotVariable[iPut++] = i + base3; } } if (iPut= 0) { printf("var %d status %d ", jj, dynamicStatus_[jj]); jj = next_[jj]; } printf("\n"); #endif #ifdef CLP_DEBUG { // debug coding to analyze set int i; int kSet = -6; if (kSet >= 0) { int * back = new int [numberGubColumns_]; for (i = 0; i < numberGubColumns_; i++) back[i] = -1; for (i = firstDynamic_; i < firstAvailable_; i++) back[id_[i-firstDynamic_]] = i; int sequence = startSet_[kSet]; if (toIndex_[kSet] < 0) { printf("Not in - Set %d - slack status %d, key %d\n", kSet, status_[kSet], keyVariable_[kSet]); while (sequence >= 0) { printf("( %d status %d ) ", sequence, getDynamicStatus(sequence)); sequence = next_[sequence]; } } else { int iRow = numberStaticRows_ + toIndex_[kSet]; printf("In - Set %d - slack status %d, key %d offset %g slack %g\n", kSet, status_[kSet], keyVariable_[kSet], rhsOffset_[iRow], model->solutionRegion(0)[iRow]); while (sequence >= 0) { int iBack = back[sequence]; printf("( %d status %d value %g) ", sequence, getDynamicStatus(sequence), model->solutionRegion()[iBack]); sequence = next_[sequence]; } } printf("\n"); delete [] back; } } #endif int n = numberActiveSets_; for (i = 0; i < numberSets_; i++) { if (toIndex_[i] < 0) { //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]); n++; } int k=0; for (int j=startSet_[i];jspecialOptions() != 16) { ClpPackedMatrix::times(scalar, x, y); } else { int iRow; const double * element = matrix_->getElements(); const int * row = matrix_->getIndices(); const CoinBigIndex * startColumn = matrix_->getVectorStarts(); const int * length = matrix_->getVectorLengths(); int * pivotVariable = model_->pivotVariable(); for (iRow = 0; iRow < numberStaticRows_ + numberActiveSets_; iRow++) { y[iRow] -= scalar * rhsOffset_[iRow]; int iColumn = pivotVariable[iRow]; if (iColumn < lastDynamic_) { CoinBigIndex j; double value = scalar * x[iColumn]; if (value) { for (j = startColumn[iColumn]; j < startColumn[iColumn] + length[iColumn]; j++) { int jRow = row[j]; y[jRow] += value * element[j]; } } } } } } // Modifies rhs offset void ClpDynamicMatrix::modifyOffset(int sequence, double amount) { if (amount) { assert (rhsOffset_); CoinBigIndex j; for (j = startColumn_[sequence]; j < startColumn_[sequence+1]; j++) { int iRow = row_[j]; rhsOffset_[iRow] += amount * element_[j]; } } } // Gets key value when none in small double ClpDynamicMatrix::keyValue(int iSet) const { double value = 0.0; if (toIndex_[iSet] < 0) { int key = keyVariable_[iSet]; if (key < maximumGubColumns_) { if (getStatus(iSet) == ClpSimplex::atLowerBound) value = lowerSet_[iSet]; else value = upperSet_[iSet]; int numberKey = 0; int j = startSet_[iSet]; while (j >= 0) { DynamicStatus status = getDynamicStatus(j); assert (status != inSmall); if (status == soloKey) { numberKey++; } else if (status == atUpperBound) { value -= columnUpper_[j]; } else if (columnLower_) { value -= columnLower_[j]; } j = next_[j]; //onto next in set } assert (numberKey == 1); } else { int j = startSet_[iSet]; while (j >= 0) { DynamicStatus status = getDynamicStatus(j); assert (status != inSmall); assert (status != soloKey); if (status == atUpperBound) { value += columnUpper_[j]; } else if (columnLower_) { value += columnLower_[j]; } j = next_[j]; //onto next in set } #if 0 // slack must be feasible double oldValue=value; value = CoinMax(value,lowerSet_[iSet]); value = CoinMin(value,upperSet_[iSet]); if (value!=oldValue) printf("using %g (not %g) for slack on set %d (%g,%g)\n", value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]); #endif } } return value; } // Switches off dj checking each factorization (for BIG models) void ClpDynamicMatrix::switchOffCheck() { noCheck_ = 0; infeasibilityWeight_ = 0.0; } /* Creates a variable. This is called after partial pricing and may modify matrix. May update bestSequence. */ void ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence) { int numberRows = model->numberRows(); int slackOffset = lastDynamic_ + numberRows; int structuralOffset = slackOffset + numberSets_; int bestSequence2 = savedBestSequence_ - structuralOffset; if (bestSequence >= slackOffset) { double * columnLower = model->lowerRegion(); double * columnUpper = model->upperRegion(); double * solution = model->solutionRegion(); double * reducedCost = model->djRegion(); const double * duals = model->dualRowSolution(); if (toIndex_[savedBestSet_] < 0) { // need to put key into basis int newRow = numberActiveSets_ + numberStaticRows_; model->dualRowSolution()[newRow] = savedBestGubDual_; double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set toIndex_[savedBestSet_] = numberActiveSets_; fromIndex_[numberActiveSets_++] = savedBestSet_; int iSequence = lastDynamic_ + newRow; // we need to get lower and upper correct double shift = 0.0; int j = startSet_[savedBestSet_]; while (j >= 0) { if (getDynamicStatus(j) == atUpperBound) shift += columnUpper_[j]; else if (getDynamicStatus(j) == atLowerBound && columnLower_) shift += columnLower_[j]; j = next_[j]; //onto next in set } if (lowerSet_[savedBestSet_] > -1.0e20) columnLower[iSequence] = lowerSet_[savedBestSet_]; else columnLower[iSequence] = -COIN_DBL_MAX; if (upperSet_[savedBestSet_] < 1.0e20) columnUpper[iSequence] = upperSet_[savedBestSet_]; else columnUpper[iSequence] = COIN_DBL_MAX; #ifdef CLP_DEBUG if (model_->logLevel() == 63) { printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n", bestSequence2, savedBestSet_, keyVariable_[savedBestSet_], columnLower[iSequence], columnUpper[iSequence], valueOfKey); int j = startSet_[savedBestSet_]; while (j >= 0) { if (getDynamicStatus(j) == atUpperBound) printf("%d atup ", j); else if (getDynamicStatus(j) == atLowerBound) printf("%d atlo ", j); else if (getDynamicStatus(j) == soloKey) printf("%d solo ", j); else abort(); j = next_[j]; //onto next in set } printf("\n"); } #endif if (keyVariable_[savedBestSet_] < maximumGubColumns_) { // slack not key model_->pivotVariable()[newRow] = firstAvailable_; backToPivotRow_[firstAvailable_] = newRow; model->setStatus(iSequence, getStatus(savedBestSet_)); model->djRegion()[iSequence] = savedBestGubDual_; solution[iSequence] = valueOfKey; // create variable and pivot in int key = keyVariable_[savedBestSet_]; setDynamicStatus(key, inSmall); double * element = matrix_->getMutableElements(); int * row = matrix_->getMutableIndices(); CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); int * length = matrix_->getMutableVectorLengths(); CoinBigIndex numberElements = startColumn[firstAvailable_]; int numberThis = startColumn_[key+1] - startColumn_[key] + 1; if (numberElements + numberThis > numberElements_) { // need to redo numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); matrix_->reserve(lastDynamic_, numberElements_); element = matrix_->getMutableElements(); row = matrix_->getMutableIndices(); // these probably okay but be safe startColumn = matrix_->getMutableVectorStarts(); length = matrix_->getMutableVectorLengths(); } // already set startColumn[firstAvailable_]=numberElements; length[firstAvailable_] = numberThis; model->costRegion()[firstAvailable_] = cost_[key]; CoinBigIndex base = startColumn_[key]; for (int j = 0; j < numberThis - 1; j++) { row[numberElements] = row_[base+j]; element[numberElements++] = element_[base+j]; } row[numberElements] = newRow; element[numberElements++] = 1.0; id_[firstAvailable_-firstDynamic_] = key; model->setObjectiveOffset(model->objectiveOffset() + cost_[key]* valueOfKey); model->solutionRegion()[firstAvailable_] = valueOfKey; model->setStatus(firstAvailable_, ClpSimplex::basic); // ***** need to adjust effective rhs if (!columnLower_ && !columnUpper_) { columnLower[firstAvailable_] = 0.0; columnUpper[firstAvailable_] = COIN_DBL_MAX; } else { if (columnLower_) columnLower[firstAvailable_] = columnLower_[key]; else columnLower[firstAvailable_] = 0.0; if (columnUpper_) columnUpper[firstAvailable_] = columnUpper_[key]; else columnUpper[firstAvailable_] = COIN_DBL_MAX; } model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_], columnLower[firstAvailable_], columnUpper[firstAvailable_], cost_[key]); startColumn[firstAvailable_+1] = numberElements; reducedCost[firstAvailable_] = 0.0; modifyOffset(key, valueOfKey); rhsOffset_[newRow] = -shift; // sign? #ifdef CLP_DEBUG model->rowArray(1)->checkClear(); #endif // now pivot in unpack(model, model->rowArray(1), firstAvailable_); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(1)); double alpha = model->rowArray(1)->denseVector()[newRow]; int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(1), newRow, alpha); model->rowArray(1)->clear(); if (updateStatus) { if (updateStatus == 3) { // out of memory // increase space if not many iterations if (model->factorization()->pivots() < 0.5 * model->factorization()->maximumPivots() && model->factorization()->pivots() < 400) model->factorization()->areaFactor( model->factorization()->areaFactor() * 1.1); } else { printf("Bad returncode %d from replaceColumn\n", updateStatus); } bestSequence = -1; return; } // firstAvailable_ only finally updated if good pivot (in updatePivot) // otherwise it reverts to firstAvailableBefore_ firstAvailable_++; } else { // slack key model->setStatus(iSequence, ClpSimplex::basic); model->djRegion()[iSequence] = 0.0; solution[iSequence] = valueOfKey+shift; rhsOffset_[newRow] = -shift; // sign? } // correct slack model->costRegion()[iSequence] = 0.0; model->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence], columnUpper[iSequence], 0.0); } if (savedBestSequence_ >= structuralOffset) { // recompute dj and create double value = cost_[bestSequence2] - savedBestGubDual_; for (CoinBigIndex jBigIndex = startColumn_[bestSequence2]; jBigIndex < startColumn_[bestSequence2+1]; jBigIndex++) { int jRow = row_[jBigIndex]; value -= duals[jRow] * element_[jBigIndex]; } int gubRow = toIndex_[savedBestSet_] + numberStaticRows_; double * element = matrix_->getMutableElements(); int * row = matrix_->getMutableIndices(); CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); int * length = matrix_->getMutableVectorLengths(); CoinBigIndex numberElements = startColumn[firstAvailable_]; int numberThis = startColumn_[bestSequence2+1] - startColumn_[bestSequence2] + 1; if (numberElements + numberThis > numberElements_) { // need to redo numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); matrix_->reserve(lastDynamic_, numberElements_); element = matrix_->getMutableElements(); row = matrix_->getMutableIndices(); // these probably okay but be safe startColumn = matrix_->getMutableVectorStarts(); length = matrix_->getMutableVectorLengths(); } // already set startColumn[firstAvailable_]=numberElements; length[firstAvailable_] = numberThis; model->costRegion()[firstAvailable_] = cost_[bestSequence2]; CoinBigIndex base = startColumn_[bestSequence2]; for (int j = 0; j < numberThis - 1; j++) { row[numberElements] = row_[base+j]; element[numberElements++] = element_[base+j]; } row[numberElements] = gubRow; element[numberElements++] = 1.0; id_[firstAvailable_-firstDynamic_] = bestSequence2; //printf("best %d\n",bestSequence2); model->solutionRegion()[firstAvailable_] = 0.0; model->clearFlagged(firstAvailable_); if (!columnLower_ && !columnUpper_) { model->setStatus(firstAvailable_, ClpSimplex::atLowerBound); columnLower[firstAvailable_] = 0.0; columnUpper[firstAvailable_] = COIN_DBL_MAX; } else { DynamicStatus status = getDynamicStatus(bestSequence2); if (columnLower_) columnLower[firstAvailable_] = columnLower_[bestSequence2]; else columnLower[firstAvailable_] = 0.0; if (columnUpper_) columnUpper[firstAvailable_] = columnUpper_[bestSequence2]; else columnUpper[firstAvailable_] = COIN_DBL_MAX; if (status == atLowerBound) { solution[firstAvailable_] = columnLower[firstAvailable_]; model->setStatus(firstAvailable_, ClpSimplex::atLowerBound); } else { solution[firstAvailable_] = columnUpper[firstAvailable_]; model->setStatus(firstAvailable_, ClpSimplex::atUpperBound); } } model->setObjectiveOffset(model->objectiveOffset() + cost_[bestSequence2]* solution[firstAvailable_]); model->nonLinearCost()->setOne(firstAvailable_, solution[firstAvailable_], columnLower[firstAvailable_], columnUpper[firstAvailable_], cost_[bestSequence2]); bestSequence = firstAvailable_; // firstAvailable_ only updated if good pivot (in updatePivot) startColumn[firstAvailable_+1] = numberElements; //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_); reducedCost[bestSequence] = value; } else { // slack - row must just have been created assert (toIndex_[savedBestSet_] == numberActiveSets_ - 1); int newRow = numberStaticRows_ + numberActiveSets_ - 1; bestSequence = lastDynamic_ + newRow; reducedCost[bestSequence] = savedBestGubDual_; } } // clear for next iteration savedBestSequence_ = -1; } // Returns reduced cost of a variable double ClpDynamicMatrix::reducedCost(ClpSimplex * model, int sequence) const { int numberRows = model->numberRows(); int slackOffset = lastDynamic_ + numberRows; if (sequence < slackOffset) return model->djRegion()[sequence]; else return savedBestDj_; } // Does gub crash void ClpDynamicMatrix::gubCrash() { // Do basis - cheapest or slack if feasible int longestSet = 0; int iSet; for (iSet = 0; iSet < numberSets_; iSet++) { int n = 0; int j = startSet_[iSet]; while (j >= 0) { n++; j = next_[j]; } longestSet = CoinMax(longestSet, n); } double * upper = new double[longestSet+1]; double * cost = new double[longestSet+1]; double * lower = new double[longestSet+1]; double * solution = new double[longestSet+1]; int * back = new int[longestSet+1]; double tolerance = model_->primalTolerance(); double objectiveOffset = 0.0; for (iSet = 0; iSet < numberSets_; iSet++) { int iBasic = -1; double value = 0.0; // find cheapest int numberInSet = 0; int j = startSet_[iSet]; while (j >= 0) { if (!columnLower_) lower[numberInSet] = 0.0; else lower[numberInSet] = columnLower_[j]; if (!columnUpper_) upper[numberInSet] = COIN_DBL_MAX; else upper[numberInSet] = columnUpper_[j]; back[numberInSet++] = j; j = next_[j]; } CoinFillN(solution, numberInSet, 0.0); // and slack iBasic = numberInSet; solution[iBasic] = -value; lower[iBasic] = -upperSet_[iSet]; upper[iBasic] = -lowerSet_[iSet]; int kphase; if (value < lowerSet_[iSet] - tolerance || value > upperSet_[iSet] + tolerance) { // infeasible kphase = 0; // remember bounds are flipped so opposite to natural if (value < lowerSet_[iSet] - tolerance) cost[iBasic] = 1.0; else cost[iBasic] = -1.0; CoinZeroN(cost, numberInSet); double dualTolerance = model_->dualTolerance(); for (int iphase = kphase; iphase < 2; iphase++) { if (iphase) { cost[numberInSet] = 0.0; for (int j = 0; j < numberInSet; j++) cost[j] = cost_[back[j]]; } // now do one row lp bool improve = true; while (improve) { improve = false; double dual = cost[iBasic]; int chosen = -1; double best = dualTolerance; int way = 0; for (int i = 0; i <= numberInSet; i++) { double dj = cost[i] - dual; double improvement = 0.0; if (iphase || i < numberInSet) assert (solution[i] >= lower[i] && solution[i] <= upper[i]); if (dj > dualTolerance) improvement = dj * (solution[i] - lower[i]); else if (dj < -dualTolerance) improvement = dj * (solution[i] - upper[i]); if (improvement > best) { best = improvement; chosen = i; if (dj < 0.0) { way = 1; } else { way = -1; } } } if (chosen >= 0) { improve = true; // now see how far if (way > 0) { // incoming increasing so basic decreasing // if phase 0 then go to nearest bound double distance = upper[chosen] - solution[chosen]; double basicDistance; if (!iphase) { assert (iBasic == numberInSet); assert (solution[iBasic] > upper[iBasic]); basicDistance = solution[iBasic] - upper[iBasic]; } else { basicDistance = solution[iBasic] - lower[iBasic]; } // need extra coding for unbounded assert (CoinMin(distance, basicDistance) < 1.0e20); if (distance > basicDistance) { // incoming becomes basic solution[chosen] += basicDistance; if (!iphase) solution[iBasic] = upper[iBasic]; else solution[iBasic] = lower[iBasic]; iBasic = chosen; } else { // flip solution[chosen] = upper[chosen]; solution[iBasic] -= distance; } } else { // incoming decreasing so basic increasing // if phase 0 then go to nearest bound double distance = solution[chosen] - lower[chosen]; double basicDistance; if (!iphase) { assert (iBasic == numberInSet); assert (solution[iBasic] < lower[iBasic]); basicDistance = lower[iBasic] - solution[iBasic]; } else { basicDistance = upper[iBasic] - solution[iBasic]; } // need extra coding for unbounded - for now just exit if (CoinMin(distance, basicDistance) > 1.0e20) { printf("unbounded on set %d\n", iSet); iphase = 1; iBasic = numberInSet; break; } if (distance > basicDistance) { // incoming becomes basic solution[chosen] -= basicDistance; if (!iphase) solution[iBasic] = lower[iBasic]; else solution[iBasic] = upper[iBasic]; iBasic = chosen; } else { // flip solution[chosen] = lower[chosen]; solution[iBasic] += distance; } } if (!iphase) { if(iBasic < numberInSet) break; // feasible else if (solution[iBasic] >= lower[iBasic] && solution[iBasic] <= upper[iBasic]) break; // feasible (on flip) } } } } } // do solution i.e. bounds if (columnLower_ || columnUpper_) { for (int j = 0; j < numberInSet; j++) { if (j != iBasic) { objectiveOffset += solution[j] * cost[j]; if (columnLower_ && columnUpper_) { if (fabs(solution[j] - columnLower_[back[j]]) > fabs(solution[j] - columnUpper_[back[j]])) setDynamicStatus(back[j], atUpperBound); } else if (columnUpper_ && solution[j] > 0.0) { setDynamicStatus(back[j], atUpperBound); } else { setDynamicStatus(back[j], atLowerBound); assert(!solution[j]); } } } } // convert iBasic back and do bounds if (iBasic == numberInSet) { // slack basic setStatus(iSet, ClpSimplex::basic); iBasic = iSet + maximumGubColumns_; } else { iBasic = back[iBasic]; setDynamicStatus(iBasic, soloKey); // remember bounds flipped if (upper[numberInSet] == lower[numberInSet]) setStatus(iSet, ClpSimplex::isFixed); else if (solution[numberInSet] == upper[numberInSet]) setStatus(iSet, ClpSimplex::atLowerBound); else if (solution[numberInSet] == lower[numberInSet]) setStatus(iSet, ClpSimplex::atUpperBound); else abort(); } keyVariable_[iSet] = iBasic; } model_->setObjectiveOffset(objectiveOffset_ - objectiveOffset); delete [] lower; delete [] solution; delete [] upper; delete [] cost; delete [] back; // make sure matrix is in good shape matrix_->orderMatrix(); } // Populates initial matrix from dynamic status void ClpDynamicMatrix::initialProblem() { int iSet; double * element = matrix_->getMutableElements(); int * row = matrix_->getMutableIndices(); CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); int * length = matrix_->getMutableVectorLengths(); double * cost = model_->objective(); double * solution = model_->primalColumnSolution(); double * columnLower = model_->columnLower(); double * columnUpper = model_->columnUpper(); double * rowSolution = model_->primalRowSolution(); double * rowLower = model_->rowLower(); double * rowUpper = model_->rowUpper(); CoinBigIndex numberElements = startColumn[firstDynamic_]; firstAvailable_ = firstDynamic_; numberActiveSets_ = 0; for (iSet = 0; iSet < numberSets_; iSet++) { toIndex_[iSet] = -1; int numberActive = 0; int whichKey = -1; if (getStatus(iSet) == ClpSimplex::basic) { whichKey = maximumGubColumns_ + iSet; numberActive = 1; } else { whichKey = -1; } int j = startSet_[iSet]; while (j >= 0) { assert (getDynamicStatus(j) != soloKey || whichKey == -1); if (getDynamicStatus(j) == inSmall) { numberActive++; } else if (getDynamicStatus(j) == soloKey) { whichKey = j; numberActive++; } j = next_[j]; //onto next in set } if (numberActive > 1) { int iRow = numberActiveSets_ + numberStaticRows_; rowSolution[iRow] = 0.0; double lowerValue; if (lowerSet_[iSet] > -1.0e20) lowerValue = lowerSet_[iSet]; else lowerValue = -COIN_DBL_MAX; double upperValue; if (upperSet_[iSet] < 1.0e20) upperValue = upperSet_[iSet]; else upperValue = COIN_DBL_MAX; rowLower[iRow] = lowerValue; rowUpper[iRow] = upperValue; if (getStatus(iSet) == ClpSimplex::basic) { model_->setRowStatus(iRow, ClpSimplex::basic); rowSolution[iRow] = 0.0; } else if (getStatus(iSet) == ClpSimplex::atLowerBound) { model_->setRowStatus(iRow, ClpSimplex::atLowerBound); rowSolution[iRow] = lowerValue; } else { model_->setRowStatus(iRow, ClpSimplex::atUpperBound); rowSolution[iRow] = upperValue; } j = startSet_[iSet]; while (j >= 0) { DynamicStatus status = getDynamicStatus(j); if (status == inSmall) { int numberThis = startColumn_[j+1] - startColumn_[j] + 1; if (numberElements + numberThis > numberElements_) { // need to redo numberElements_ = CoinMax(3 * numberElements_ / 2, numberElements + numberThis); matrix_->reserve(lastDynamic_, numberElements_); element = matrix_->getMutableElements(); row = matrix_->getMutableIndices(); // these probably okay but be safe startColumn = matrix_->getMutableVectorStarts(); length = matrix_->getMutableVectorLengths(); } length[firstAvailable_] = numberThis; cost[firstAvailable_] = cost_[j]; CoinBigIndex base = startColumn_[j]; for (int k = 0; k < numberThis - 1; k++) { row[numberElements] = row_[base+k]; element[numberElements++] = element_[base+k]; } row[numberElements] = iRow; element[numberElements++] = 1.0; id_[firstAvailable_-firstDynamic_] = j; solution[firstAvailable_] = 0.0; model_->setStatus(firstAvailable_, ClpSimplex::basic); if (!columnLower_ && !columnUpper_) { columnLower[firstAvailable_] = 0.0; columnUpper[firstAvailable_] = COIN_DBL_MAX; } else { if (columnLower_) columnLower[firstAvailable_] = columnLower_[j]; else columnLower[firstAvailable_] = 0.0; if (columnUpper_) columnUpper[firstAvailable_] = columnUpper_[j]; else columnUpper[firstAvailable_] = COIN_DBL_MAX; if (status != atUpperBound) { solution[firstAvailable_] = columnLower[firstAvailable_]; } else { solution[firstAvailable_] = columnUpper[firstAvailable_]; } } firstAvailable_++; startColumn[firstAvailable_] = numberElements; } j = next_[j]; //onto next in set } model_->setRowStatus(numberActiveSets_ + numberStaticRows_, getStatus(iSet)); toIndex_[iSet] = numberActiveSets_; fromIndex_[numberActiveSets_++] = iSet; } else { // solo key bool needKey=false; if (numberActive) { if (whichKeyupperSet_[iSet]+1.0e-8) needKey=true; } } else { needKey = true; } if (needKey) { // all to lb then up some (slack/null if possible) int length=99999999; int which=-1; double sum=0.0; for (int iColumn=startSet_[iSet];iColumnstartColumn_[iColumn+1]-startColumn_[iColumn]) { which=iColumn; length=startColumn_[iColumn+1]-startColumn_[iColumn]; } } if (sum>lowerSet_[iSet]-1.0e-8) { // slack can be basic setStatus(iSet,ClpSimplex::basic); keyVariable_[iSet] = maximumGubColumns_ + iSet; } else { // use shortest setDynamicStatus(which,soloKey); keyVariable_[iSet] = which; setStatus(iSet,ClpSimplex::atLowerBound); } } } assert (toIndex_[iSet] >= 0 || whichKey >= 0); keyVariable_[iSet] = whichKey; } // clean up pivotVariable int numberColumns = model_->numberColumns(); int numberRows = model_->numberRows(); int * pivotVariable = model_->pivotVariable(); if (pivotVariable) { for (int i=0; igetRowStatus(i)!=ClpSimplex::basic) pivotVariable[i]=-1; else pivotVariable[i]=numberColumns+i; } for (int i=numberStaticRows_+numberActiveSets_;igetColumnStatus(i)==ClpSimplex::basic) { while(putcostRegion(); double * columnLower = model_->lowerRegion(); double * columnUpper = model_->upperRegion(); double * solution = model_->solutionRegion(); int numberRows = model_->numberRows(); for (int i = numberActiveSets_; i < numberRows-numberStaticRows_; i++) { int iSequence = i + numberStaticRows_ + numberColumns; solution[iSequence] = 0.0; columnLower[iSequence] = -COIN_DBL_MAX; columnUpper[iSequence] = COIN_DBL_MAX; cost[iSequence] = 0.0; model_->nonLinearCost()->setOne(iSequence, solution[iSequence], columnLower[iSequence], columnUpper[iSequence], 0.0); model_->setStatus(iSequence, ClpSimplex::basic); rhsOffset_[i+numberStaticRows_] = 0.0; } #if 0 for (int i=0;icolumnLower(); const double * columnUpperA = model_->columnUpper(); const double * costA = model_->objective(); const double * rowLowerA = model_->rowLower(); const double * rowUpperA = model_->rowUpper(); start[0]=0; numberElements=0; for (int i=0;i= 0) { if (startColumn_[j+1] - startColumn_[j] == numberEntries) { const int * row2 = row_ + startColumn_[j]; const double * element2 = element_ + startColumn_[j]; bool same = true; for (int k = 0; k < numberEntries; k++) { if (row[k] != row2[k] || element[k] != element2[k]) { same = false; break; } } if (same) { bool odd = false; if (cost != cost_[j]) odd = true; if (columnLower_ && lower != columnLower_[j]) odd = true; if (columnUpper_ && upper != columnUpper_[j]) odd = true; if (odd) { printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n", cost, lower, upper, cost_[j], columnLower_ ? columnLower_[j] : 0.0, columnUpper_ ? columnUpper_[j] : 1.0e100); } else { setDynamicStatus(j, status); return j; } } } j = next_[j]; } if (numberGubColumns_ == maximumGubColumns_ || startColumn_[numberGubColumns_] + numberEntries > maximumElements_) { CoinBigIndex j; int i; int put = 0; int numberElements = 0; CoinBigIndex start = 0; // compress - leave ones at ub and basic int * which = new int [numberGubColumns_]; for (i = 0; i < numberGubColumns_; i++) { CoinBigIndex end = startColumn_[i+1]; // what about ubs if column generation? if (getDynamicStatus(i) != atLowerBound) { // keep in for (j = start; j < end; j++) { row_[numberElements] = row_[j]; element_[numberElements++] = element_[j]; } startColumn_[put+1] = numberElements; cost_[put] = cost_[i]; if (columnLower_) columnLower_[put] = columnLower_[i]; if (columnUpper_) columnUpper_[put] = columnUpper_[i]; dynamicStatus_[put] = dynamicStatus_[i]; id_[put] = id_[i]; which[i] = put; put++; } else { // out which[i] = -1; } start = end; } // now redo startSet_ and next_ int * newNext = new int [maximumGubColumns_]; for (int jSet = 0; jSet < numberSets_; jSet++) { int sequence = startSet_[jSet]; while (which[sequence] < 0) { // out assert (next_[sequence] >= 0); sequence = next_[sequence]; } startSet_[jSet] = which[sequence]; int last = which[sequence]; while (next_[sequence] >= 0) { sequence = next_[sequence]; if(which[sequence] >= 0) { // keep newNext[last] = which[sequence]; last = which[sequence]; } } newNext[last] = -jSet - 1; } delete [] next_; next_ = newNext; delete [] which; abort(); } CoinBigIndex start = startColumn_[numberGubColumns_]; CoinMemcpyN(row, numberEntries, row_ + start); CoinMemcpyN(element, numberEntries, element_ + start); startColumn_[numberGubColumns_+1] = start + numberEntries; cost_[numberGubColumns_] = cost; if (columnLower_) columnLower_[numberGubColumns_] = lower; else assert (!lower); if (columnUpper_) columnUpper_[numberGubColumns_] = upper; else assert (upper > 1.0e20); setDynamicStatus(numberGubColumns_, status); // Do next_ j = startSet_[iSet]; startSet_[iSet] = numberGubColumns_; next_[numberGubColumns_] = j; numberGubColumns_++; return numberGubColumns_ - 1; } // Returns which set a variable is in int ClpDynamicMatrix::whichSet (int sequence) const { while (next_[sequence] >= 0) sequence = next_[sequence]; int iSet = - next_[sequence] - 1; return iSet; }