// Copyright (C) 2002, International Business Machines // Corporation and others. All Rights Reserved. #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=0) longest[iSet]++; } int length = 0; for (j=0;jmatrix()), 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) 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]=0) { firstGub_ = CoinMin(firstGub_,i); lastGub_ = CoinMax(lastGub_,i); } } gubType_=0; // adjust lastGub_ if (lastGub_>0) lastGub_++; for (i=firstGub_;i=0) longest[iSet]++; } int length = 0; for (j=0;jgetNumCols(); int * array = new int [ numberColumnsOld]; int i; for (i=0;i=0) { firstGub_ = CoinMin(firstGub_,i); lastGub_ = CoinMax(lastGub_,i); } } if (lastGub_>0) lastGub_++; gubType_=0; for (i=firstGub_;ix * 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->factorization()->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*10numberIterations()%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=0) { int iBasic = keyVariable_[iSet]; if (iBasicgetStatus(iBasic)==ClpSimplex::basic); djMod=0.0; for (CoinBigIndex j=columnStart[iBasic]; jzeroTolerance) { array[numberNonZero]=value; index[numberNonZero++]=iColumn; } } } else { // scaled // modify pi so can collapse to one loop for (i=0;i=0) { int iBasic = keyVariable_[iSet]; if (iBasicgetStatus(iBasic)==ClpSimplex::basic); djMod=0.0; // scaled for (CoinBigIndex j=columnStart[iBasic]; jcolumnScale(); for (j=columnStart[iColumn]; jzeroTolerance) { array[numberNonZero]=value; index[numberNonZero++]=iColumn; } } } // zero out for (i=0;izeroTolerance) { index[numberNonZero++]=iColumn; array[iColumn]=-value; } } } else if (scalar==1.0) { for (iColumn=0;iColumnzeroTolerance) { index[numberNonZero++]=iColumn; array[iColumn]=value; } } } else { for (iColumn=0;iColumnzeroTolerance) { index[numberNonZero++]=iColumn; array[iColumn]=value; } } } } else { // scaled if (scalar==-1.0) { for (iColumn=0;iColumncolumnScale(); for (j=columnStart[iColumn]; jzeroTolerance) { index[numberNonZero++]=iColumn; array[iColumn]=-value; } } } else if (scalar==1.0) { for (iColumn=0;iColumncolumnScale(); for (j=columnStart[iColumn]; jzeroTolerance) { index[numberNonZero++]=iColumn; array[iColumn]=value; } } } else { for (iColumn=0;iColumncolumnScale(); for (j=columnStart[iColumn]; jzeroTolerance) { 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;ix * 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=0) { int iBasic = keyVariable_[iSet]; if (iBasic==iColumn) { toIndex_[iSet]=jColumn; fromIndex_[numberTouched++]=iSet; } } } } } else { // scaled for (jColumn=0;jColumncolumnScale(); for (j=columnStart[iColumn]; j=0) { int iBasic = keyVariable_[iSet]; if (iBasic==iColumn) { toIndex_[iSet]=jColumn; fromIndex_[numberTouched++]=iSet; } } } } } // adjust djs for (jColumn=0;jColumn=0) { int kColumn = toIndex_[iSet]; if (kColumn>=0) array[jColumn] -= array[kColumn]; } } // and clear basic for (int j=0;jgetVectorLengths(); int numberRows = getNumRows(); int saveNumberBasic=numberBasic; 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=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=0) { for (j=columnStart[key];j1.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-saveNumberBasic; return numberElements; } void ClpGubMatrix::fillBasis(ClpSimplex * model, const int * whichColumn, int & numberColumnBasic, int * indexRowU, int * start, int * rowCount, int * columnCount, double * 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=numberColumns) { for (j=columnStart[iColumn];j1.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];j1.0e-20) { indexRowU[numberElements]=iRow; rowCount[iRow]++; elementU[numberElements++]=value; } } for (j=columnStart[key];j1.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=numberColumns) { double scale = columnScale[iColumn]; for (j=columnStart[iColumn];j1.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=0) { for (j=columnStart[key];j1.0e-20) { indexRowU[numberElements]=iRow; rowCount[iRow]++; elementU[numberElements++]=value; } } for (j=columnStart[key];j1.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 (iColumnnumberColumns()); // Do packed part ClpPackedMatrix::unpack(model,rowArray,iColumn); int iSet = backward_[iColumn]; if (iSet>=0) { int iBasic = keyVariable_[iSet]; if (iBasic 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=0) { // columns are in order int iBasic = keyVariable_[iSet]; if (iBasic 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]; inext) { lastIndex++; if (lastIndex==numberOld) next=matrix_->getNumRows(); else next=index[lastIndex]; } if (iRowcolumnScale()[iBasic]; for (i=columnStart[iBasic]; inext) { lastIndex++; if (lastIndex==numberOld) next=matrix_->getNumRows(); else next=index[lastIndex]; } if (iRowsetNumElements(number); } } } else { // key slack entering int iBasic = keyVariable_[gubSlackIn_]; assert (iBasic 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]; icolumnScale()[iBasic]; for (i=columnStart[iBasic]; isetNumElements(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 (iColumnnumberColumns()); // 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 (iColumnnumberColumns()); // Do packed part ClpPackedMatrix::add(model,array,iColumn,multiplier); if (iColumnnumberColumns()) { int iSet = backward_[iColumn]; if (iSet>=0&&iColumn!=keyVariable_[iSet]&&keyVariable_[iSet]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 = ((double) firstGub_)/((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_+ (int) (startFraction* (lastGub_-firstGub_)); int endG = firstGub_+ (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;iSequenceminSet) { // 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]; jtolerance) { 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]; jFREE_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]; jtolerance) { 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]; jtolerance) { 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 (bestSequencesetStatus(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;iSequenceminSet) { // 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]; jtolerance) { 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]; jFREE_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]; jtolerance) { 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]; jtolerance) { 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 (bestSequencesetStatus(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 = ((double) lastGub_)/((double) numberColumns); double ratio = ((double) numberColumns)/((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 (sequenceInlowerRegion(); 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=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=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=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=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 (solsaveNumber_) { // clear double theta = model->theta(); double * solution = model->solutionRegion(); for (i=saveNumber_;itheta()*array[i],model->theta()); #endif double value = array[i]; array[i]=0.0; int iSet=fromIndex_[i-saveNumber_]; toIndex_[iSet]=-1; if (iSet==iSetIn&&iColumnsetNumElements(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;isolutionRegion(); ClpSimplex::Status iStatus; //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=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 (kColumnsetStatus(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 iStatus = getStatus(i); assert (iStatus==ClpSimplex::basic); double infeasibility=0.0; if (value>upper_[i]+primalTolerance) { infeasibility=value-upper_[i]-primalTolerance; setAbove(i); } else if (value0.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(); ClpSimplex::Status iStatus; // 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=0) { int kColumn = keyVariable_[iSet]; double costValue; if (kColumnadd(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=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=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;igetStatus(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 (keynonLinearCost()->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=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=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=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 * model, int mode) { 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;igetColumnStatus(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;isetColumnStatus(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=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=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;ifabs(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=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 (keycostRegion(); 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&&iRowrowArray(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=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); memcpy(saveStatus_,status_,numberSets_); memcpy(savedKeyVariable_,keyVariable_,numberSets_*sizeof(int)); } break; // restore status case 6: { memcpy(status_,saveStatus_,numberSets_); memcpy(keyVariable_,savedKeyVariable_,numberSets_*sizeof(int)); // restore firstAvailable_ synchronize(model,7); // redo next_ int i; int * last = new int[numberSets_]; for (i=0;i=numberColumns||backward_[iKey]==i); last[i]=iKey; // make sure basic //if (iKeysetStatus(iKey,ClpSimplex::basic); } for (i=0;i=0) { next_[last[iSet]]=i; last[iSet]=i; } } for (i=0;ipivotVariable(); int iRow; int numberBasic=0; int numberRows = model->numberRows(); for (iRow=0;iRowgetRowStatus(iRow)==ClpSimplex::basic) { numberBasic++; pivotVariable[iRow]=iRow+numberColumns; } else { pivotVariable[iRow]=-1; } } i=0; int iColumn; for (iColumn=0;iColumngetStatus(iColumn)==ClpSimplex::basic) { int iSet = backward_[iColumn]; if (iSet<0||keyVariable_[iSet]!=iColumn) { while (pivotVariable[i]>=0) { i++; assert (isequenceIn()); 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;iSetgetVectorLengths(); 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;iSetprimalTolerance(); bool noNormalBounds=true; gubType_ &= ~8; bool gotBasis=false; for (iSet=0;iSet-1.0e20) noNormalBounds=false; if (columnUpper[j]&&columnUpper[j]<1.0e20) noNormalBounds=false; } } if (noNormalBounds) gubType_ |= 8; if (!gotBasis) { for (iSet=0;iSetgetStatus(j)== ClpSimplex::basic) { if (columnLength[j]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-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)-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)=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 (valuedualTolerance(); 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; double distance=0.0; if (iphase||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; distance = upper[i]-solution[i]; } else { way = -1; distance = solution[i]-lower[i]; } } } 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]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=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;jgetStatus(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;iColumngetStatus(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=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 (solnumberColumns(); 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;jgetStatus(iColumn)==ClpSimplex::basic) { next_[last]=iColumn; next_[iColumn]=lastMarker; last = iColumn; } } } // now add in non-basic for ( j=0;jgetStatus(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 check) { //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;iColumngetColumnStatus(iColumn)==ClpSimplex::basic) solution[iColumn]=0.0; } for (int iSet=0;iSetsolutionRegion(); // and now subtract out non basic ClpSimplex::Status iStatus; for (int iSet=0;iSet=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;iRow1.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;iColumngetColumnStatus(iColumn)==ClpSimplex::basic) solution[iColumn]=0.0; } int iSet; for ( iSet=0;iSetnumberIterations(); const double * columnSolution = model->solutionRegion(); // and now subtract out non basic ClpSimplex::Status iStatus; for ( iSet=0;iSet=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=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=0) { // old contribution to rhsOffset_ int key = keyVariable_[iSetIn]; if (key=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=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=0) { // old contribution to rhsOffset_ int key = keyVariable_[iSetOut]; if (key=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); } } } 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+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=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 (pivotRowsetStatus(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 (pivotRowsetStatus(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=0) { assert (possiblePivotKey_=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=0) { // new contribution to rhsOffset_ int key = keyVariable_[iSetIn]; if (key=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=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;igetStatus(iPivot)==ClpSimplex::basic); if (iPivot=0) xxxx[iPivot]=1; } double primalTolerance = model->primalTolerance(); for (i=0;i1.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=0) { assert (xxxx[iColumn]); xxxx[iColumn]=0; iColumn=next_[iColumn]; } } for (i=0;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_; } }