/* $Id$ */ // 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 #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 "ClpGubMatrix.hpp" //#include "ClpGubDynamicMatrix.hpp" #include "ClpMessage.hpp" //#define CLP_DEBUG //#define CLP_DEBUG_PRINT //############################################################################# // Constructors / Destructor / Assignment //############################################################################# //------------------------------------------------------------------- // Default Constructor //------------------------------------------------------------------- ClpGubMatrix::ClpGubMatrix() : ClpPackedMatrix() , sumDualInfeasibilities_(0.0) , sumPrimalInfeasibilities_(0.0) , sumOfRelaxedDualInfeasibilities_(0.0) , sumOfRelaxedPrimalInfeasibilities_(0.0) , infeasibilityWeight_(0.0) , start_(NULL) , end_(NULL) , lower_(NULL) , upper_(NULL) , status_(NULL) , saveStatus_(NULL) , savedKeyVariable_(NULL) , backward_(NULL) , backToPivotRow_(NULL) , changeCost_(NULL) , keyVariable_(NULL) , next_(NULL) , toIndex_(NULL) , fromIndex_(NULL) , model_(NULL) , numberDualInfeasibilities_(0) , numberPrimalInfeasibilities_(0) , noCheck_(-1) , numberSets_(0) , saveNumber_(0) , possiblePivotKey_(0) , gubSlackIn_(-1) , firstGub_(0) , lastGub_(0) , gubType_(0) { setType(16); } //------------------------------------------------------------------- // Copy constructor //------------------------------------------------------------------- ClpGubMatrix::ClpGubMatrix(const ClpGubMatrix &rhs) : ClpPackedMatrix(rhs) { numberSets_ = rhs.numberSets_; saveNumber_ = rhs.saveNumber_; possiblePivotKey_ = rhs.possiblePivotKey_; gubSlackIn_ = rhs.gubSlackIn_; start_ = ClpCopyOfArray(rhs.start_, numberSets_); end_ = ClpCopyOfArray(rhs.end_, numberSets_); lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); status_ = ClpCopyOfArray(rhs.status_, numberSets_); saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); int numberColumns = getNumCols(); backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); // find longest set int *longest = new int[numberSets_]; CoinZeroN(longest, numberSets_); int j; for (j = 0; j < numberColumns; j++) { int iSet = backward_[j]; if (iSet >= 0) longest[iSet]++; } int length = 0; for (j = 0; j < numberSets_; j++) length = CoinMax(length, longest[j]); next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; infeasibilityWeight_ = rhs.infeasibilityWeight_; numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; noCheck_ = rhs.noCheck_; firstGub_ = rhs.firstGub_; lastGub_ = rhs.lastGub_; gubType_ = rhs.gubType_; model_ = rhs.model_; } //------------------------------------------------------------------- // assign matrix (for space reasons) //------------------------------------------------------------------- ClpGubMatrix::ClpGubMatrix(CoinPackedMatrix *rhs) : ClpPackedMatrix(rhs) , sumDualInfeasibilities_(0.0) , sumPrimalInfeasibilities_(0.0) , sumOfRelaxedDualInfeasibilities_(0.0) , sumOfRelaxedPrimalInfeasibilities_(0.0) , infeasibilityWeight_(0.0) , start_(NULL) , end_(NULL) , lower_(NULL) , upper_(NULL) , status_(NULL) , saveStatus_(NULL) , savedKeyVariable_(NULL) , backward_(NULL) , backToPivotRow_(NULL) , changeCost_(NULL) , keyVariable_(NULL) , next_(NULL) , toIndex_(NULL) , fromIndex_(NULL) , model_(NULL) , numberDualInfeasibilities_(0) , numberPrimalInfeasibilities_(0) , noCheck_(-1) , numberSets_(0) , saveNumber_(0) , possiblePivotKey_(0) , gubSlackIn_(-1) , firstGub_(0) , lastGub_(0) , gubType_(0) { setType(16); } /* This takes over ownership (for space reasons) and is the real constructor*/ ClpGubMatrix::ClpGubMatrix(ClpPackedMatrix *matrix, int numberSets, const int *start, const int *end, const double *lower, const double *upper, const unsigned char *status) : ClpPackedMatrix(matrix->matrix()) , sumDualInfeasibilities_(0.0) , sumPrimalInfeasibilities_(0.0) , sumOfRelaxedDualInfeasibilities_(0.0) , sumOfRelaxedPrimalInfeasibilities_(0.0) , numberDualInfeasibilities_(0) , numberPrimalInfeasibilities_(0) , saveNumber_(0) , possiblePivotKey_(0) , gubSlackIn_(-1) { model_ = NULL; numberSets_ = numberSets; start_ = ClpCopyOfArray(start, numberSets_); end_ = ClpCopyOfArray(end, numberSets_); lower_ = ClpCopyOfArray(lower, numberSets_); upper_ = ClpCopyOfArray(upper, numberSets_); // Check valid and ordered int last = -1; int numberColumns = matrix_->getNumCols(); int numberRows = matrix_->getNumRows(); backward_ = new int[numberColumns]; backToPivotRow_ = new int[numberColumns]; changeCost_ = new double[numberRows + numberSets_]; keyVariable_ = new int[numberSets_]; // signal to need new ordering next_ = NULL; for (int iColumn = 0; iColumn < numberColumns; iColumn++) backward_[iColumn] = -1; int iSet; for (iSet = 0; iSet < numberSets_; iSet++) { // set key variable as slack keyVariable_[iSet] = iSet + numberColumns; if (start_[iSet] < 0 || start_[iSet] >= numberColumns) throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); if (end_[iSet] < 0 || end_[iSet] > numberColumns) throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); if (end_[iSet] <= start_[iSet]) throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix"); if (start_[iSet] < last) throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix"); last = end_[iSet]; int j; for (j = start_[iSet]; j < end_[iSet]; j++) backward_[j] = iSet; } // Find type of gub firstGub_ = numberColumns + 1; lastGub_ = -1; int i; for (i = 0; i < numberColumns; i++) { if (backward_[i] >= 0) { firstGub_ = CoinMin(firstGub_, i); lastGub_ = CoinMax(lastGub_, i); } } gubType_ = 0; // adjust lastGub_ if (lastGub_ > 0) lastGub_++; for (i = firstGub_; i < lastGub_; i++) { if (backward_[i] < 0) { gubType_ = 1; printf("interior non gub %d\n", i); break; } } if (status) { status_ = ClpCopyOfArray(status, numberSets_); } else { status_ = new unsigned char[numberSets_]; memset(status_, 0, numberSets_); int i; for (i = 0; i < numberSets_; i++) { // make slack key setStatus(i, ClpSimplex::basic); } } saveStatus_ = new unsigned char[numberSets_]; memset(saveStatus_, 0, numberSets_); savedKeyVariable_ = new int[numberSets_]; memset(savedKeyVariable_, 0, numberSets_ * sizeof(int)); noCheck_ = -1; infeasibilityWeight_ = 0.0; } ClpGubMatrix::ClpGubMatrix(const CoinPackedMatrix &rhs) : ClpPackedMatrix(rhs) , sumDualInfeasibilities_(0.0) , sumPrimalInfeasibilities_(0.0) , sumOfRelaxedDualInfeasibilities_(0.0) , sumOfRelaxedPrimalInfeasibilities_(0.0) , infeasibilityWeight_(0.0) , start_(NULL) , end_(NULL) , lower_(NULL) , upper_(NULL) , status_(NULL) , saveStatus_(NULL) , savedKeyVariable_(NULL) , backward_(NULL) , backToPivotRow_(NULL) , changeCost_(NULL) , keyVariable_(NULL) , next_(NULL) , toIndex_(NULL) , fromIndex_(NULL) , model_(NULL) , numberDualInfeasibilities_(0) , numberPrimalInfeasibilities_(0) , noCheck_(-1) , numberSets_(0) , saveNumber_(0) , possiblePivotKey_(0) , gubSlackIn_(-1) , firstGub_(0) , lastGub_(0) , gubType_(0) { setType(16); } //------------------------------------------------------------------- // Destructor //------------------------------------------------------------------- ClpGubMatrix::~ClpGubMatrix() { delete[] start_; delete[] end_; delete[] lower_; delete[] upper_; delete[] status_; delete[] saveStatus_; delete[] savedKeyVariable_; delete[] backward_; delete[] backToPivotRow_; delete[] changeCost_; delete[] keyVariable_; delete[] next_; delete[] toIndex_; delete[] fromIndex_; } //---------------------------------------------------------------- // Assignment operator //------------------------------------------------------------------- ClpGubMatrix & ClpGubMatrix::operator=(const ClpGubMatrix &rhs) { if (this != &rhs) { ClpPackedMatrix::operator=(rhs); delete[] start_; delete[] end_; delete[] lower_; delete[] upper_; delete[] status_; delete[] saveStatus_; delete[] savedKeyVariable_; delete[] backward_; delete[] backToPivotRow_; delete[] changeCost_; delete[] keyVariable_; delete[] next_; delete[] toIndex_; delete[] fromIndex_; numberSets_ = rhs.numberSets_; saveNumber_ = rhs.saveNumber_; possiblePivotKey_ = rhs.possiblePivotKey_; gubSlackIn_ = rhs.gubSlackIn_; start_ = ClpCopyOfArray(rhs.start_, numberSets_); end_ = ClpCopyOfArray(rhs.end_, numberSets_); lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); status_ = ClpCopyOfArray(rhs.status_, numberSets_); saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); int numberColumns = getNumCols(); backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); // find longest set int *longest = new int[numberSets_]; CoinZeroN(longest, numberSets_); int j; for (j = 0; j < numberColumns; j++) { int iSet = backward_[j]; if (iSet >= 0) longest[iSet]++; } int length = 0; for (j = 0; j < numberSets_; j++) length = CoinMax(length, longest[j]); next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; infeasibilityWeight_ = rhs.infeasibilityWeight_; numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; noCheck_ = rhs.noCheck_; firstGub_ = rhs.firstGub_; lastGub_ = rhs.lastGub_; gubType_ = rhs.gubType_; model_ = rhs.model_; } return *this; } //------------------------------------------------------------------- // Clone //------------------------------------------------------------------- ClpMatrixBase *ClpGubMatrix::clone() const { return new ClpGubMatrix(*this); } /* Subset clone (without gaps). Duplicates are allowed and order is as given */ ClpMatrixBase * ClpGubMatrix::subsetClone(int numberRows, const int *whichRows, int numberColumns, const int *whichColumns) const { return new ClpGubMatrix(*this, numberRows, whichRows, numberColumns, whichColumns); } /* Returns a new matrix in reverse order without gaps Is allowed to return NULL if doesn't want to have row copy */ ClpMatrixBase * ClpGubMatrix::reverseOrderedCopy() const { return NULL; } int ClpGubMatrix::hiddenRows() const { return numberSets_; } /* Subset constructor (without gaps). Duplicates are allowed and order is as given */ ClpGubMatrix::ClpGubMatrix( const ClpGubMatrix &rhs, int numberRows, const int *whichRows, int numberColumns, const int *whichColumns) : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns) { // Assuming no gub rows deleted // We also assume all sets in same order // Get array with backward pointers int numberColumnsOld = rhs.matrix_->getNumCols(); int *array = new int[numberColumnsOld]; int i; for (i = 0; i < numberColumnsOld; i++) array[i] = -1; for (int iSet = 0; iSet < numberSets_; iSet++) { for (int j = start_[iSet]; j < end_[iSet]; j++) array[j] = iSet; } numberSets_ = -1; int lastSet = -1; bool inSet = false; for (i = 0; i < numberColumns; i++) { int iColumn = whichColumns[i]; int iSet = array[iColumn]; if (iSet < 0) { inSet = false; } else { if (!inSet) { // start of new set but check okay if (iSet <= lastSet) throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); lastSet = iSet; numberSets_++; start_[numberSets_] = i; end_[numberSets_] = i + 1; lower_[numberSets_] = lower_[iSet]; upper_[numberSets_] = upper_[iSet]; inSet = true; } else { if (iSet < lastSet) { throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); } else if (iSet == lastSet) { end_[numberSets_] = i + 1; } else { // new set lastSet = iSet; numberSets_++; start_[numberSets_] = i; end_[numberSets_] = i + 1; lower_[numberSets_] = lower_[iSet]; upper_[numberSets_] = upper_[iSet]; } } } } delete[] array; numberSets_++; // adjust // Find type of gub firstGub_ = numberColumns + 1; lastGub_ = -1; for (i = 0; i < numberColumns; i++) { if (backward_[i] >= 0) { firstGub_ = CoinMin(firstGub_, i); lastGub_ = CoinMax(lastGub_, i); } } if (lastGub_ > 0) lastGub_++; gubType_ = 0; for (i = firstGub_; i < lastGub_; i++) { if (backward_[i] < 0) { gubType_ = 1; break; } } // Make sure key is feasible if only key in set } ClpGubMatrix::ClpGubMatrix( const CoinPackedMatrix &rhs, int numberRows, const int *whichRows, int numberColumns, const int *whichColumns) : ClpPackedMatrix(rhs, numberRows, whichRows, numberColumns, whichColumns) , sumDualInfeasibilities_(0.0) , sumPrimalInfeasibilities_(0.0) , sumOfRelaxedDualInfeasibilities_(0.0) , sumOfRelaxedPrimalInfeasibilities_(0.0) , start_(NULL) , end_(NULL) , lower_(NULL) , upper_(NULL) , backward_(NULL) , backToPivotRow_(NULL) , changeCost_(NULL) , keyVariable_(NULL) , next_(NULL) , toIndex_(NULL) , fromIndex_(NULL) , numberDualInfeasibilities_(0) , numberPrimalInfeasibilities_(0) , numberSets_(0) , saveNumber_(0) , possiblePivotKey_(0) , gubSlackIn_(-1) , firstGub_(0) , lastGub_(0) , gubType_(0) { setType(16); } /* Return x * A + y in z. Squashes small elements and knows about ClpSimplex */ void ClpGubMatrix::transposeTimes(const ClpSimplex *model, double scalar, const CoinIndexedVector *rowArray, CoinIndexedVector *y, CoinIndexedVector *columnArray) const { columnArray->clear(); double *pi = rowArray->denseVector(); int numberNonZero = 0; int *index = columnArray->getIndices(); double *array = columnArray->denseVector(); int numberInRowArray = rowArray->getNumElements(); // maybe I need one in OsiSimplex double zeroTolerance = model->zeroTolerance(); int numberRows = model->numberRows(); ClpPackedMatrix *rowCopy = dynamic_cast< ClpPackedMatrix * >(model->rowCopy()); bool packed = rowArray->packedMode(); double factor = 0.3; // We may not want to do by row if there may be cache problems int numberColumns = model->numberColumns(); // It would be nice to find L2 cache size - for moment 512K // Be slightly optimistic if (numberColumns * sizeof(double) > 1000000) { if (numberRows * 10 < numberColumns) factor = 0.1; else if (numberRows * 4 < numberColumns) factor = 0.15; else if (numberRows * 2 < numberColumns) factor = 0.2; //if (model->numberIterations()%50==0) //printf("%d nonzero\n",numberInRowArray); } // reduce for gub factor *= 0.5; assert(!y->getNumElements()); if (numberInRowArray > factor * numberRows || !rowCopy) { // do by column int iColumn; // get matrix data pointers const int *row = matrix_->getIndices(); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *columnLength = matrix_->getVectorLengths(); const double *elementByColumn = matrix_->getElements(); const double *rowScale = model->rowScale(); int numberColumns = model->numberColumns(); int iSet = -1; double djMod = 0.0; if (packed) { // need to expand pi into y assert(y->capacity() >= numberRows); double *piOld = pi; pi = y->denseVector(); const int *whichRow = rowArray->getIndices(); int i; if (!rowScale) { // modify pi so can collapse to one loop for (i = 0; i < numberInRowArray; i++) { int iRow = whichRow[i]; pi[iRow] = scalar * piOld[i]; } for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (backward_[iColumn] != iSet) { // get pi on gub row iSet = backward_[iColumn]; if (iSet >= 0) { int iBasic = keyVariable_[iSet]; if (iBasic < numberColumns) { // get dj without assert(model->getStatus(iBasic) == ClpSimplex::basic); djMod = 0.0; for (CoinBigIndex j = columnStart[iBasic]; j < columnStart[iBasic] + columnLength[iBasic]; j++) { int jRow = row[j]; djMod -= pi[jRow] * elementByColumn[j]; } } else { djMod = 0.0; } } else { djMod = 0.0; } } double value = -djMod; CoinBigIndex j; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } if (fabs(value) > zeroTolerance) { array[numberNonZero] = value; index[numberNonZero++] = iColumn; } } } else { // scaled // modify pi so can collapse to one loop for (i = 0; i < numberInRowArray; i++) { int iRow = whichRow[i]; pi[iRow] = scalar * piOld[i] * rowScale[iRow]; } for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (backward_[iColumn] != iSet) { // get pi on gub row iSet = backward_[iColumn]; if (iSet >= 0) { int iBasic = keyVariable_[iSet]; if (iBasic < numberColumns) { // get dj without assert(model->getStatus(iBasic) == ClpSimplex::basic); djMod = 0.0; // scaled for (CoinBigIndex j = columnStart[iBasic]; j < columnStart[iBasic] + columnLength[iBasic]; j++) { int jRow = row[j]; djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow]; } } else { djMod = 0.0; } } else { djMod = 0.0; } } double value = -djMod; CoinBigIndex j; const double *columnScale = model->columnScale(); for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } value *= columnScale[iColumn]; if (fabs(value) > zeroTolerance) { array[numberNonZero] = value; index[numberNonZero++] = iColumn; } } } // zero out for (i = 0; i < numberInRowArray; i++) { int iRow = whichRow[i]; pi[iRow] = 0.0; } } else { // code later assert(packed); if (!rowScale) { if (scalar == -1.0) { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = -value; } } } else if (scalar == 1.0) { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = value; } } } else { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } value *= scalar; if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = value; } } } } else { // scaled if (scalar == -1.0) { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; const double *columnScale = model->columnScale(); for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; } value *= columnScale[iColumn]; if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = -value; } } } else if (scalar == 1.0) { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; const double *columnScale = model->columnScale(); for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; } value *= columnScale[iColumn]; if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = value; } } } else { for (iColumn = 0; iColumn < numberColumns; iColumn++) { double value = 0.0; CoinBigIndex j; const double *columnScale = model->columnScale(); for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; } value *= scalar * columnScale[iColumn]; if (fabs(value) > zeroTolerance) { index[numberNonZero++] = iColumn; array[iColumn] = value; } } } } } columnArray->setNumElements(numberNonZero); y->setNumElements(0); } else { // do by row transposeTimesByRow(model, scalar, rowArray, y, columnArray); } if (packed) columnArray->setPackedMode(true); if (0) { columnArray->checkClean(); int numberNonZero = columnArray->getNumElements(); ; int *index = columnArray->getIndices(); double *array = columnArray->denseVector(); int i; for (i = 0; i < numberNonZero; i++) { int j = index[i]; double value; if (packed) value = array[i]; else value = array[j]; printf("Ti %d %d %g\n", i, j, value); } } } /* Return x * A + y in z. Squashes small elements and knows about ClpSimplex */ void ClpGubMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar, const CoinIndexedVector *rowArray, CoinIndexedVector *y, CoinIndexedVector *columnArray) const { // Do packed part ClpPackedMatrix::transposeTimesByRow(model, scalar, rowArray, y, columnArray); if (numberSets_) { /* what we need to do is do by row as normal but get list of sets touched and then update those ones */ abort(); } } /* Return x *A in z but just for indices in y. */ void ClpGubMatrix::subsetTransposeTimes(const ClpSimplex *model, const CoinIndexedVector *rowArray, const CoinIndexedVector *y, CoinIndexedVector *columnArray) const { columnArray->clear(); double *pi = rowArray->denseVector(); double *array = columnArray->denseVector(); int jColumn; // get matrix data pointers const int *row = matrix_->getIndices(); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *columnLength = matrix_->getVectorLengths(); const double *elementByColumn = matrix_->getElements(); const double *rowScale = model->rowScale(); int numberToDo = y->getNumElements(); const int *which = y->getIndices(); assert(!rowArray->packedMode()); columnArray->setPacked(); int numberTouched = 0; if (!rowScale) { for (jColumn = 0; jColumn < numberToDo; jColumn++) { int iColumn = which[jColumn]; double value = 0.0; CoinBigIndex j; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j]; } array[jColumn] = value; if (value) { int iSet = backward_[iColumn]; if (iSet >= 0) { int iBasic = keyVariable_[iSet]; if (iBasic == iColumn) { toIndex_[iSet] = jColumn; fromIndex_[numberTouched++] = iSet; } } } } } else { // scaled for (jColumn = 0; jColumn < numberToDo; jColumn++) { int iColumn = which[jColumn]; double value = 0.0; CoinBigIndex j; const double *columnScale = model->columnScale(); for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; } value *= columnScale[iColumn]; array[jColumn] = value; if (value) { int iSet = backward_[iColumn]; if (iSet >= 0) { int iBasic = keyVariable_[iSet]; if (iBasic == iColumn) { toIndex_[iSet] = jColumn; fromIndex_[numberTouched++] = iSet; } } } } } // adjust djs for (jColumn = 0; jColumn < numberToDo; jColumn++) { int iColumn = which[jColumn]; int iSet = backward_[iColumn]; if (iSet >= 0) { int kColumn = toIndex_[iSet]; if (kColumn >= 0) array[jColumn] -= array[kColumn]; } } // and clear basic for (int j = 0; j < numberTouched; j++) { int iSet = fromIndex_[j]; int kColumn = toIndex_[iSet]; toIndex_[iSet] = -1; array[kColumn] = 0.0; } } /// returns number of elements in column part of basis, int ClpGubMatrix::countBasis(const int *whichColumn, int &numberColumnBasic) { int i; int numberColumns = getNumCols(); const int *columnLength = matrix_->getVectorLengths(); int numberRows = getNumRows(); int numberBasic = 0; CoinBigIndex numberElements = 0; int lastSet = -1; int key = -1; int keyLength = -1; double *work = new double[numberRows]; CoinZeroN(work, numberRows); char *mark = new char[numberRows]; CoinZeroN(mark, numberRows); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *row = matrix_->getIndices(); const double *elementByColumn = matrix_->getElements(); //ClpGubDynamicMatrix* gubx = //dynamic_cast< ClpGubDynamicMatrix*>(this); //int * id = gubx->id(); // just count for (i = 0; i < numberColumnBasic; i++) { int iColumn = whichColumn[i]; int iSet = backward_[iColumn]; int length = columnLength[iColumn]; if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { numberElements += length; numberBasic++; //printf("non gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,length); } else { // in gub set if (iColumn != keyVariable_[iSet]) { numberBasic++; CoinBigIndex j; // not key if (lastSet < iSet) { // erase work if (key >= 0) { for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) work[row[j]] = 0.0; } key = keyVariable_[iSet]; lastSet = iSet; keyLength = columnLength[key]; for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) work[row[j]] = elementByColumn[j]; } int extra = keyLength; for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { int iRow = row[j]; double keyValue = work[iRow]; double value = elementByColumn[j]; if (!keyValue) { if (fabs(value) > 1.0e-20) extra++; } else { value -= keyValue; if (fabs(value) <= 1.0e-20) extra--; } } numberElements += extra; //printf("gub - set %d id %d (column %d) nel %d\n",iSet,id[iColumn-20],iColumn,extra); } } } delete[] work; delete[] mark; // update number of column basic numberColumnBasic = numberBasic; if (numberElements > COIN_INT_MAX) { printf("Factorization too large\n"); abort(); } return static_cast< int >(numberElements); } void ClpGubMatrix::fillBasis(ClpSimplex *model, const int *whichColumn, int &numberColumnBasic, int *indexRowU, int *start, int *rowCount, int *columnCount, CoinFactorizationDouble *elementU) { int i; int numberColumns = getNumCols(); const int *columnLength = matrix_->getVectorLengths(); int numberRows = getNumRows(); assert(next_ || !elementU); CoinBigIndex numberElements = start[0]; int lastSet = -1; int key = -1; int keyLength = -1; double *work = new double[numberRows]; CoinZeroN(work, numberRows); char *mark = new char[numberRows]; CoinZeroN(mark, numberRows); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *row = matrix_->getIndices(); const double *elementByColumn = matrix_->getElements(); const double *rowScale = model->rowScale(); int numberBasic = 0; if (0) { printf("%d basiccolumns\n", numberColumnBasic); int i; for (i = 0; i < numberSets_; i++) { int k = keyVariable_[i]; if (k < numberColumns) { printf("key %d on set %d, %d elements\n", k, i, columnStart[k + 1] - columnStart[k]); for (int j = columnStart[k]; j < columnStart[k + 1]; j++) printf("row %d el %g\n", row[j], elementByColumn[j]); } else { printf("slack key on set %d\n", i); } } } // fill if (!rowScale) { // no scaling for (i = 0; i < numberColumnBasic; i++) { int iColumn = whichColumn[i]; int iSet = backward_[iColumn]; int length = columnLength[iColumn]; if (0) { int k = iColumn; printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k + 1] - columnStart[k]); for (int j = columnStart[k]; j < columnStart[k + 1]; j++) printf("row %d el %g\n", row[j], elementByColumn[j]); } CoinBigIndex j; if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { double value = elementByColumn[j]; if (fabs(value) > 1.0e-20) { int iRow = row[j]; indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } // end of column columnCount[numberBasic] = numberElements - start[numberBasic]; numberBasic++; start[numberBasic] = numberElements; } else { // in gub set if (iColumn != keyVariable_[iSet]) { // not key if (lastSet != iSet) { // erase work if (key >= 0) { for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; work[iRow] = 0.0; mark[iRow] = 0; } } key = keyVariable_[iSet]; lastSet = iSet; keyLength = columnLength[key]; for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; work[iRow] = elementByColumn[j]; mark[iRow] = 1; } } for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { int iRow = row[j]; double value = elementByColumn[j]; if (mark[iRow]) { mark[iRow] = 0; double keyValue = work[iRow]; value -= keyValue; } if (fabs(value) > 1.0e-20) { indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; if (mark[iRow]) { double value = -work[iRow]; if (fabs(value) > 1.0e-20) { indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } else { // just put back mark mark[iRow] = 1; } } // end of column columnCount[numberBasic] = numberElements - start[numberBasic]; numberBasic++; start[numberBasic] = numberElements; } } } } else { // scaling const double *columnScale = model->columnScale(); for (i = 0; i < numberColumnBasic; i++) { int iColumn = whichColumn[i]; int iSet = backward_[iColumn]; int length = columnLength[iColumn]; CoinBigIndex j; if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { double scale = columnScale[iColumn]; for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { int iRow = row[j]; double value = elementByColumn[j] * scale * rowScale[iRow]; if (fabs(value) > 1.0e-20) { indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } // end of column columnCount[numberBasic] = numberElements - start[numberBasic]; numberBasic++; start[numberBasic] = numberElements; } else { // in gub set if (iColumn != keyVariable_[iSet]) { double scale = columnScale[iColumn]; // not key if (lastSet < iSet) { // erase work if (key >= 0) { for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; work[iRow] = 0.0; mark[iRow] = 0; } } key = keyVariable_[iSet]; lastSet = iSet; keyLength = columnLength[key]; double scale = columnScale[key]; for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; work[iRow] = elementByColumn[j] * scale * rowScale[iRow]; mark[iRow] = 1; } } for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { int iRow = row[j]; double value = elementByColumn[j] * scale * rowScale[iRow]; if (mark[iRow]) { mark[iRow] = 0; double keyValue = work[iRow]; value -= keyValue; } if (fabs(value) > 1.0e-20) { indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { int iRow = row[j]; if (mark[iRow]) { double value = -work[iRow]; if (fabs(value) > 1.0e-20) { indexRowU[numberElements] = iRow; rowCount[iRow]++; elementU[numberElements++] = value; } } else { // just put back mark mark[iRow] = 1; } } // end of column columnCount[numberBasic] = numberElements - start[numberBasic]; numberBasic++; start[numberBasic] = numberElements; } } } } delete[] work; delete[] mark; // update number of column basic numberColumnBasic = numberBasic; } /* Unpacks a column into an CoinIndexedvector */ void ClpGubMatrix::unpack(const ClpSimplex *model, CoinIndexedVector *rowArray, int iColumn) const { assert(iColumn < model->numberColumns()); // Do packed part ClpPackedMatrix::unpack(model, rowArray, iColumn); int iSet = backward_[iColumn]; if (iSet >= 0) { int iBasic = keyVariable_[iSet]; if (iBasic < model->numberColumns()) { add(model, rowArray, iBasic, -1.0); } } } /* Unpacks a column into a CoinIndexedVector ** in packed format Note that model is NOT const. Bounds and objective could be modified if doing column generation (just for this variable) */ void ClpGubMatrix::unpackPacked(ClpSimplex *model, CoinIndexedVector *rowArray, int iColumn) const { int numberColumns = model->numberColumns(); if (iColumn < numberColumns) { // Do packed part ClpPackedMatrix::unpackPacked(model, rowArray, iColumn); int iSet = backward_[iColumn]; if (iSet >= 0) { // columns are in order int iBasic = keyVariable_[iSet]; if (iBasic < numberColumns) { int number = rowArray->getNumElements(); const double *rowScale = model->rowScale(); const int *row = matrix_->getIndices(); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *columnLength = matrix_->getVectorLengths(); const double *elementByColumn = matrix_->getElements(); double *array = rowArray->denseVector(); int *index = rowArray->getIndices(); CoinBigIndex i; int numberOld = number; int lastIndex = 0; int next = index[lastIndex]; if (!rowScale) { for (i = columnStart[iBasic]; i < columnStart[iBasic] + columnLength[iBasic]; i++) { int iRow = row[i]; while (iRow > next) { lastIndex++; if (lastIndex == numberOld) next = matrix_->getNumRows(); else next = index[lastIndex]; } if (iRow < next) { array[number] = -elementByColumn[i]; index[number++] = iRow; } else { assert(iRow == next); array[lastIndex] -= elementByColumn[i]; if (!array[lastIndex]) array[lastIndex] = 1.0e-100; } } } else { // apply scaling double scale = model->columnScale()[iBasic]; for (i = columnStart[iBasic]; i < columnStart[iBasic] + columnLength[iBasic]; i++) { int iRow = row[i]; while (iRow > next) { lastIndex++; if (lastIndex == numberOld) next = matrix_->getNumRows(); else next = index[lastIndex]; } if (iRow < next) { array[number] = -elementByColumn[i] * scale * rowScale[iRow]; index[number++] = iRow; } else { assert(iRow == next); array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow]; if (!array[lastIndex]) array[lastIndex] = 1.0e-100; } } } rowArray->setNumElements(number); } } } else { // key slack entering int iBasic = keyVariable_[gubSlackIn_]; assert(iBasic < numberColumns); int number = 0; const double *rowScale = model->rowScale(); const int *row = matrix_->getIndices(); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *columnLength = matrix_->getVectorLengths(); const double *elementByColumn = matrix_->getElements(); double *array = rowArray->denseVector(); int *index = rowArray->getIndices(); CoinBigIndex i; if (!rowScale) { for (i = columnStart[iBasic]; i < columnStart[iBasic] + columnLength[iBasic]; i++) { int iRow = row[i]; array[number] = elementByColumn[i]; index[number++] = iRow; } } else { // apply scaling double scale = model->columnScale()[iBasic]; for (i = columnStart[iBasic]; i < columnStart[iBasic] + columnLength[iBasic]; i++) { int iRow = row[i]; array[number] = elementByColumn[i] * scale * rowScale[iRow]; index[number++] = iRow; } } rowArray->setNumElements(number); rowArray->setPacked(); } } /* Adds multiple of a column into an CoinIndexedvector You can use quickAdd to add to vector */ void ClpGubMatrix::add(const ClpSimplex *model, CoinIndexedVector *rowArray, int iColumn, double multiplier) const { assert(iColumn < model->numberColumns()); // Do packed part ClpPackedMatrix::add(model, rowArray, iColumn, multiplier); int iSet = backward_[iColumn]; if (iSet >= 0 && iColumn != keyVariable_[iSet]) { ClpPackedMatrix::add(model, rowArray, keyVariable_[iSet], -multiplier); } } /* Adds multiple of a column into an array */ void ClpGubMatrix::add(const ClpSimplex *model, double *array, int iColumn, double multiplier) const { assert(iColumn < model->numberColumns()); // Do packed part ClpPackedMatrix::add(model, array, iColumn, multiplier); if (iColumn < model->numberColumns()) { int iSet = backward_[iColumn]; if (iSet >= 0 && iColumn != keyVariable_[iSet] && keyVariable_[iSet] < model->numberColumns()) { ClpPackedMatrix::add(model, array, keyVariable_[iSet], -multiplier); } } } // Partial pricing void ClpGubMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction, int &bestSequence, int &numberWanted) { numberWanted = currentWanted_; if (numberSets_) { // Do packed part before gub int numberColumns = matrix_->getNumCols(); double ratio = static_cast< double >(firstGub_) / static_cast< double >(numberColumns); ClpPackedMatrix::partialPricing(model, startFraction * ratio, endFraction * ratio, bestSequence, numberWanted); if (numberWanted || minimumGoodReducedCosts_ < -1) { // do gub const double *element = matrix_->getElements(); const int *row = matrix_->getIndices(); const CoinBigIndex *startColumn = matrix_->getVectorStarts(); const int *length = matrix_->getVectorLengths(); const double *rowScale = model->rowScale(); const double *columnScale = model->columnScale(); int iSequence; CoinBigIndex j; double tolerance = model->currentDualTolerance(); double *reducedCost = model->djRegion(); const double *duals = model->dualRowSolution(); const double *cost = model->costRegion(); double bestDj; int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); if (bestSequence >= 0) bestDj = fabs(this->reducedCost(model, bestSequence)); else bestDj = tolerance; int sequenceOut = model->sequenceOut(); int saveSequence = bestSequence; int startG = firstGub_ + static_cast< int >(startFraction * (lastGub_ - firstGub_)); int endG = firstGub_ + static_cast< int >(endFraction * (lastGub_ - firstGub_)); endG = CoinMin(lastGub_, endG + 1); // If nothing found yet can go all the way to end int endAll = endG; if (bestSequence < 0 && !startG) endAll = lastGub_; int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_; int nSets = 0; int iSet = -1; double djMod = 0.0; double infeasibilityCost = model->infeasibilityCost(); if (rowScale) { double bestDjMod = 0.0; // scaled for (iSequence = startG; iSequence < endAll; iSequence++) { if (numberWanted + minNeg < originalWanted_ && nSets > minSet) { // give up numberWanted = 0; break; } else if (iSequence == endG && bestSequence >= 0) { break; } if (backward_[iSequence] != iSet) { // get pi on gub row iSet = backward_[iSequence]; if (iSet >= 0) { nSets++; int iBasic = keyVariable_[iSet]; if (iBasic >= numberColumns) { djMod = -weight(iSet) * infeasibilityCost; } else { // get dj without assert(model->getStatus(iBasic) == ClpSimplex::basic); djMod = 0.0; // scaled for (j = startColumn[iBasic]; j < startColumn[iBasic] + length[iBasic]; j++) { int jRow = row[j]; djMod -= duals[jRow] * element[j] * rowScale[jRow]; } // allow for scaling djMod += cost[iBasic] / columnScale[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 = numberRows + numberColumns + iSet; bestDjMod = djMod; } 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 = numberRows + numberColumns + iSet; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; abort(); } } } } } } else { // not in set djMod = 0.0; } } if (iSequence != sequenceOut) { double value; ClpSimplex::Status status = model->getStatus(iSequence); switch (status) { case ClpSimplex::basic: case ClpSimplex::isFixed: break; case ClpSimplex::isFree: case ClpSimplex::superBasic: value = -djMod; // scaled for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j] * rowScale[jRow]; } value = fabs(cost[iSequence] + value * columnScale[iSequence]); if (value > FREE_ACCEPT * tolerance) { numberWanted--; // we are going to bias towards free (but only if reasonable) value *= FREE_BIAS; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; case ClpSimplex::atUpperBound: value = -djMod; // scaled for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j] * rowScale[jRow]; } value = cost[iSequence] + value * columnScale[iSequence]; if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; case ClpSimplex::atLowerBound: value = -djMod; // scaled for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j] * rowScale[jRow]; } value = -(cost[iSequence] + value * columnScale[iSequence]); if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; } } if (!numberWanted) break; } if (bestSequence != saveSequence) { if (bestSequence < numberRows + numberColumns) { // recompute dj double value = bestDjMod; // scaled for (j = startColumn[bestSequence]; j < startColumn[bestSequence] + length[bestSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j] * rowScale[jRow]; } reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence]; gubSlackIn_ = -1; } else { // slack - make last column gubSlackIn_ = bestSequence - numberRows - numberColumns; bestSequence = numberColumns + 2 * numberRows; reducedCost[bestSequence] = bestDjMod; model->setStatus(bestSequence, getStatus(gubSlackIn_)); if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound) model->solutionRegion()[bestSequence] = upper_[gubSlackIn_]; else model->solutionRegion()[bestSequence] = lower_[gubSlackIn_]; model->lowerRegion()[bestSequence] = lower_[gubSlackIn_]; model->upperRegion()[bestSequence] = upper_[gubSlackIn_]; model->costRegion()[bestSequence] = 0.0; } savedBestSequence_ = bestSequence; savedBestDj_ = reducedCost[savedBestSequence_]; } } else { double bestDjMod = 0.0; //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(), // startG,endG,numberWanted); for (iSequence = startG; iSequence < endG; iSequence++) { if (numberWanted + minNeg < originalWanted_ && nSets > minSet) { // give up numberWanted = 0; break; } else if (iSequence == endG && bestSequence >= 0) { break; } if (backward_[iSequence] != iSet) { // get pi on gub row iSet = backward_[iSequence]; if (iSet >= 0) { nSets++; int iBasic = keyVariable_[iSet]; if (iBasic >= numberColumns) { djMod = -weight(iSet) * infeasibilityCost; } else { // get dj without assert(model->getStatus(iBasic) == ClpSimplex::basic); djMod = 0.0; for (j = startColumn[iBasic]; j < startColumn[iBasic] + length[iBasic]; 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 = numberRows + numberColumns + iSet; bestDjMod = djMod; } 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 = numberRows + numberColumns + iSet; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; abort(); } } } } } } else { // not in set djMod = 0.0; } } if (iSequence != sequenceOut) { double value; ClpSimplex::Status status = model->getStatus(iSequence); switch (status) { case ClpSimplex::basic: case ClpSimplex::isFixed: break; case ClpSimplex::isFree: case ClpSimplex::superBasic: value = cost[iSequence] - djMod; for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j]; } value = fabs(value); if (value > FREE_ACCEPT * tolerance) { numberWanted--; // we are going to bias towards free (but only if reasonable) value *= FREE_BIAS; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; case ClpSimplex::atUpperBound: value = cost[iSequence] - djMod; for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j]; } if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; case ClpSimplex::atLowerBound: value = cost[iSequence] - djMod; for (j = startColumn[iSequence]; j < startColumn[iSequence] + length[iSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j]; } value = -value; if (value > tolerance) { numberWanted--; if (value > bestDj) { // check flagged variable and correct dj if (!model->flagged(iSequence)) { bestDj = value; bestSequence = iSequence; bestDjMod = djMod; } else { // just to make sure we don't exit before got something numberWanted++; } } } break; } } if (!numberWanted) break; } if (bestSequence != saveSequence) { if (bestSequence < numberRows + numberColumns) { // recompute dj double value = cost[bestSequence] - bestDjMod; for (j = startColumn[bestSequence]; j < startColumn[bestSequence] + length[bestSequence]; j++) { int jRow = row[j]; value -= duals[jRow] * element[j]; } //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod); reducedCost[bestSequence] = value; gubSlackIn_ = -1; } else { // slack - make last column gubSlackIn_ = bestSequence - numberRows - numberColumns; bestSequence = numberColumns + 2 * numberRows; reducedCost[bestSequence] = bestDjMod; //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod); model->setStatus(bestSequence, getStatus(gubSlackIn_)); if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound) model->solutionRegion()[bestSequence] = upper_[gubSlackIn_]; else model->solutionRegion()[bestSequence] = lower_[gubSlackIn_]; model->lowerRegion()[bestSequence] = lower_[gubSlackIn_]; model->upperRegion()[bestSequence] = upper_[gubSlackIn_]; model->costRegion()[bestSequence] = 0.0; } } } // See if may be finished if (startG == firstGub_ && bestSequence < 0) infeasibilityWeight_ = model_->infeasibilityCost(); else if (bestSequence >= 0) infeasibilityWeight_ = -1.0; } if (numberWanted) { // Do packed part after gub double offset = static_cast< double >(lastGub_) / static_cast< double >(numberColumns); double ratio = static_cast< double >(numberColumns) / static_cast< double >(numberColumns) - offset; double start2 = offset + ratio * startFraction; double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6); ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted); } } else { // no gub ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); } if (bestSequence >= 0) infeasibilityWeight_ = -1.0; // not optimal currentWanted_ = numberWanted; } /* expands an updated column to allow for extra rows which the main solver does not know about and returns number added. */ int ClpGubMatrix::extendUpdated(ClpSimplex *model, CoinIndexedVector *update, int mode) { // I think we only need to bother about sets with two in basis or incoming set int number = update->getNumElements(); double *array = update->denseVector(); int *index = update->getIndices(); int i; assert(!number || update->packedMode()); int *pivotVariable = model->pivotVariable(); int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); int numberTotal = numberRows + numberColumns; int sequenceIn = model->sequenceIn(); int returnCode = 0; int iSetIn; if (sequenceIn < numberColumns) { iSetIn = backward_[sequenceIn]; gubSlackIn_ = -1; // in case set } else if (sequenceIn < numberRows + numberColumns) { iSetIn = -1; gubSlackIn_ = -1; // in case set } else { iSetIn = gubSlackIn_; } double *lower = model->lowerRegion(); double *upper = model->upperRegion(); double *cost = model->costRegion(); double *solution = model->solutionRegion(); int number2 = number; if (!mode) { double primalTolerance = model->primalTolerance(); double infeasibilityCost = model->infeasibilityCost(); // extend saveNumber_ = number; for (i = 0; i < number; i++) { int iRow = index[i]; int iPivot = pivotVariable[iRow]; if (iPivot < numberColumns) { int iSet = backward_[iPivot]; if (iSet >= 0) { // two (or more) in set int iIndex = toIndex_[iSet]; double otherValue = array[i]; double value; if (iIndex < 0) { toIndex_[iSet] = number2; int iNew = number2 - number; fromIndex_[number2 - number] = iSet; iIndex = number2; index[number2] = numberRows + iNew; // do key stuff int iKey = keyVariable_[iSet]; if (iKey < numberColumns) { // Save current cost of key changeCost_[number2 - number] = cost[iKey]; if (iSet != iSetIn) value = 0.0; else if (iSetIn != gubSlackIn_) value = 1.0; else value = -1.0; pivotVariable[numberRows + iNew] = iKey; // Do I need to recompute? double sol; assert(getStatus(iSet) != ClpSimplex::basic); if (getStatus(iSet) == ClpSimplex::atLowerBound) sol = lower_[iSet]; else sol = upper_[iSet]; if ((gubType_ & 8) != 0) { int iColumn = next_[iKey]; // sum all non-key variables while (iColumn >= 0) { sol -= solution[iColumn]; iColumn = next_[iColumn]; } } else { int stop = -(iKey + 1); int iColumn = next_[iKey]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; sol -= solution[iColumn]; iColumn = next_[iColumn]; } } solution[iKey] = sol; if (model->algorithm() > 0) model->nonLinearCost()->setOne(iKey, sol); //assert (fabs(sol-solution[iKey])<1.0e-3); } else { // gub slack is basic // Save current cost of key changeCost_[number2 - number] = -weight(iSet) * infeasibilityCost; otherValue = -otherValue; //allow for - sign on slack if (iSet != iSetIn) value = 0.0; else value = -1.0; pivotVariable[numberRows + iNew] = iNew + numberTotal; model->djRegion()[iNew + numberTotal] = 0.0; double sol = 0.0; if ((gubType_ & 8) != 0) { int iColumn = next_[iKey]; // sum all non-key variables while (iColumn >= 0) { sol += solution[iColumn]; iColumn = next_[iColumn]; } } else { int stop = -(iKey + 1); int iColumn = next_[iKey]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; sol += solution[iColumn]; iColumn = next_[iColumn]; } } solution[iNew + numberTotal] = sol; // and do cost in nonLinearCost if (model->algorithm() > 0) model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]); if (sol > upper_[iSet] + primalTolerance) { setAbove(iSet); lower[iNew + numberTotal] = upper_[iSet]; upper[iNew + numberTotal] = COIN_DBL_MAX; } else if (sol < lower_[iSet] - primalTolerance) { setBelow(iSet); lower[iNew + numberTotal] = -COIN_DBL_MAX; upper[iNew + numberTotal] = lower_[iSet]; } else { setFeasible(iSet); lower[iNew + numberTotal] = lower_[iSet]; upper[iNew + numberTotal] = upper_[iSet]; } cost[iNew + numberTotal] = weight(iSet) * infeasibilityCost; } number2++; } else { value = array[iIndex]; int iKey = keyVariable_[iSet]; if (iKey >= numberColumns) otherValue = -otherValue; //allow for - sign on slack } value -= otherValue; array[iIndex] = value; } } } if (iSetIn >= 0 && toIndex_[iSetIn] < 0) { // Do incoming update->setPacked(); // just in case no elements toIndex_[iSetIn] = number2; int iNew = number2 - number; fromIndex_[number2 - number] = iSetIn; // Save current cost of key double currentCost; int key = keyVariable_[iSetIn]; if (key < numberColumns) currentCost = cost[key]; else currentCost = -weight(iSetIn) * infeasibilityCost; changeCost_[number2 - number] = currentCost; index[number2] = numberRows + iNew; // do key stuff int iKey = keyVariable_[iSetIn]; if (iKey < numberColumns) { if (gubSlackIn_ < 0) array[number2] = 1.0; else array[number2] = -1.0; pivotVariable[numberRows + iNew] = iKey; // Do I need to recompute? double sol; assert(getStatus(iSetIn) != ClpSimplex::basic); if (getStatus(iSetIn) == ClpSimplex::atLowerBound) sol = lower_[iSetIn]; else sol = upper_[iSetIn]; if ((gubType_ & 8) != 0) { int iColumn = next_[iKey]; // sum all non-key variables while (iColumn >= 0) { sol -= solution[iColumn]; iColumn = next_[iColumn]; } } else { // bounds exist - sum over all except key int stop = -(iKey + 1); int iColumn = next_[iKey]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; sol -= solution[iColumn]; iColumn = next_[iColumn]; } } solution[iKey] = sol; if (model->algorithm() > 0) model->nonLinearCost()->setOne(iKey, sol); //assert (fabs(sol-solution[iKey])<1.0e-3); } else { // gub slack is basic array[number2] = -1.0; pivotVariable[numberRows + iNew] = iNew + numberTotal; model->djRegion()[iNew + numberTotal] = 0.0; double sol = 0.0; if ((gubType_ & 8) != 0) { int iColumn = next_[iKey]; // sum all non-key variables while (iColumn >= 0) { sol += solution[iColumn]; iColumn = next_[iColumn]; } } else { // bounds exist - sum over all except key int stop = -(iKey + 1); int iColumn = next_[iKey]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; sol += solution[iColumn]; iColumn = next_[iColumn]; } } solution[iNew + numberTotal] = sol; // and do cost in nonLinearCost if (model->algorithm() > 0) model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]); if (sol > upper_[iSetIn] + primalTolerance) { setAbove(iSetIn); lower[iNew + numberTotal] = upper_[iSetIn]; upper[iNew + numberTotal] = COIN_DBL_MAX; } else if (sol < lower_[iSetIn] - primalTolerance) { setBelow(iSetIn); lower[iNew + numberTotal] = -COIN_DBL_MAX; upper[iNew + numberTotal] = lower_[iSetIn]; } else { setFeasible(iSetIn); lower[iNew + numberTotal] = lower_[iSetIn]; upper[iNew + numberTotal] = upper_[iSetIn]; } cost[iNew + numberTotal] = weight(iSetIn) * infeasibilityCost; } number2++; } // mark end fromIndex_[number2 - number] = -1; returnCode = number2 - number; // make sure lower_ upper_ adjusted synchronize(model, 9); } else { // take off? if (number > saveNumber_) { // clear double theta = model->theta(); double *solution = model->solutionRegion(); for (i = saveNumber_; i < number; i++) { int iRow = index[i]; int iColumn = pivotVariable[iRow]; #ifdef CLP_DEBUG_PRINT printf("Column %d (set %d) lower %g, upper %g - alpha %g - old value %g, new %g (theta %g)\n", iColumn, fromIndex_[i - saveNumber_], lower[iColumn], upper[iColumn], array[i], solution[iColumn], solution[iColumn] - model->theta() * array[i], model->theta()); #endif double value = array[i]; array[i] = 0.0; int iSet = fromIndex_[i - saveNumber_]; toIndex_[iSet] = -1; if (iSet == iSetIn && iColumn < numberColumns) { // update as may need value solution[iColumn] -= theta * value; } } } #ifdef CLP_DEBUG for (i = 0; i < numberSets_; i++) assert(toIndex_[i] == -1); #endif number2 = saveNumber_; } update->setNumElements(number2); return returnCode; } /* utility primal function for dealing with dynamic constraints mode=n see ClpGubMatrix.hpp for definition Remember to update here when settled down */ void ClpGubMatrix::primalExpanded(ClpSimplex *model, int mode) { int numberColumns = model->numberColumns(); switch (mode) { // If key variable then slot in gub rhs so will get correct contribution case 0: { int i; double *solution = model->solutionRegion(); ClpSimplex::Status iStatus; for (i = 0; i < numberSets_; i++) { int iColumn = keyVariable_[i]; if (iColumn < numberColumns) { // key is structural - where is slack iStatus = getStatus(i); assert(iStatus != ClpSimplex::basic); if (iStatus == ClpSimplex::atLowerBound) solution[iColumn] = lower_[i]; else solution[iColumn] = upper_[i]; } } } break; // Compute values of key variables case 1: { int i; double *solution = model->solutionRegion(); //const int * columnLength = matrix_->getVectorLengths(); //const CoinBigIndex * columnStart = matrix_->getVectorStarts(); //const int * row = matrix_->getIndices(); //const double * elementByColumn = matrix_->getElements(); //int * pivotVariable = model->pivotVariable(); sumPrimalInfeasibilities_ = 0.0; numberPrimalInfeasibilities_ = 0; double primalTolerance = model->primalTolerance(); double relaxedTolerance = primalTolerance; // we can't really trust infeasibilities if there is primal error double error = CoinMin(1.0e-2, model->largestPrimalError()); // allow tolerance at least slightly bigger than standard relaxedTolerance = relaxedTolerance + error; // but we will be using difference relaxedTolerance -= primalTolerance; sumOfRelaxedPrimalInfeasibilities_ = 0.0; for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds) int kColumn = keyVariable_[i]; double value = 0.0; if ((gubType_ & 8) != 0) { int iColumn = next_[kColumn]; // sum all non-key variables while (iColumn >= 0) { value += solution[iColumn]; iColumn = next_[iColumn]; } } else { // bounds exist - sum over all except key int stop = -(kColumn + 1); int iColumn = next_[kColumn]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; value += solution[iColumn]; iColumn = next_[iColumn]; } } if (kColumn < numberColumns) { // make sure key is basic - so will be skipped in values pass model->setStatus(kColumn, ClpSimplex::basic); // feasibility will be done later assert(getStatus(i) != ClpSimplex::basic); if (getStatus(i) == ClpSimplex::atUpperBound) solution[kColumn] = upper_[i] - value; else solution[kColumn] = lower_[i] - value; //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]); } else { // slack is key assert(getStatus(i) == ClpSimplex::basic); double infeasibility = 0.0; if (value > upper_[i] + primalTolerance) { infeasibility = value - upper_[i] - primalTolerance; setAbove(i); } else if (value < lower_[i] - primalTolerance) { infeasibility = lower_[i] - value - primalTolerance; setBelow(i); } else { setFeasible(i); } //printf("Value of key slack for set %d is %g\n",i,value); if (infeasibility > 0.0) { sumPrimalInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedPrimalInfeasibilities_ += infeasibility; numberPrimalInfeasibilities_++; } } } } break; // Report on infeasibilities of key variables case 2: { model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + sumPrimalInfeasibilities_); model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + numberPrimalInfeasibilities_); model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + sumOfRelaxedPrimalInfeasibilities_); } break; } } /* utility dual function for dealing with dynamic constraints mode=n see ClpGubMatrix.hpp for definition Remember to update here when settled down */ void ClpGubMatrix::dualExpanded(ClpSimplex *model, CoinIndexedVector *array, double * /*other*/, int mode) { switch (mode) { // modify costs before transposeUpdate case 0: { int i; double *cost = model->costRegion(); // not dual values yet //assert (!other); //double * work = array->denseVector(); double infeasibilityCost = model->infeasibilityCost(); int *pivotVariable = model->pivotVariable(); int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns) { int iSet = backward_[iPivot]; if (iSet >= 0) { int kColumn = keyVariable_[iSet]; double costValue; if (kColumn < numberColumns) { // structural has cost costValue = cost[kColumn]; } else { // slack is key assert(getStatus(iSet) == ClpSimplex::basic); // negative as -1.0 for slack costValue = -weight(iSet) * infeasibilityCost; } array->add(i, -costValue); // was work[i]-costValue } } } } break; // create duals for key variables (without check on dual infeasible) case 1: { // If key slack then dual 0.0 (if feasible) // dj for key is zero so that defines dual on set int i; double *dj = model->djRegion(); int numberColumns = model->numberColumns(); double infeasibilityCost = model->infeasibilityCost(); for (i = 0; i < numberSets_; i++) { int kColumn = keyVariable_[i]; if (kColumn < numberColumns) { // dj without set double value = dj[kColumn]; // Now subtract out from all dj[kColumn] = 0.0; int iColumn = next_[kColumn]; // modify all non-key variables while (iColumn >= 0) { dj[iColumn] -= value; iColumn = next_[iColumn]; } } else { // slack key - may not be feasible assert(getStatus(i) == ClpSimplex::basic); // negative as -1.0 for slack double value = -weight(i) * infeasibilityCost; if (value) { int iColumn = next_[kColumn]; // modify all non-key variables basic while (iColumn >= 0) { dj[iColumn] -= value; iColumn = next_[iColumn]; } } } } } break; // as 1 but check slacks and compute djs case 2: { // 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; // make sure fromIndex will not confuse pricing fromIndex_[0] = -1; possiblePivotKey_ = -1; // Create array int numberColumns = model->numberColumns(); int *pivotVariable = model->pivotVariable(); int numberRows = model->numberRows(); for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; if (iPivot < numberColumns) backToPivotRow_[iPivot] = i; } if (noCheck_ >= 0) { if (infeasibilityWeight_ != model->infeasibilityCost()) { // don't bother checking sumDualInfeasibilities_ = 100.0; numberDualInfeasibilities_ = 1; sumOfRelaxedDualInfeasibilities_ = 100.0; return; } } double *dj = model->djRegion(); double *dual = model->dualRowSolution(); double *cost = model->costRegion(); ClpSimplex::Status iStatus; const int *columnLength = matrix_->getVectorLengths(); const CoinBigIndex *columnStart = matrix_->getVectorStarts(); const int *row = matrix_->getIndices(); const double *elementByColumn = matrix_->getElements(); double infeasibilityCost = model->infeasibilityCost(); sumDualInfeasibilities_ = 0.0; numberDualInfeasibilities_ = 0; 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; sumOfRelaxedDualInfeasibilities_ = 0.0; for (i = 0; i < numberSets_; i++) { int kColumn = keyVariable_[i]; if (kColumn < numberColumns) { // dj without set double value = cost[kColumn]; for (CoinBigIndex j = columnStart[kColumn]; j < columnStart[kColumn] + columnLength[kColumn]; j++) { int iRow = row[j]; value -= dual[iRow] * elementByColumn[j]; } // Now subtract out from all dj[kColumn] -= value; int stop = -(kColumn + 1); kColumn = next_[kColumn]; while (kColumn != stop) { if (kColumn < 0) kColumn = -kColumn - 1; double djValue = dj[kColumn] - value; dj[kColumn] = djValue; ; double infeasibility = 0.0; iStatus = model->getStatus(kColumn); if (iStatus == ClpSimplex::atLowerBound) { if (djValue < -dualTolerance) infeasibility = -djValue - dualTolerance; } else if (iStatus == ClpSimplex::atUpperBound) { // at upper bound if (djValue > dualTolerance) infeasibility = djValue - dualTolerance; } if (infeasibility > 0.0) { sumDualInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += infeasibility; numberDualInfeasibilities_++; } kColumn = next_[kColumn]; } // check slack iStatus = getStatus(i); assert(iStatus != ClpSimplex::basic); double infeasibility = 0.0; // dj of slack is -(-1.0)value if (iStatus == ClpSimplex::atLowerBound) { if (value < -dualTolerance) infeasibility = -value - dualTolerance; } else if (iStatus == ClpSimplex::atUpperBound) { // at upper bound if (value > dualTolerance) infeasibility = value - dualTolerance; } if (infeasibility > 0.0) { sumDualInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += infeasibility; numberDualInfeasibilities_++; } } else { // slack key - may not be feasible assert(getStatus(i) == ClpSimplex::basic); // negative as -1.0 for slack double value = -weight(i) * infeasibilityCost; if (value) { // Now subtract out from all int kColumn = i + numberColumns; int stop = -(kColumn + 1); kColumn = next_[kColumn]; while (kColumn != stop) { if (kColumn < 0) kColumn = -kColumn - 1; double djValue = dj[kColumn] - value; dj[kColumn] = djValue; ; double infeasibility = 0.0; iStatus = model->getStatus(kColumn); if (iStatus == ClpSimplex::atLowerBound) { if (djValue < -dualTolerance) infeasibility = -djValue - dualTolerance; } else if (iStatus == ClpSimplex::atUpperBound) { // at upper bound if (djValue > dualTolerance) infeasibility = djValue - dualTolerance; } if (infeasibility > 0.0) { sumDualInfeasibilities_ += infeasibility; if (infeasibility > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += infeasibility; numberDualInfeasibilities_++; } kColumn = next_[kColumn]; } } } } // and get statistics for column generation synchronize(model, 4); 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: { // First compute new costs etc for interesting gubs int iLook = 0; int iSet = fromIndex_[0]; double primalTolerance = model->primalTolerance(); const double *cost = model->costRegion(); double *solution = model->solutionRegion(); double infeasibilityCost = model->infeasibilityCost(); int numberColumns = model->numberColumns(); int numberChanged = 0; int *pivotVariable = model->pivotVariable(); while (iSet >= 0) { int key = keyVariable_[iSet]; double value = 0.0; // sum over all except key if ((gubType_ & 8) != 0) { int iColumn = next_[key]; // sum all non-key variables while (iColumn >= 0) { value += solution[iColumn]; iColumn = next_[iColumn]; } } else { // bounds exist - sum over all except key int stop = -(key + 1); int iColumn = next_[key]; // sum all non-key variables while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; value += solution[iColumn]; iColumn = next_[iColumn]; } } double costChange; double oldCost = changeCost_[iLook]; if (key < numberColumns) { assert(getStatus(iSet) != ClpSimplex::basic); double sol; if (getStatus(iSet) == ClpSimplex::atUpperBound) sol = upper_[iSet] - value; else sol = lower_[iSet] - value; solution[key] = sol; // fix up cost model->nonLinearCost()->setOne(key, sol); #ifdef CLP_DEBUG_PRINT printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol, cost[key], oldCost); #endif costChange = cost[key] - oldCost; } else { // slack is key if (value > upper_[iSet] + primalTolerance) { setAbove(iSet); } else if (value < lower_[iSet] - primalTolerance) { setBelow(iSet); } else { setFeasible(iSet); } // negative as -1.0 for slack costChange = -weight(iSet) * infeasibilityCost - oldCost; #ifdef CLP_DEBUG_PRINT printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value, weight(iSet) * infeasibilityCost, oldCost); #endif } if (costChange) { fromIndex_[numberChanged] = iSet; toIndex_[iSet] = numberChanged; changeCost_[numberChanged++] = costChange; } iSet = fromIndex_[++iLook]; } if (numberChanged || possiblePivotKey_ >= 0) { // first do those in list already int number = array->getNumElements(); array->setPacked(); int i; double *work = array->denseVector(); int *which = array->getIndices(); for (i = 0; i < number; i++) { int iRow = which[i]; int iPivot = pivotVariable[iRow]; if (iPivot < numberColumns) { int iSet = backward_[iPivot]; if (iSet >= 0 && toIndex_[iSet] >= 0) { double newValue = work[i] + changeCost_[toIndex_[iSet]]; if (!newValue) newValue = 1.0e-100; work[i] = newValue; // mark as done backward_[iPivot] = -1; } } if (possiblePivotKey_ == iRow) { double newValue = work[i] - model->dualIn(); if (!newValue) newValue = 1.0e-100; work[i] = newValue; possiblePivotKey_ = -1; } } // now do rest and clean up for (i = 0; i < numberChanged; i++) { int iSet = fromIndex_[i]; int key = keyVariable_[iSet]; int iColumn = next_[key]; double change = changeCost_[i]; while (iColumn >= 0) { if (backward_[iColumn] >= 0) { int iRow = backToPivotRow_[iColumn]; assert(iRow >= 0); work[number] = change; if (possiblePivotKey_ == iRow) { double newValue = work[number] - model->dualIn(); if (!newValue) newValue = 1.0e-100; work[number] = newValue; possiblePivotKey_ = -1; } which[number++] = iRow; } else { // reset backward_[iColumn] = iSet; } iColumn = next_[iColumn]; } toIndex_[iSet] = -1; } if (possiblePivotKey_ >= 0) { work[number] = -model->dualIn(); which[number++] = possiblePivotKey_; possiblePivotKey_ = -1; } fromIndex_[0] = -1; array->setNumElements(number); } } break; } } // This is local to Gub to allow synchronization when status is good int ClpGubMatrix::synchronize(ClpSimplex *, int) { return 0; } /* general utility function for dealing with dynamic constraints mode=n see ClpGubMatrix.hpp for definition Remember to update here when settled down */ int ClpGubMatrix::generalExpanded(ClpSimplex *model, int mode, int &number) { int returnCode = 0; int numberColumns = model->numberColumns(); switch (mode) { // Fill in pivotVariable but not for key variables case 0: { if (!next_) { // do ordering assert(!rhsOffset_); // create and do gub crash useEffectiveRhs(model, false); } int i; int numberBasic = number; // 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) { int iSet = backward_[i]; if (iSet < 0 || i != keyVariable_[iSet]) pivotVariable[numberBasic++] = i; } } number = numberBasic; if (model->numberIterations()) assert(number == model->numberRows()); } break; // Make all key variables basic case 1: { int i; for (i = 0; i < numberSets_; i++) { int iColumn = keyVariable_[i]; if (iColumn < numberColumns) model->setColumnStatus(iColumn, ClpSimplex::basic); } } break; // Do initial extra rows + maximum basic case 2: { returnCode = getNumRows() + 1; number = model->numberRows() + numberSets_; } break; // Before normal replaceColumn case 3: { int sequenceIn = model->sequenceIn(); int sequenceOut = model->sequenceOut(); int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); int pivotRow = model->pivotRow(); if (gubSlackIn_ >= 0) assert(sequenceIn > numberRows + numberColumns); if (sequenceIn == sequenceOut) return -1; int iSetIn = -1; int iSetOut = -1; if (sequenceOut < numberColumns) { iSetOut = backward_[sequenceOut]; } else if (sequenceOut >= numberRows + numberColumns) { assert(pivotRow >= numberRows); int iExtra = pivotRow - numberRows; assert(iExtra >= 0); if (iSetOut < 0) iSetOut = fromIndex_[iExtra]; else assert(iSetOut == fromIndex_[iExtra]); } if (sequenceIn < numberColumns) { iSetIn = backward_[sequenceIn]; } else if (gubSlackIn_ >= 0) { iSetIn = gubSlackIn_; } possiblePivotKey_ = -1; number = 0; // say do ordinary int *pivotVariable = model->pivotVariable(); if (pivotRow >= numberRows) { int iExtra = pivotRow - numberRows; //const int * length = matrix_->getVectorLengths(); assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]); int incomingColumn = sequenceIn; // to be used in updates if (iSetIn != iSetOut) { // We need to find a possible pivot for incoming // look through rowArray_[1] int n = model->rowArray(1)->getNumElements(); int *which = model->rowArray(1)->getIndices(); double *array = model->rowArray(1)->denseVector(); double bestAlpha = 1.0e-5; //int shortest=numberRows+1; for (int i = 0; i < n; i++) { int iRow = which[i]; int iPivot = pivotVariable[iRow]; if (iPivot < numberColumns && backward_[iPivot] == iSetOut) { if (fabs(array[i]) > fabs(bestAlpha)) { bestAlpha = array[i]; possiblePivotKey_ = iRow; } } } assert(possiblePivotKey_ >= 0); // could set returnCode=4 number = 1; if (sequenceIn >= numberRows + numberColumns) { number = 3; // need swap as gub slack in and must become key // is this best way int key = keyVariable_[iSetIn]; assert(key < numberColumns); // check other basic int iColumn = next_[key]; // set new key to be used by unpack keyVariable_[iSetIn] = iSetIn + numberColumns; // change cost in changeCost { int iLook = 0; int iSet = fromIndex_[0]; while (iSet >= 0) { if (iSet == iSetIn) { changeCost_[iLook] = 0.0; break; } iSet = fromIndex_[++iLook]; } } while (iColumn >= 0) { if (iColumn != sequenceOut) { // need partial ftran and skip accuracy check in replaceColumn #ifdef CLP_DEBUG_PRINT printf("TTTTTry 5\n"); #endif int iRow = backToPivotRow_[iColumn]; assert(iRow >= 0); unpack(model, model->rowArray(3), iColumn); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3)); double alpha = model->rowArray(3)->denseVector()[iRow]; //if (!alpha) //printf("zero alpha a\n"); int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(3), iRow, alpha); returnCode = CoinMax(updateStatus, returnCode); model->rowArray(3)->clear(); if (returnCode) break; } iColumn = next_[iColumn]; } if (!returnCode) { // now factorization looks as if key is out // pivot back in #ifdef CLP_DEBUG_PRINT printf("TTTTTry 6\n"); #endif unpack(model, model->rowArray(3), key); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3)); pivotRow = possiblePivotKey_; double alpha = model->rowArray(3)->denseVector()[pivotRow]; //if (!alpha) //printf("zero alpha b\n"); int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(3), pivotRow, alpha); returnCode = CoinMax(updateStatus, returnCode); model->rowArray(3)->clear(); } // restore key keyVariable_[iSetIn] = key; // now alternate column can replace key on out incomingColumn = pivotVariable[possiblePivotKey_]; } else { #ifdef CLP_DEBUG_PRINT printf("TTTTTTry 4 %d\n", possiblePivotKey_); #endif int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(1), possiblePivotKey_, bestAlpha); returnCode = CoinMax(updateStatus, returnCode); incomingColumn = pivotVariable[possiblePivotKey_]; } //returnCode=4; // need swap } else { // key swap number = -1; } int key = keyVariable_[iSetOut]; if (key < numberColumns) assert(key == sequenceOut); // check if any other basic int iColumn = next_[key]; if (returnCode) iColumn = -1; // skip if error on previous // set new key to be used by unpack if (incomingColumn < numberColumns) keyVariable_[iSetOut] = incomingColumn; else keyVariable_[iSetOut] = iSetIn + numberColumns; double *cost = model->costRegion(); if (possiblePivotKey_ < 0) { double dj = model->djRegion()[sequenceIn] - cost[sequenceIn]; changeCost_[iExtra] = -dj; #ifdef CLP_DEBUG_PRINT printf("modifying changeCost %d by %g - cost %g\n", iExtra, dj, cost[sequenceIn]); #endif } while (iColumn >= 0) { if (iColumn != incomingColumn) { number = -2; // need partial ftran and skip accuracy check in replaceColumn #ifdef CLP_DEBUG_PRINT printf("TTTTTTry 1\n"); #endif int iRow = backToPivotRow_[iColumn]; assert(iRow >= 0 && iRow < numberRows); unpack(model, model->rowArray(3), iColumn); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3)); double *array = model->rowArray(3)->denseVector(); double alpha = array[iRow]; //if (!alpha) //printf("zero alpha d\n"); int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(3), iRow, alpha); returnCode = CoinMax(updateStatus, returnCode); model->rowArray(3)->clear(); if (returnCode) break; } iColumn = next_[iColumn]; } // restore key keyVariable_[iSetOut] = key; } else if (sequenceIn >= numberRows + numberColumns) { number = 2; //returnCode=4; // need swap as gub slack in and must become key // is this best way int key = keyVariable_[iSetIn]; assert(key < numberColumns); // check other basic int iColumn = next_[key]; // set new key to be used by unpack keyVariable_[iSetIn] = iSetIn + numberColumns; // change cost in changeCost { int iLook = 0; int iSet = fromIndex_[0]; while (iSet >= 0) { if (iSet == iSetIn) { changeCost_[iLook] = 0.0; break; } iSet = fromIndex_[++iLook]; } } while (iColumn >= 0) { if (iColumn != sequenceOut) { // need partial ftran and skip accuracy check in replaceColumn #ifdef CLP_DEBUG_PRINT printf("TTTTTry 2\n"); #endif int iRow = backToPivotRow_[iColumn]; assert(iRow >= 0); unpack(model, model->rowArray(3), iColumn); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3)); double alpha = model->rowArray(3)->denseVector()[iRow]; //if (!alpha) //printf("zero alpha e\n"); int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(3), iRow, alpha); returnCode = CoinMax(updateStatus, returnCode); model->rowArray(3)->clear(); if (returnCode) break; } iColumn = next_[iColumn]; } if (!returnCode) { // now factorization looks as if key is out // pivot back in #ifdef CLP_DEBUG_PRINT printf("TTTTTry 3\n"); #endif unpack(model, model->rowArray(3), key); model->factorization()->updateColumnFT(model->rowArray(2), model->rowArray(3)); double alpha = model->rowArray(3)->denseVector()[pivotRow]; //if (!alpha) //printf("zero alpha f\n"); int updateStatus = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(3), pivotRow, alpha); returnCode = CoinMax(updateStatus, returnCode); model->rowArray(3)->clear(); } // restore key keyVariable_[iSetIn] = key; } else { // normal - but might as well do here returnCode = model->factorization()->replaceColumn(model, model->rowArray(2), model->rowArray(1), model->pivotRow(), model->alpha()); } } #ifdef CLP_DEBUG_PRINT printf("Update type after %d - status %d - pivot row %d\n", number, returnCode, model->pivotRow()); #endif // see if column generation says time to re-factorize returnCode = CoinMax(returnCode, synchronize(model, 5)); 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: { synchronize(model, 0); CoinMemcpyN(status_, numberSets_, saveStatus_); CoinMemcpyN(keyVariable_, numberSets_, savedKeyVariable_); } break; // restore status case 6: { CoinMemcpyN(saveStatus_, numberSets_, status_); CoinMemcpyN(savedKeyVariable_, numberSets_, keyVariable_); // restore firstAvailable_ synchronize(model, 7); // redo next_ int i; int *last = new int[numberSets_]; for (i = 0; i < numberSets_; i++) { int iKey = keyVariable_[i]; assert(iKey >= numberColumns || backward_[iKey] == i); last[i] = iKey; // make sure basic //if (iKeysetStatus(iKey,ClpSimplex::basic); } for (i = 0; i < numberColumns; i++) { int iSet = backward_[i]; if (iSet >= 0) { next_[last[iSet]] = i; last[iSet] = i; } } for (i = 0; i < numberSets_; i++) { next_[last[i]] = -(keyVariable_[i] + 1); redoSet(model, keyVariable_[i], keyVariable_[i], i); } delete[] last; // redo pivotVariable int *pivotVariable = model->pivotVariable(); int iRow; int numberBasic = 0; int numberRows = model->numberRows(); for (iRow = 0; iRow < numberRows; iRow++) { if (model->getRowStatus(iRow) == ClpSimplex::basic) { numberBasic++; pivotVariable[iRow] = iRow + numberColumns; } else { pivotVariable[iRow] = -1; } } i = 0; int iColumn; for (iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getStatus(iColumn) == ClpSimplex::basic) { int iSet = backward_[iColumn]; if (iSet < 0 || keyVariable_[iSet] != iColumn) { while (pivotVariable[i] >= 0) { i++; assert(i < numberRows); } pivotVariable[i] = iColumn; backToPivotRow_[iColumn] = i; numberBasic++; } } } assert(numberBasic == numberRows); rhsOffset(model, true); } break; // flag a variable case 7: { assert(number == model->sequenceIn()); synchronize(model, 1); synchronize(model, 8); } break; // unflag all variables case 8: { returnCode = synchronize(model, 2); } break; // redo costs in primal case 9: { returnCode = synchronize(model, 3); } break; // return 1 if there may be changing bounds on variable (column generation) case 10: { returnCode = synchronize(model, 6); } break; // make sure set is clean case 11: { assert(number == model->sequenceIn()); returnCode = synchronize(model, 8); } break; default: break; } return returnCode; } // Sets up an effective RHS and does gub crash if needed void ClpGubMatrix::useEffectiveRhs(ClpSimplex *model, bool cheapest) { // Do basis - cheapest or slack if feasible (unless cheapest set) int longestSet = 0; int iSet; for (iSet = 0; iSet < numberSets_; iSet++) longestSet = CoinMax(longestSet, end_[iSet] - start_[iSet]); double *upper = new double[longestSet + 1]; double *cost = new double[longestSet + 1]; double *lower = new double[longestSet + 1]; double *solution = new double[longestSet + 1]; assert(!next_); int numberColumns = getNumCols(); const int *columnLength = matrix_->getVectorLengths(); const double *columnLower = model->lowerRegion(); const double *columnUpper = model->upperRegion(); double *columnSolution = model->solutionRegion(); const double *objective = model->costRegion(); int numberRows = getNumRows(); toIndex_ = new int[numberSets_]; for (iSet = 0; iSet < numberSets_; iSet++) toIndex_[iSet] = -1; fromIndex_ = new int[getNumRows() + 1]; double tolerance = model->primalTolerance(); bool noNormalBounds = true; gubType_ &= ~8; bool gotBasis = false; for (iSet = 0; iSet < numberSets_; iSet++) { if (keyVariable_[iSet] < numberColumns) gotBasis = true; CoinBigIndex j; CoinBigIndex iStart = start_[iSet]; CoinBigIndex iEnd = end_[iSet]; for (j = iStart; j < iEnd; j++) { if (columnLower[j] && columnLower[j] > -1.0e20) noNormalBounds = false; if (columnUpper[j] && columnUpper[j] < 1.0e20) noNormalBounds = false; } } if (noNormalBounds) gubType_ |= 8; if (!gotBasis) { for (iSet = 0; iSet < numberSets_; iSet++) { CoinBigIndex j; int numberBasic = 0; int iBasic = -1; CoinBigIndex iStart = start_[iSet]; CoinBigIndex iEnd = end_[iSet]; // find one with smallest length int smallest = numberRows + 1; double value = 0.0; for (j = iStart; j < iEnd; j++) { if (model->getStatus(j) == ClpSimplex::basic) { if (columnLength[j] < smallest) { smallest = columnLength[j]; iBasic = j; } numberBasic++; } value += columnSolution[j]; } bool done = false; if (numberBasic > 1 || (numberBasic == 1 && getStatus(iSet) == ClpSimplex::basic)) { if (getStatus(iSet) == ClpSimplex::basic) iBasic = iSet + numberColumns; // slack key - use done = true; } else if (numberBasic == 1) { // see if can be key double thisSolution = columnSolution[iBasic]; if (thisSolution > columnUpper[iBasic]) { value -= thisSolution - columnUpper[iBasic]; thisSolution = columnUpper[iBasic]; columnSolution[iBasic] = thisSolution; } if (thisSolution < columnLower[iBasic]) { value -= thisSolution - columnLower[iBasic]; thisSolution = columnLower[iBasic]; columnSolution[iBasic] = thisSolution; } // try setting slack to a bound assert(upper_[iSet] < 1.0e20 || lower_[iSet] > -1.0e20); double cost1 = COIN_DBL_MAX; int whichBound = -1; if (upper_[iSet] < 1.0e20) { // try slack at ub double newBasic = thisSolution + upper_[iSet] - value; if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) { // can go whichBound = 1; cost1 = newBasic * objective[iBasic]; // But if exact then may be good solution if (fabs(upper_[iSet] - value) < tolerance) cost1 = -COIN_DBL_MAX; } } if (lower_[iSet] > -1.0e20) { // try slack at lb double newBasic = thisSolution + lower_[iSet] - value; if (newBasic >= columnLower[iBasic] - tolerance && newBasic <= columnUpper[iBasic] + tolerance) { // can go but is it cheaper double cost2 = newBasic * objective[iBasic]; // But if exact then may be good solution if (fabs(lower_[iSet] - value) < tolerance) cost2 = -COIN_DBL_MAX; if (cost2 < cost1) whichBound = 0; } } if (whichBound != -1) { // key done = true; if (whichBound) { // slack to upper columnSolution[iBasic] = thisSolution + upper_[iSet] - value; setStatus(iSet, ClpSimplex::atUpperBound); } else { // slack to lower columnSolution[iBasic] = thisSolution + lower_[iSet] - value; setStatus(iSet, ClpSimplex::atLowerBound); } } } if (!done) { if (!cheapest) { // see if slack can be key if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) { done = true; setStatus(iSet, ClpSimplex::basic); iBasic = iSet + numberColumns; } } if (!done) { // set non basic if there was one if (iBasic >= 0) model->setStatus(iBasic, ClpSimplex::atLowerBound); // find cheapest int numberInSet = iEnd - iStart; CoinMemcpyN(columnLower + iStart, numberInSet, lower); CoinMemcpyN(columnUpper + iStart, numberInSet, upper); CoinMemcpyN(columnSolution + iStart, numberInSet, solution); // and slack iBasic = numberInSet; solution[iBasic] = -value; lower[iBasic] = -upper_[iSet]; upper[iBasic] = -lower_[iSet]; int kphase; if (value >= lower_[iSet] - tolerance && value <= upper_[iSet] + tolerance) { // feasible kphase = 1; cost[iBasic] = 0.0; CoinMemcpyN(objective + iStart, numberInSet, cost); } else { // infeasible kphase = 0; // remember bounds are flipped so opposite to natural if (value < lower_[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; CoinMemcpyN(objective + iStart, numberInSet, cost); } // 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) } } } } // convert iBasic back and do bounds if (iBasic == numberInSet) { // slack basic setStatus(iSet, ClpSimplex::basic); iBasic = iSet + numberColumns; } else { iBasic += start_[iSet]; model->setStatus(iBasic, ClpSimplex::basic); // 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(); } for (j = iStart; j < iEnd; j++) { if (model->getStatus(j) != ClpSimplex::basic) { int inSet = j - iStart; columnSolution[j] = solution[inSet]; if (upper[inSet] == lower[inSet]) model->setStatus(j, ClpSimplex::isFixed); else if (solution[inSet] == upper[inSet]) model->setStatus(j, ClpSimplex::atUpperBound); else if (solution[inSet] == lower[inSet]) model->setStatus(j, ClpSimplex::atLowerBound); } } } } keyVariable_[iSet] = iBasic; } } delete[] lower; delete[] solution; delete[] upper; delete[] cost; // make sure matrix is in good shape matrix_->orderMatrix(); // create effective rhs delete[] rhsOffset_; rhsOffset_ = new double[numberRows]; delete[] next_; next_ = new int[numberColumns + numberSets_ + 2 * longestSet]; char *mark = new char[numberColumns]; memset(mark, 0, numberColumns); for (int iColumn = 0; iColumn < numberColumns; iColumn++) next_[iColumn] = COIN_INT_MAX; int i; int *keys = new int[numberSets_]; for (i = 0; i < numberSets_; i++) keys[i] = COIN_INT_MAX; // set up chains for (i = 0; i < numberColumns; i++) { if (model->getStatus(i) == ClpSimplex::basic) mark[i] = 1; int iSet = backward_[i]; if (iSet >= 0) { int iNext = keys[iSet]; next_[i] = iNext; keys[iSet] = i; } } for (i = 0; i < numberSets_; i++) { int j; if (getStatus(i) != ClpSimplex::basic) { // make sure fixed if it is if (upper_[i] == lower_[i]) setStatus(i, ClpSimplex::isFixed); // slack not key - choose one with smallest length int smallest = numberRows + 1; int key = -1; j = keys[i]; if (j != COIN_INT_MAX) { while (1) { if (mark[j] && columnLength[j] < smallest && !gotBasis) { key = j; smallest = columnLength[j]; } if (next_[j] != COIN_INT_MAX) { j = next_[j]; } else { // correct end next_[j] = -(keys[i] + 1); break; } } } else { next_[i + numberColumns] = -(numberColumns + i + 1); } if (gotBasis) key = keyVariable_[i]; if (key >= 0) { keyVariable_[i] = key; } else { // nothing basic - make slack key //((ClpGubMatrix *)this)->setStatus(i,ClpSimplex::basic); // fudge to avoid const problem status_[i] = 1; } } else { // slack key keyVariable_[i] = numberColumns + i; int j; double sol = 0.0; j = keys[i]; if (j != COIN_INT_MAX) { while (1) { sol += columnSolution[j]; if (next_[j] != COIN_INT_MAX) { j = next_[j]; } else { // correct end next_[j] = -(keys[i] + 1); break; } } } else { next_[i + numberColumns] = -(numberColumns + i + 1); } if (sol > upper_[i] + tolerance) { setAbove(i); } else if (sol < lower_[i] - tolerance) { setBelow(i); } else { setFeasible(i); } } // Create next_ int key = keyVariable_[i]; redoSet(model, key, keys[i], i); } delete[] keys; delete[] mark; rhsOffset(model, true); } // redoes next_ for a set. void ClpGubMatrix::redoSet(ClpSimplex *model, int newKey, int oldKey, int iSet) { int numberColumns = model->numberColumns(); int *save = next_ + numberColumns + numberSets_; int number = 0; int stop = -(oldKey + 1); int j = next_[oldKey]; while (j != stop) { if (j < 0) j = -j - 1; if (j != newKey) save[number++] = j; j = next_[j]; } // and add oldkey if (newKey != oldKey) save[number++] = oldKey; // now do basic int lastMarker = -(newKey + 1); keyVariable_[iSet] = newKey; next_[newKey] = lastMarker; int last = newKey; for (j = 0; j < number; j++) { int iColumn = save[j]; if (iColumn < numberColumns) { if (model->getStatus(iColumn) == ClpSimplex::basic) { next_[last] = iColumn; next_[iColumn] = lastMarker; last = iColumn; } } } // now add in non-basic for (j = 0; j < number; j++) { int iColumn = save[j]; if (iColumn < numberColumns) { if (model->getStatus(iColumn) != ClpSimplex::basic) { next_[last] = -(iColumn + 1); next_[iColumn] = lastMarker; last = iColumn; } } } } /* 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 * ClpGubMatrix::rhsOffset(ClpSimplex *model, bool forceRefresh, bool #ifdef CLP_DEBUG check #endif ) { //forceRefresh=true; if (rhsOffset_) { #ifdef CLP_DEBUG if (check) { // no need - but check anyway // zero out basic int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); double *solution = new double[numberColumns]; double *rhs = new double[numberRows]; CoinMemcpyN(model->solutionRegion(), numberColumns, solution); CoinZeroN(rhs, numberRows); int iRow; for (int iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) == ClpSimplex::basic) solution[iColumn] = 0.0; } for (int iSet = 0; iSet < numberSets_; iSet++) { int iColumn = keyVariable_[iSet]; if (iColumn < numberColumns) solution[iColumn] = 0.0; } times(-1.0, solution, rhs); delete[] solution; const double *columnSolution = model->solutionRegion(); // and now subtract out non basic ClpSimplex::Status iStatus; for (int iSet = 0; iSet < numberSets_; iSet++) { int iColumn = keyVariable_[iSet]; if (iColumn < numberColumns) { double b = 0.0; // key is structural - where is slack iStatus = getStatus(iSet); assert(iStatus != ClpSimplex::basic); if (iStatus == ClpSimplex::atLowerBound) b = lower_[iSet]; else b = upper_[iSet]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(iColumn + 1); int jColumn = next_[iColumn]; // sum all non-basic variables - first skip basic while (jColumn >= 0) jColumn = next_[jColumn]; while (jColumn != stop) { assert(jColumn < 0); jColumn = -jColumn - 1; b -= columnSolution[jColumn]; jColumn = next_[jColumn]; } } // subtract out ClpPackedMatrix::add(model, rhs, iColumn, -b); } } 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]); } delete[] rhs; } #endif if (forceRefresh || (refreshFrequency_ && model->numberIterations() >= lastRefresh_ + refreshFrequency_)) { // zero out basic int numberRows = model->numberRows(); int numberColumns = model->numberColumns(); double *solution = new double[numberColumns]; CoinMemcpyN(model->solutionRegion(), numberColumns, solution); CoinZeroN(rhsOffset_, numberRows); for (int iColumn = 0; iColumn < numberColumns; iColumn++) { if (model->getColumnStatus(iColumn) == ClpSimplex::basic) solution[iColumn] = 0.0; } int iSet; for (iSet = 0; iSet < numberSets_; iSet++) { int iColumn = keyVariable_[iSet]; if (iColumn < numberColumns) solution[iColumn] = 0.0; } times(-1.0, solution, rhsOffset_); delete[] solution; lastRefresh_ = model->numberIterations(); const double *columnSolution = model->solutionRegion(); // and now subtract out non basic ClpSimplex::Status iStatus; for (iSet = 0; iSet < numberSets_; iSet++) { int iColumn = keyVariable_[iSet]; if (iColumn < numberColumns) { double b = 0.0; // key is structural - where is slack iStatus = getStatus(iSet); assert(iStatus != ClpSimplex::basic); if (iStatus == ClpSimplex::atLowerBound) b = lower_[iSet]; else b = upper_[iSet]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(iColumn + 1); int jColumn = next_[iColumn]; // sum all non-basic variables - first skip basic while (jColumn >= 0) jColumn = next_[jColumn]; while (jColumn != stop) { assert(jColumn < 0); jColumn = -jColumn - 1; b -= columnSolution[jColumn]; jColumn = next_[jColumn]; } } // subtract out if (b) ClpPackedMatrix::add(model, rhsOffset_, iColumn, -b); } } } } return rhsOffset_; } /* update information for a pivot (and effective rhs) */ int ClpGubMatrix::updatePivot(ClpSimplex *model, double oldInValue, double /*oldOutValue*/) { int sequenceIn = model->sequenceIn(); int sequenceOut = model->sequenceOut(); double *solution = model->solutionRegion(); int numberColumns = model->numberColumns(); int numberRows = model->numberRows(); int pivotRow = model->pivotRow(); int iSetIn; // Correct sequence in trueSequenceIn_ = sequenceIn; if (sequenceIn < numberColumns) { iSetIn = backward_[sequenceIn]; } else if (sequenceIn < numberColumns + numberRows) { iSetIn = -1; } else { iSetIn = gubSlackIn_; trueSequenceIn_ = numberColumns + numberRows + iSetIn; } int iSetOut = -1; trueSequenceOut_ = sequenceOut; if (sequenceOut < numberColumns) { iSetOut = backward_[sequenceOut]; } else if (sequenceOut >= numberRows + numberColumns) { assert(pivotRow >= numberRows); int iExtra = pivotRow - numberRows; assert(iExtra >= 0); if (iSetOut < 0) iSetOut = fromIndex_[iExtra]; else assert(iSetOut == fromIndex_[iExtra]); trueSequenceOut_ = numberColumns + numberRows + iSetOut; } if (rhsOffset_) { // update effective rhs if (sequenceIn == sequenceOut) { assert(sequenceIn < numberRows + numberColumns); // should be easy to deal with if (sequenceIn < numberColumns) add(model, rhsOffset_, sequenceIn, oldInValue - solution[sequenceIn]); } else { if (sequenceIn < numberColumns) { // we need to test if WILL be key ClpPackedMatrix::add(model, rhsOffset_, sequenceIn, oldInValue); if (iSetIn >= 0) { // old contribution to rhsOffset_ int key = keyVariable_[iSetIn]; if (key < numberColumns) { double oldB = 0.0; ClpSimplex::Status iStatus = getStatus(iSetIn); if (iStatus == ClpSimplex::atLowerBound) oldB = lower_[iSetIn]; else oldB = upper_[iSetIn]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(key + 1); int iColumn = next_[key]; // skip basic while (iColumn >= 0) iColumn = next_[iColumn]; // sum all non-key variables while (iColumn != stop) { assert(iColumn < 0); iColumn = -iColumn - 1; if (iColumn == sequenceIn) oldB -= oldInValue; else if (iColumn != sequenceOut) oldB -= solution[iColumn]; iColumn = next_[iColumn]; } } if (oldB) ClpPackedMatrix::add(model, rhsOffset_, key, oldB); } } } else if (sequenceIn < numberRows + numberColumns) { //rhsOffset_[sequenceIn-numberColumns] -= oldInValue; } else { #ifdef CLP_DEBUG_PRINT printf("** in is key slack %d\n", sequenceIn); #endif // old contribution to rhsOffset_ int key = keyVariable_[iSetIn]; if (key < numberColumns) { double oldB = 0.0; ClpSimplex::Status iStatus = getStatus(iSetIn); if (iStatus == ClpSimplex::atLowerBound) oldB = lower_[iSetIn]; else oldB = upper_[iSetIn]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(key + 1); int iColumn = next_[key]; // skip basic while (iColumn >= 0) iColumn = next_[iColumn]; // sum all non-key variables while (iColumn != stop) { assert(iColumn < 0); iColumn = -iColumn - 1; if (iColumn != sequenceOut) oldB -= solution[iColumn]; iColumn = next_[iColumn]; } } if (oldB) ClpPackedMatrix::add(model, rhsOffset_, key, oldB); } } if (sequenceOut < numberColumns) { ClpPackedMatrix::add(model, rhsOffset_, sequenceOut, -solution[sequenceOut]); if (iSetOut >= 0) { // old contribution to rhsOffset_ int key = keyVariable_[iSetOut]; if (key < numberColumns && iSetIn != iSetOut) { double oldB = 0.0; ClpSimplex::Status iStatus = getStatus(iSetOut); if (iStatus == ClpSimplex::atLowerBound) oldB = lower_[iSetOut]; else oldB = upper_[iSetOut]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(key + 1); int iColumn = next_[key]; // skip basic while (iColumn >= 0) iColumn = next_[iColumn]; // sum all non-key variables while (iColumn != stop) { assert(iColumn < 0); iColumn = -iColumn - 1; if (iColumn == sequenceIn) oldB -= oldInValue; else if (iColumn != sequenceOut) oldB -= solution[iColumn]; iColumn = next_[iColumn]; } } if (oldB) ClpPackedMatrix::add(model, rhsOffset_, key, oldB); } } } else if (sequenceOut < numberRows + numberColumns) { //rhsOffset_[sequenceOut-numberColumns] -= -solution[sequenceOut]; } else { #ifdef CLP_DEBUG_PRINT printf("** out is key slack %d\n", sequenceOut); #endif assert(pivotRow >= numberRows); } } } int *pivotVariable = model->pivotVariable(); // may need to deal with key // Also need coding to mark/allow key slack entering if (pivotRow >= numberRows) { assert(sequenceOut >= numberRows + numberColumns || sequenceOut == keyVariable_[iSetOut]); #ifdef CLP_DEBUG_PRINT if (sequenceIn >= numberRows + numberColumns) printf("key slack %d in, set out %d\n", gubSlackIn_, iSetOut); printf("** danger - key out for set %d in %d (set %d)\n", iSetOut, sequenceIn, iSetIn); #endif // if slack out mark correctly if (sequenceOut >= numberRows + numberColumns) { double value = model->valueOut(); if (value == upper_[iSetOut]) { setStatus(iSetOut, ClpSimplex::atUpperBound); } else if (value == lower_[iSetOut]) { setStatus(iSetOut, ClpSimplex::atLowerBound); } else { if (fabs(value - upper_[iSetOut]) < fabs(value - lower_[iSetOut])) { setStatus(iSetOut, ClpSimplex::atUpperBound); } else { setStatus(iSetOut, ClpSimplex::atLowerBound); } } if (upper_[iSetOut] == lower_[iSetOut]) setStatus(iSetOut, ClpSimplex::isFixed); setFeasible(iSetOut); } if (iSetOut == iSetIn) { // key swap int key; if (sequenceIn >= numberRows + numberColumns) { key = numberColumns + iSetIn; setStatus(iSetIn, ClpSimplex::basic); } else { key = sequenceIn; } redoSet(model, key, keyVariable_[iSetIn], iSetIn); } else { // key was chosen assert(possiblePivotKey_ >= 0 && possiblePivotKey_ < numberRows); int key = pivotVariable[possiblePivotKey_]; // and set incoming here if (sequenceIn >= numberRows + numberColumns) { // slack in - so use old key sequenceIn = keyVariable_[iSetIn]; model->setStatus(sequenceIn, ClpSimplex::basic); setStatus(iSetIn, ClpSimplex::basic); redoSet(model, iSetIn + numberColumns, keyVariable_[iSetIn], iSetIn); } //? do not do if iSetIn<0 ? as will be done later pivotVariable[possiblePivotKey_] = sequenceIn; if (sequenceIn < numberColumns) backToPivotRow_[sequenceIn] = possiblePivotKey_; redoSet(model, key, keyVariable_[iSetOut], iSetOut); } } else { if (sequenceOut < numberColumns) { if (iSetIn >= 0 && iSetOut == iSetIn) { // key not out - only problem is if slack in int key; if (sequenceIn >= numberRows + numberColumns) { key = numberColumns + iSetIn; setStatus(iSetIn, ClpSimplex::basic); assert(pivotRow < numberRows); // must swap with current key int key = keyVariable_[iSetIn]; model->setStatus(key, ClpSimplex::basic); pivotVariable[pivotRow] = key; backToPivotRow_[key] = pivotRow; } else { key = keyVariable_[iSetIn]; } redoSet(model, key, keyVariable_[iSetIn], iSetIn); } else if (iSetOut >= 0) { // just redo set int key = keyVariable_[iSetOut]; ; redoSet(model, key, keyVariable_[iSetOut], iSetOut); } } } if (iSetIn >= 0 && iSetIn != iSetOut) { int key = keyVariable_[iSetIn]; if (sequenceIn == numberColumns + 2 * numberRows) { // key slack in assert(pivotRow < numberRows); // must swap with current key model->setStatus(key, ClpSimplex::basic); pivotVariable[pivotRow] = key; backToPivotRow_[key] = pivotRow; setStatus(iSetIn, ClpSimplex::basic); key = iSetIn + numberColumns; } // redo set to allow for new one redoSet(model, key, keyVariable_[iSetIn], iSetIn); } // update pivot if (sequenceIn < numberColumns) { if (pivotRow < numberRows) { backToPivotRow_[sequenceIn] = pivotRow; } else { if (possiblePivotKey_ >= 0) { assert(possiblePivotKey_ < numberRows); backToPivotRow_[sequenceIn] = possiblePivotKey_; pivotVariable[possiblePivotKey_] = sequenceIn; } } } else if (sequenceIn >= numberRows + numberColumns) { // key in - something should have been done before int key = keyVariable_[iSetIn]; assert(key == numberColumns + iSetIn); //pivotVariable[pivotRow]=key; //backToPivotRow_[key]=pivotRow; //model->setStatus(key,ClpSimplex::basic); //key=numberColumns+iSetIn; setStatus(iSetIn, ClpSimplex::basic); redoSet(model, key, keyVariable_[iSetIn], iSetIn); } #ifdef CLP_DEBUG { char *xx = new char[numberColumns + numberRows]; memset(xx, 0, numberRows + numberColumns); for (int i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; assert(iPivot < numberRows + numberColumns); assert(!xx[iPivot]); xx[iPivot] = 1; if (iPivot < numberColumns) { int iBack = backToPivotRow_[iPivot]; assert(i == iBack); } } delete[] xx; } #endif if (rhsOffset_) { // update effective rhs if (sequenceIn != sequenceOut) { if (sequenceIn < numberColumns) { if (iSetIn >= 0) { // new contribution to rhsOffset_ int key = keyVariable_[iSetIn]; if (key < numberColumns) { double newB = 0.0; ClpSimplex::Status iStatus = getStatus(iSetIn); if (iStatus == ClpSimplex::atLowerBound) newB = lower_[iSetIn]; else newB = upper_[iSetIn]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(key + 1); int iColumn = next_[key]; // skip basic while (iColumn >= 0) iColumn = next_[iColumn]; // sum all non-key variables while (iColumn != stop) { assert(iColumn < 0); iColumn = -iColumn - 1; newB -= solution[iColumn]; iColumn = next_[iColumn]; } } if (newB) ClpPackedMatrix::add(model, rhsOffset_, key, -newB); } } } if (iSetOut >= 0) { // new contribution to rhsOffset_ int key = keyVariable_[iSetOut]; if (key < numberColumns && iSetIn != iSetOut) { double newB = 0.0; ClpSimplex::Status iStatus = getStatus(iSetOut); if (iStatus == ClpSimplex::atLowerBound) newB = lower_[iSetOut]; else newB = upper_[iSetOut]; // subtract out others at bounds if ((gubType_ & 8) == 0) { int stop = -(key + 1); int iColumn = next_[key]; // skip basic while (iColumn >= 0) iColumn = next_[iColumn]; // sum all non-key variables while (iColumn != stop) { assert(iColumn < 0); iColumn = -iColumn - 1; newB -= solution[iColumn]; iColumn = next_[iColumn]; } } if (newB) ClpPackedMatrix::add(model, rhsOffset_, key, -newB); } } } } #ifdef CLP_DEBUG // debug { int i; char *xxxx = new char[numberColumns]; memset(xxxx, 0, numberColumns); for (i = 0; i < numberRows; i++) { int iPivot = pivotVariable[i]; assert(model->getStatus(iPivot) == ClpSimplex::basic); if (iPivot < numberColumns && backward_[iPivot] >= 0) xxxx[iPivot] = 1; } double primalTolerance = model->primalTolerance(); for (i = 0; i < numberSets_; i++) { int key = keyVariable_[i]; double value = 0.0; // sum over all except key int iColumn = next_[key]; // sum all non-key variables int k = 0; int stop = -(key + 1); while (iColumn != stop) { if (iColumn < 0) iColumn = -iColumn - 1; value += solution[iColumn]; k++; assert(k < 100); assert(backward_[iColumn] == i); iColumn = next_[iColumn]; } iColumn = next_[key]; if (key < numberColumns) { // feasibility will be done later assert(getStatus(i) != ClpSimplex::basic); double sol; if (getStatus(i) == ClpSimplex::atUpperBound) sol = upper_[i] - value; else sol = lower_[i] - value; //printf("xx Value of key structural %d for set %d is %g - cost %g\n",key,i,sol, // cost[key]); //if (fabs(sol-solution[key])>1.0e-3) //printf("** stored value was %g\n",solution[key]); } else { // slack is key double infeasibility = 0.0; if (value > upper_[i] + primalTolerance) { infeasibility = value - upper_[i] - primalTolerance; //setAbove(i); } else if (value < lower_[i] - primalTolerance) { infeasibility = lower_[i] - value - primalTolerance; //setBelow(i); } else { //setFeasible(i); } //printf("xx Value of key slack for set %d is %g\n",i,value); } while (iColumn >= 0) { assert(xxxx[iColumn]); xxxx[iColumn] = 0; iColumn = next_[iColumn]; } } for (i = 0; i < numberColumns; i++) { if (i < numberColumns && backward_[i] >= 0) { assert(!xxxx[i] || i == keyVariable_[backward_[i]]); } } delete[] xxxx; } #endif return 0; } // Switches off dj checking each factorization (for BIG models) void ClpGubMatrix::switchOffCheck() { noCheck_ = 0; infeasibilityWeight_ = 0.0; } // Correct sequence in and out to give true value void ClpGubMatrix::correctSequence(const ClpSimplex * /*model*/, int &sequenceIn, int &sequenceOut) { if (sequenceIn != -999) { sequenceIn = trueSequenceIn_; sequenceOut = trueSequenceOut_; } } /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 */