// Copyright (C) 2002, International Business Machines // Corporation and others. All Rights Reserved. //#undef NDEBUG #include "ClpConfig.h" #include "CoinPragma.hpp" #include #if SLIM_CLP==2 #define SLIM_NOIO #endif #include "CoinHelperFunctions.hpp" #include "ClpSimplex.hpp" #include "ClpFactorization.hpp" #include "ClpPackedMatrix.hpp" #include "CoinIndexedVector.hpp" #include "ClpDualRowDantzig.hpp" #include "ClpDualRowSteepest.hpp" #include "ClpPrimalColumnDantzig.hpp" #include "ClpPrimalColumnSteepest.hpp" #include "ClpNonLinearCost.hpp" #include "ClpMessage.hpp" #include "ClpEventHandler.hpp" #include "ClpLinearObjective.hpp" #include "ClpHelperFunctions.hpp" #include "CoinModel.hpp" #include "CoinLpIO.hpp" #include #include #include #include //############################################################################# ClpSimplex::ClpSimplex (bool emptyMessages) : ClpModel(emptyMessages), columnPrimalInfeasibility_(0.0), rowPrimalInfeasibility_(0.0), columnPrimalSequence_(-2), rowPrimalSequence_(-2), columnDualInfeasibility_(0.0), rowDualInfeasibility_(0.0), moreSpecialOptions_(2), rowDualSequence_(-2), primalToleranceToGetOptimal_(-1.0), remainingDualInfeasibility_(0.0), largeValue_(1.0e15), largestPrimalError_(0.0), largestDualError_(0.0), alphaAccuracy_(-1.0), dualBound_(1.0e10), alpha_(0.0), theta_(0.0), lowerIn_(0.0), valueIn_(0.0), upperIn_(-COIN_DBL_MAX), dualIn_(0.0), lowerOut_(-1), valueOut_(-1), upperOut_(-1), dualOut_(-1), dualTolerance_(0.0), primalTolerance_(0.0), sumDualInfeasibilities_(0.0), sumPrimalInfeasibilities_(0.0), infeasibilityCost_(1.0e10), sumOfRelaxedDualInfeasibilities_(0.0), sumOfRelaxedPrimalInfeasibilities_(0.0), acceptablePivot_(1.0e-8), lower_(NULL), rowLowerWork_(NULL), columnLowerWork_(NULL), upper_(NULL), rowUpperWork_(NULL), columnUpperWork_(NULL), cost_(NULL), rowObjectiveWork_(NULL), objectiveWork_(NULL), sequenceIn_(-1), directionIn_(-1), sequenceOut_(-1), directionOut_(-1), pivotRow_(-1), lastGoodIteration_(-100), dj_(NULL), rowReducedCost_(NULL), reducedCostWork_(NULL), solution_(NULL), rowActivityWork_(NULL), columnActivityWork_(NULL), auxiliaryModel_(NULL), numberDualInfeasibilities_(0), numberDualInfeasibilitiesWithoutFree_(0), numberPrimalInfeasibilities_(100), numberRefinements_(0), pivotVariable_(NULL), factorization_(NULL), savedSolution_(NULL), numberTimesOptimal_(0), disasterArea_(NULL), changeMade_(1), algorithm_(0), forceFactorization_(-1), perturbation_(100), nonLinearCost_(NULL), lastBadIteration_(-999999), lastFlaggedIteration_(-999999), numberFake_(0), numberChanged_(0), progressFlag_(0), firstFree_(-1), numberExtraRows_(0), maximumBasic_(0), incomingInfeasibility_(1.0), allowedInfeasibility_(10.0), automaticScale_(0), progress_(NULL) { int i; for (i=0;i<6;i++) { rowArray_[i]=NULL; columnArray_[i]=NULL; } for (i=0;i<4;i++) { spareIntArray_[i]=0; spareDoubleArray_[i]=0.0; } saveStatus_=NULL; // get an empty factorization so we can set tolerances etc getEmptyFactorization(); // Say sparse factorization_->sparseThreshold(1); // say Steepest pricing dualRowPivot_ = new ClpDualRowSteepest(); // say Steepest pricing primalColumnPivot_ = new ClpPrimalColumnSteepest(); solveType_=1; // say simplex based life form } // Subproblem constructor ClpSimplex::ClpSimplex ( const ClpModel * rhs, int numberRows, const int * whichRow, int numberColumns, const int * whichColumn, bool dropNames, bool dropIntegers,bool fixOthers) : ClpModel(rhs, numberRows, whichRow, numberColumns,whichColumn,dropNames,dropIntegers), columnPrimalInfeasibility_(0.0), rowPrimalInfeasibility_(0.0), columnPrimalSequence_(-2), rowPrimalSequence_(-2), columnDualInfeasibility_(0.0), rowDualInfeasibility_(0.0), moreSpecialOptions_(2), rowDualSequence_(-2), primalToleranceToGetOptimal_(-1.0), remainingDualInfeasibility_(0.0), largeValue_(1.0e15), largestPrimalError_(0.0), largestDualError_(0.0), alphaAccuracy_(-1.0), dualBound_(1.0e10), alpha_(0.0), theta_(0.0), lowerIn_(0.0), valueIn_(0.0), upperIn_(-COIN_DBL_MAX), dualIn_(0.0), lowerOut_(-1), valueOut_(-1), upperOut_(-1), dualOut_(-1), dualTolerance_(0.0), primalTolerance_(0.0), sumDualInfeasibilities_(0.0), sumPrimalInfeasibilities_(0.0), infeasibilityCost_(1.0e10), sumOfRelaxedDualInfeasibilities_(0.0), sumOfRelaxedPrimalInfeasibilities_(0.0), acceptablePivot_(1.0e-8), lower_(NULL), rowLowerWork_(NULL), columnLowerWork_(NULL), upper_(NULL), rowUpperWork_(NULL), columnUpperWork_(NULL), cost_(NULL), rowObjectiveWork_(NULL), objectiveWork_(NULL), sequenceIn_(-1), directionIn_(-1), sequenceOut_(-1), directionOut_(-1), pivotRow_(-1), lastGoodIteration_(-100), dj_(NULL), rowReducedCost_(NULL), reducedCostWork_(NULL), solution_(NULL), rowActivityWork_(NULL), columnActivityWork_(NULL), auxiliaryModel_(NULL), numberDualInfeasibilities_(0), numberDualInfeasibilitiesWithoutFree_(0), numberPrimalInfeasibilities_(100), numberRefinements_(0), pivotVariable_(NULL), factorization_(NULL), savedSolution_(NULL), numberTimesOptimal_(0), disasterArea_(NULL), changeMade_(1), algorithm_(0), forceFactorization_(-1), perturbation_(100), nonLinearCost_(NULL), lastBadIteration_(-999999), lastFlaggedIteration_(-999999), numberFake_(0), numberChanged_(0), progressFlag_(0), firstFree_(-1), numberExtraRows_(0), maximumBasic_(0), incomingInfeasibility_(1.0), allowedInfeasibility_(10.0), automaticScale_(0), progress_(NULL) { int i; for (i=0;i<6;i++) { rowArray_[i]=NULL; columnArray_[i]=NULL; } for (i=0;i<4;i++) { spareIntArray_[i]=0; spareDoubleArray_[i]=0.0; } saveStatus_=NULL; // get an empty factorization so we can set tolerances etc getEmptyFactorization(); // say Steepest pricing dualRowPivot_ = new ClpDualRowSteepest(); // say Steepest pricing primalColumnPivot_ = new ClpPrimalColumnSteepest(); solveType_=1; // say simplex based life form if (fixOthers) { int numberOtherColumns=rhs->numberColumns(); int numberOtherRows=rhs->numberRows(); double * solution = new double [numberOtherColumns]; CoinZeroN(solution,numberOtherColumns); int i; for (i=0;iprimalColumnSolution(); const double * objective = rhs->objective(); double offset=0.0; for (i=0;imatrix()->times(solution,rhsModification) ; for ( i=0;i-1.0e20) rowLower_[i] -= rhsModification[iRow]; if (rowUpper_[i]<1.0e20) rowUpper_[i] -= rhsModification[iRow]; } delete [] rhsModification; setObjectiveOffset(rhs->objectiveOffset()-offset); // And set objective value to match setObjectiveValue(rhs->objectiveValue()); } delete [] solution; } } // Puts solution back small model void ClpSimplex::getbackSolution(const ClpSimplex & smallModel,const int * whichRow, const int * whichColumn) { setSumDualInfeasibilities(smallModel.sumDualInfeasibilities()); setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities()); setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities()); setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities()); setNumberIterations(smallModel.numberIterations()); setProblemStatus(smallModel.status()); setObjectiveValue(smallModel.objectiveValue()); const double * solution2 = smallModel.primalColumnSolution(); int i; int numberRows2 = smallModel.numberRows(); int numberColumns2 = smallModel.numberColumns(); const double * dj2 = smallModel.dualColumnSolution(); for ( i=0;itimes(columnActivity_,rowActivity_) ; } //----------------------------------------------------------------------------- ClpSimplex::~ClpSimplex () { setPersistenceFlag(0); gutsOfDelete(0); delete nonLinearCost_; } //############################################################################# void ClpSimplex::setLargeValue( double value) { if (value>0.0&&value0); // only primal at present assert(nonLinearCost_); int iRow; checkPrimalSolution( rowActivityWork_, columnActivityWork_); // get correct bounds on all variables nonLinearCost_->checkInfeasibilities(primalTolerance_); oldValue = nonLinearCost_->largestInfeasibility(); save = new double[numberRows_]; for (iRow=0;iRow0&&nonLinearCost_!=NULL) { // primal algorithm // get correct bounds on all variables // If 4 bit set - Force outgoing variables to exact bound (primal) if ((specialOptions_&4)==0) nonLinearCost_->checkInfeasibilities(primalTolerance_); else nonLinearCost_->checkInfeasibilities(0.0); objectiveModification += nonLinearCost_->changeInCost(); if (nonLinearCost_->numberInfeasibilities()) if (handler_->detail(CLP_SIMPLEX_NONLINEAR,messages_)<100) { handler_->message(CLP_SIMPLEX_NONLINEAR,messages_) <changeInCost() <numberInfeasibilities() <largestInfeasibility(); #ifdef CLP_DEBUG std::cout<<"Largest given infeasibility "<largestInfeasibility()< (CoinMax(10.0*allowedInfeasibility_,100.0*oldValue))) &&(badInfeasibility>CoinMax(incomingInfeasibility_,allowedInfeasibility_)|| useError>1.0e-3)) { //printf("Original largest infeas %g, now %g, primalError %g\n", // oldValue,nonLinearCost_->largestInfeasibility(), // largestPrimalError_); // throw out up to 1000 structurals int iRow; int * sort = new int[numberRows_]; // first put back solution and store difference for (iRow=0;iRow1.0e-4) { sort[numberOut]=iPivot; save[numberOut++]=difference; if (getStatus(iPivot)==basic) numberBasic++; } } } if (!numberBasic) { //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut); allSlackBasis(true); } CoinSort_2(save, save + numberOut, sort, CoinFirstGreater_2()); numberOut = CoinMin(1000,numberOut); for (iRow=0;iRow1.0e10) { if (upper_[iColumn]<0.0) { solution_[iColumn]=upper_[iColumn]; } else if (lower_[iColumn]>0.0) { solution_[iColumn]=lower_[iColumn]; } else { solution_[iColumn]=0.0; } } } delete [] sort; } delete [] save; if (numberOut) return numberOut; } computeDuals(givenDuals); // now check solutions //checkPrimalSolution( rowActivityWork_, columnActivityWork_); //checkDualSolution(); checkBothSolutions(); objectiveValue_ += objectiveModification/(objectiveScale_*rhsScale_); if (handler_->logLevel()>3||(largestPrimalError_>1.0e-2|| largestDualError_>1.0e-2)) handler_->message(CLP_SIMPLEX_ACCURACY,messages_) <0) firstFree_ = -1; return 0; } void ClpSimplex::computePrimals ( const double * rowActivities, const double * columnActivities) { //work space CoinIndexedVector * workSpace = rowArray_[0]; CoinIndexedVector * arrayVector = rowArray_[1]; arrayVector->clear(); CoinIndexedVector * previousVector = rowArray_[2]; previousVector->clear(); // accumulate non basic stuff int iRow; // order is this way for scaling if (columnActivities!=columnActivityWork_) ClpDisjointCopyN(columnActivities,numberColumns_,columnActivityWork_); if (rowActivities!=rowActivityWork_) ClpDisjointCopyN(rowActivities,numberRows_,rowActivityWork_); double * array = arrayVector->denseVector(); int * index = arrayVector->getIndices(); int number=0; const double * rhsOffset = matrix_->rhsOffset(this,false,true); if (!rhsOffset) { // Use whole matrix every time to make it easier for ClpMatrixBase // So zero out basic for (iRow=0;iRow=0); solution_[iPivot] = 0.0; } // Extended solution before "update" matrix_->primalExpanded(this,0); times(-1.0,columnActivityWork_,array); for (iRow=0;iRowsetNumElements(number); #ifdef CLP_DEBUG if (numberIterations_==-3840) { int i; for (i=0;iupdateColumn(workSpace,thisVector); double * work = workSpace->denseVector(); #ifdef CLP_DEBUG if (numberIterations_==-3840) { int i; for (i=0;igetNumElements(); int * indexIn = thisVector->getIndices(); double * arrayIn = thisVector->denseVector(); // put solution in correct place if (!rhsOffset) { int j; for (j=0;jprimalExpanded(this,1); // check Ax == b (for all) // signal column generated matrix to just do basic (and gub) unsigned int saveOptions = specialOptions(); setSpecialOptions(16); times(-1.0,columnActivityWork_,work); setSpecialOptions(saveOptions); largestPrimalError_=0.0; double multiplier = 131072.0; for (iRow=0;iRowlargestPrimalError_) { largestPrimalError_=fabs(value); } } if (largestPrimalError_>=lastError) { // restore CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; goodSolution=false; break; } if (iRefine1.0e-10) { // try and make better // save this CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; int * indexOut = thisVector->getIndices(); int number=0; array = thisVector->denseVector(); thisVector->clear(); for (iRow=0;iRowsetNumElements(number); lastError=largestPrimalError_; factorization_->updateColumn(workSpace,thisVector); multiplier = 1.0/multiplier; double * previous = lastVector->denseVector(); number=0; for (iRow=0;iRowsetNumElements(number); } else { break; } } // solution as accurate as we are going to get ClpFillN(work,numberRows_,0.0); if (!goodSolution) { array = thisVector->denseVector(); // put solution in correct place for (iRow=0;iRowclear(); previousVector->clear(); #ifdef CLP_DEBUG if (numberIterations_==-3840) { exit(77); } #endif } // now dual side void ClpSimplex::computeDuals(double * givenDjs) { #ifndef SLIM_CLP if (objective_->type()==1||!objective_->activated()) { #endif // Linear //work space CoinIndexedVector * workSpace = rowArray_[0]; CoinIndexedVector * arrayVector = rowArray_[1]; arrayVector->clear(); CoinIndexedVector * previousVector = rowArray_[2]; previousVector->clear(); int iRow; #ifdef CLP_DEBUG workSpace->checkClear(); #endif double * array = arrayVector->denseVector(); int * index = arrayVector->getIndices(); int number=0; if (!givenDjs) { for (iRow=0;iRowsetNumElements(number); // Extended duals before "updateTranspose" matrix_->dualExpanded(this,arrayVector,givenDjs,0); // Btran basic costs and get as accurate as possible double lastError=COIN_DBL_MAX; int iRefine; double * work = workSpace->denseVector(); CoinIndexedVector * thisVector = arrayVector; CoinIndexedVector * lastVector = previousVector; factorization_->updateColumnTranspose(workSpace,thisVector); for (iRefine=0;iRefinegetIndices(); // use reduced costs for slacks as work array double * work2 = reducedCostWork_+numberColumns_; int numberStructurals=0; for (iRow=0;iRowlistTransposeTimes(this,array,index2,numberStructurals,work2); numberStructurals=0; if (!givenDjs) { for (iRow=0;iRow=numberColumns_) { // slack value = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]; } else { // column value = objectiveWork_[iPivot]-work2[numberStructurals++]; } work[iRow]=value; if (fabs(value)>largestDualError_) { largestDualError_=fabs(value); } } } else { for (iRow=0;iRow=numberColumns_) { // slack work[iRow] = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]-givenDjs[iPivot]; } else { // column work[iRow] = objectiveWork_[iPivot]-work2[numberStructurals++] - givenDjs[iPivot]; } if (fabs(work[iRow])>largestDualError_) { largestDualError_=fabs(work[iRow]); //assert (largestDualError_<1.0e-7); //if (largestDualError_>1.0e-7) //printf("large dual error %g\n",largestDualError_); } } } } else { // extra rows - be more careful #if 1 // would be faster to do just for basic but this reduces code ClpDisjointCopyN(objectiveWork_,numberColumns_,reducedCostWork_); transposeTimes(-1.0,array,reducedCostWork_); #else // Just basic int * index2 = workSpace->getIndices(); int numberStructurals=0; for (iRow=0;iRowlistTransposeTimes(this,array,index2,numberStructurals,work); for (iRow=0;iRowdualExpanded(this,NULL,NULL,1); if (!givenDjs) { for (iRow=0;iRow=numberColumns_) { // slack value = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]; } else { // column value = reducedCostWork_[iPivot]; } work[iRow]=value; if (fabs(value)>largestDualError_) { largestDualError_=fabs(value); } } } else { for (iRow=0;iRow=numberColumns_) { // slack work[iRow] = rowObjectiveWork_[iPivot-numberColumns_] + array[iPivot-numberColumns_]-givenDjs[iPivot]; } else { // column work[iRow] = reducedCostWork_[iPivot]- givenDjs[iPivot]; } if (fabs(work[iRow])>largestDualError_) { largestDualError_=fabs(work[iRow]); //assert (largestDualError_<1.0e-7); //if (largestDualError_>1.0e-7) //printf("large dual error %g\n",largestDualError_); } } } } if (largestDualError_>=lastError) { // restore CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; break; } if (iRefine1.0e-10 &&!givenDjs) { // try and make better // save this CoinIndexedVector * temp = thisVector; thisVector = lastVector; lastVector=temp; int * indexOut = thisVector->getIndices(); int number=0; array = thisVector->denseVector(); thisVector->clear(); double multiplier = 131072.0; for (iRow=0;iRowsetNumElements(number); lastError=largestDualError_; factorization_->updateColumnTranspose(workSpace,thisVector); multiplier = 1.0/multiplier; double * previous = lastVector->denseVector(); number=0; for (iRow=0;iRowsetNumElements(number); } else { break; } } // now look at dual solution array = thisVector->denseVector(); for (iRow=0;iRow10000) matrix_->transposeTimes(-1.0,dual_,reducedCostWork_, rowScale_,columnScale_,work); else matrix_->transposeTimes(-1.0,dual_,reducedCostWork_, rowScale_,columnScale_,NULL); ClpFillN(work,numberRows_,0.0); // Extended duals and check dual infeasibility if (!matrix_->skipDualCheck()||algorithm_<0||problemStatus_!=-2) matrix_->dualExpanded(this,NULL,NULL,2); // If necessary - override results if (givenDjs) { // restore accurate duals memcpy(givenDjs,dj_,(numberRows_+numberColumns_)*sizeof(double)); } arrayVector->clear(); previousVector->clear(); #ifndef SLIM_CLP } else { // Nonlinear objective_->reducedGradient(this,dj_,false); // get dual_ by moving from reduced costs for slacks memcpy(dual_,dj_+numberColumns_,numberRows_*sizeof(double)); } #endif } /* Given an existing factorization computes and checks primal and dual solutions. Uses input arrays for variables at bounds. Returns feasibility states */ int ClpSimplex::getSolution ( const double * rowActivities, const double * columnActivities) { if (!factorization_->status()) { // put in standard form createRim(7+8+16+32,false,-1); if (pivotVariable_[0]<0) internalFactorize(0); // do work gutsOfSolution ( NULL,NULL); // release extra memory deleteRim(0); } return factorization_->status(); } /* Given an existing factorization computes and checks primal and dual solutions. Uses current problem arrays for bounds. Returns feasibility states */ int ClpSimplex::getSolution ( ) { double * rowActivities = new double[numberRows_]; double * columnActivities = new double[numberColumns_]; ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities); ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities); int status = getSolution( rowActivities, columnActivities); delete [] rowActivities; delete [] columnActivities; return status; } // Factorizes using current basis. This is for external use // Return codes are as from ClpFactorization int ClpSimplex::factorize () { // put in standard form createRim(7+8+16+32,false); // do work int status = internalFactorize(-1); // release extra memory deleteRim(0); return status; } // Clean up status void ClpSimplex::cleanStatus() { int iRow,iColumn; int numberBasic=0; // make row activities correct memset(rowActivityWork_,0,numberRows_*sizeof(double)); times(1.0,columnActivityWork_,rowActivityWork_); if (!status_) createStatus(); for (iRow=0;iRow all slack basis. */ int ClpSimplex::internalFactorize ( int solveType) { int iRow,iColumn; int totalSlacks=numberRows_; if (!status_) createStatus(); bool valuesPass=false; if (solveType>=10) { valuesPass=true; solveType -= 10; } #ifdef CLP_DEBUG if (solveType>0) { int numberFreeIn=0,numberFreeOut=0; double biggestDj=0.0; for (iColumn=0;iColumnlargeValue_) numberFreeIn++; break; default: if (columnLower_[iColumn]<-largeValue_ &&columnUpper_[iColumn]>largeValue_) { numberFreeOut++; biggestDj = CoinMax(fabs(dj_[iColumn]),biggestDj); } break; } } if (numberFreeIn+numberFreeOut) printf("%d in basis, %d out - largest dj %g\n", numberFreeIn,numberFreeOut,biggestDj); } #endif if (solveType<=0) { // Make sure everything is clean for (iRow=0;iRowrowLowerWork_[iRow]) setRowStatus(iRow,atLowerBound); } else if (getRowStatus(iRow)==isFree) { // may not be free after all if (rowLowerWork_[iRow]>-largeValue_||rowUpperWork_[iRow]columnLowerWork_[iColumn]) setColumnStatus(iColumn,atLowerBound); } else if (getColumnStatus(iColumn)==isFree) { // may not be free after all if (columnLowerWork_[iColumn]>-largeValue_||columnUpperWork_[iColumn]largeValue_) { if (rowLowerWork_[iRow]>-largeValue_) { rowActivityWork_[iRow]=rowLowerWork_[iRow]; setRowStatus(iRow,atLowerBound); } else { // say free setRowStatus(iRow,isFree); rowActivityWork_[iRow]=0.0; } } break; case ClpSimplex::isFixed: case atLowerBound: rowActivityWork_[iRow]=rowLowerWork_[iRow]; if (rowActivityWork_[iRow]<-largeValue_) { if (rowUpperWork_[iRow]largeValue_) { if (rowLowerWork_[iRow]>-largeValue_) { rowActivityWork_[iRow]=rowLowerWork_[iRow]; setRowStatus(iRow,atLowerBound); } else { // say free setRowStatus(iRow,isFree); rowActivityWork_[iRow]=0.0; } } else { if (rowLowerWork_[iRow]>-largeValue_) { // set to nearest if (fabs(rowActivityWork_[iRow]-rowLowerWork_[iRow]) -largeValue_) { if (columnActivityWork_[iColumn]-columnLowerWork_[iColumn]< columnUpperWork_[iColumn]-columnActivityWork_[iColumn]) { columnActivityWork_[iColumn]=columnLowerWork_[iColumn]; setColumnStatus(iColumn,atLowerBound); } else { columnActivityWork_[iColumn]=columnUpperWork_[iColumn]; setColumnStatus(iColumn,atUpperBound); } } else if (columnUpperWork_[iColumn]largeValue_) { if (columnLowerWork_[iColumn]<-largeValue_) { columnActivityWork_[iColumn]=0.0; setColumnStatus(iColumn,isFree); } else { columnActivityWork_[iColumn]=columnLowerWork_[iColumn]; setColumnStatus(iColumn,atLowerBound); } } break; case isFixed: case atLowerBound: columnActivityWork_[iColumn]=columnLowerWork_[iColumn]; if (columnActivityWork_[iColumn]<-largeValue_) { if (columnUpperWork_[iColumn]>largeValue_) { columnActivityWork_[iColumn]=0.0; setColumnStatus(iColumn,isFree); } else { columnActivityWork_[iColumn]=columnUpperWork_[iColumn]; setColumnStatus(iColumn,atUpperBound); } } break; case isFree: break; // not really free - fall through to superbasic case superBasic: if (columnUpperWork_[iColumn]>largeValue_) { if (columnLowerWork_[iColumn]>-largeValue_) { columnActivityWork_[iColumn]=columnLowerWork_[iColumn]; setColumnStatus(iColumn,atLowerBound); } else { // say free setColumnStatus(iColumn,isFree); columnActivityWork_[iColumn]=0.0; } } else { if (columnLowerWork_[iColumn]>-largeValue_) { // set to nearest if (fabs(columnActivityWork_[iColumn]-columnLowerWork_[iColumn]) -largeValue_||upper-big_bound||upper1.0e10) { // printf("large %g at %d - status %d\n", // solution_[iRow],iRow,status_[iRow]); //} //} if (0) { static int k=0; printf("start basis\n"); int i; for (i=0;i20) exit(0); k++; } int status = factorization_->factorize(this, solveType,valuesPass); if (status) { handler_->message(CLP_SIMPLEX_BADFACTOR,messages_) <sparseThreshold()) { // get default value factorization_->sparseThreshold(0); factorization_->goSparse(); //} return status; } /* This does basis housekeeping and does values for in/out variables. Can also decide to re-factorize */ int ClpSimplex::housekeeping(double objectiveChange) { // save value of incoming and outgoing double oldIn = solution_[sequenceIn_]; double oldOut = solution_[sequenceOut_]; numberIterations_++; changeMade_++; // something has happened // incoming variable if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) { handler_->message(CLP_SIMPLEX_HOUSE1,messages_) <message(CLP_SIMPLEX_FREEIN,messages_) <=0) pivotVariable_[pivotRow_]=sequenceIn(); if (upper_[sequenceIn_]>1.0e20&&lower_[sequenceIn_]<-1.0e20) progressFlag_ |= 2; // making real progress solution_[sequenceIn_]=valueIn_; if (upper_[sequenceOut_]-lower_[sequenceOut_]<1.0e-12) progressFlag_ |= 1; // making real progress if (sequenceIn_!=sequenceOut_) { if (alphaAccuracy_>0.0) { double value = fabs(alpha_); if (value>1.0) alphaAccuracy_ *= value; else alphaAccuracy_ /= value; } //assert( getStatus(sequenceOut_)== basic); setStatus(sequenceIn_,basic); if (upper_[sequenceOut_]-lower_[sequenceOut_]>0) { // As Nonlinear costs may have moved bounds (to more feasible) // Redo using value if (fabs(valueOut_-lower_[sequenceOut_])type()<2) //assert (fabs(theta_)>1.0e-13); // flip from bound to bound // As Nonlinear costs may have moved bounds (to more feasible) // Redo using value if (fabs(valueIn_-lower_[sequenceIn_])