/* $Id$ */ // Copyright (C) 2002, International Business Machines // Corporation and others, Copyright (C) 2012, FasterCoin. All Rights Reserved. // This code is licensed under the terms of the Eclipse Public License (EPL). /* Notes on implementation of dual simplex algorithm. When dual feasible: If primal feasible, we are optimal. Otherwise choose an infeasible basic variable to leave basis (normally going to nearest bound) (B). We now need to find an incoming variable which will leave problem dual feasible so we get the row of the tableau corresponding to the basic variable (with the correct sign depending if basic variable above or below feasibility region - as that affects whether reduced cost on outgoing variable has to be positive or negative). We now perform a ratio test to determine which incoming variable will preserve dual feasibility (C). If no variable found then problem is infeasible (in primal sense). If there is a variable, we then perform pivot and repeat. Trivial? ------------------------------------------- A) How do we get dual feasible? If all variables have bounds then it is trivial to get feasible by putting non-basic variables to correct bounds. OSL did not have a phase 1/phase 2 approach but instead effectively put fake bounds on variables and this is the approach here, although I had hoped to make it cleaner. If there is a weight of X on getting dual feasible: Non-basic variables with negative reduced costs are put to lesser of their upper bound and their lower bound + X. Similarly, mutatis mutandis, for positive reduced costs. Free variables should normally be in basis, otherwise I have coding which may be able to come out (and may not be correct). In OSL, this weight was changed heuristically, here at present it is only increased if problem looks finished. If problem is feasible I check for unboundedness. If not unbounded we could play with going into primal. As long as weights increase any algorithm would be finite. B) Which outgoing variable to choose is a virtual base class. For difficult problems steepest edge is preferred while for very easy (large) problems we will need partial scan. C) Sounds easy, but this is hardest part of algorithm. 1) Instead of stopping at first choice, we may be able to flip that variable to other bound and if objective still improving choose again. These mini iterations can increase speed by orders of magnitude but we may need to go to more of a bucket choice of variable rather than looking at them one by one (for speed). 2) Accuracy. Reduced costs may be of wrong sign but less than tolerance. Pivoting on these makes objective go backwards. OSL modified cost so a zero move was made, Gill et al (in primal analogue) modified so a strictly positive move was made. It is not quite as neat in dual but that is what we try and do. The two problems are that re-factorizations can change reduced costs above and below tolerances and that when finished we need to reset costs and try again. 3) Degeneracy. Gill et al helps but may not be enough. We may need more. Also it can improve speed a lot if we perturb the costs significantly. References: Forrest and Goldfarb, Steepest-edge simplex algorithms for linear programming - Mathematical Programming 1992 Forrest and Tomlin, Implementing the simplex method for the Optimization Subroutine Library - IBM Systems Journal 1992 Gill, Murray, Saunders, Wright A Practical Anti-Cycling Procedure for Linear and Nonlinear Programming SOL report 1988 TODO: a) Better recovery procedures. At present I never check on forward progress. There is checkpoint/restart with reducing re-factorization frequency, but this is only on singular factorizations. b) Fast methods for large easy problems (and also the option for the code to automatically choose which method). c) We need to be able to stop in various ways for OSI - this is fairly easy. */ #include "CoinPragma.hpp" #include #include "CoinHelperFunctions.hpp" #include "CoinAbcHelperFunctions.hpp" #include "AbcSimplexDual.hpp" #include "ClpEventHandler.hpp" #include "AbcSimplexFactorization.hpp" #include "AbcNonLinearCost.hpp" #include "CoinPackedMatrix.hpp" #include "CoinIndexedVector.hpp" #include "CoinFloatEqual.hpp" #include "AbcDualRowDantzig.hpp" #include "AbcDualRowSteepest.hpp" #include "ClpMessage.hpp" #include "ClpLinearObjective.hpp" class ClpSimplex; #include #include #include #include #include //#define CLP_DEBUG 1 #ifdef NDEBUG #define NDEBUG_CLP #endif #ifndef CLP_INVESTIGATE #define NDEBUG_CLP #endif // dual /* *** Method This is a vanilla version of dual simplex. It tries to be a single phase approach with a weight of 1.0 being given to getting optimal and a weight of dualBound_ being given to getting dual feasible. In this version I have used the idea that this weight can be thought of as a fake bound. If the distance between the lower and upper bounds on a variable is less than the feasibility weight then we are always better off flipping to other bound to make dual feasible. If the distance is greater then we make up a fake bound dualBound_ away from one bound. If we end up optimal or primal infeasible, we check to see if bounds okay. If so we have finished, if not we increase dualBound_ and continue (after checking if unbounded). I am undecided about free variables - there is coding but I am not sure about it. At present I put them in basis anyway. The code is designed to take advantage of sparsity so arrays are seldom zeroed out from scratch or gone over in their entirety. The only exception is a full scan to find outgoing variable. This will be changed to keep an updated list of infeasibilities (or squares if steepest edge). Also on easy problems we don't need full scan - just pick first reasonable. One problem is how to tackle degeneracy and accuracy. At present I am using the modification of costs which I put in OSL and which was extended by Gill et al. I am still not sure of the exact details. The flow of dual is three while loops as follows: while (not finished) { while (not clean solution) { Factorize and/or clean up solution by flipping variables so dual feasible. If looks finished check fake dual bounds. Repeat until status is iterating (-1) or finished (0,1,2) } while (status==-1) { Iterate until no pivot in or out or time to re-factorize. Flow is: choose pivot row (outgoing variable). if none then we are primal feasible so looks as if done but we need to break and check bounds etc. Get pivot row in tableau Choose incoming column. If we don't find one then we look primal infeasible so break and check bounds etc. (Also the pivot tolerance is larger after any iterations so that may be reason) If we do find incoming column, we may have to adjust costs to keep going forwards (anti-degeneracy). Check pivot will be stable and if unstable throw away iteration (we will need to implement flagging of basic variables sometime) and break to re-factorize. If minor error re-factorize after iteration. Update everything (this may involve flipping variables to stay dual feasible. } } At present we never check we are going forwards. I overdid that in OSL so will try and make a last resort. Needs partial scan pivot out option. Needs dantzig, uninitialized and full steepest edge options (can still use partial scan) May need other anti-degeneracy measures, especially if we try and use loose tolerances as a way to solve in fewer iterations. I like idea of dynamic scaling. This gives opportunity to decouple different implications of scaling for accuracy, iteration count and feasibility tolerance. */ #define CLEAN_FIXED 0 // Startup part of dual (may be extended to other algorithms) // To force to follow another run put logfile name here and define //#define FORCE_FOLLOW #ifdef FORCE_FOLLOW static FILE *fpFollow = NULL; static const char *forceFile = NULL; static int force_in = -1; static int force_out = -1; static int force_way_in = -1; static int force_way_out = -1; static int force_iterations = 0; int force_argc = 0; const char **force_argv = NULL; #endif void AbcSimplexDual::startupSolve() { #ifdef FORCE_FOLLOW int ifld; for (ifld = 1; ifld < force_argc; ifld++) { if (strstr(argv[ifld], ".log")) { forceFile = argv[ifld]; break; } } if (!fpFollow && strlen(forceFile)) { fpFollow = fopen(forceFile, "r"); assert(fpFollow); } if (fpFollow) { int numberRead = atoi(argv[ifld + 1]); force_iterations = atoi(argv[ifld + 1]); printf("skipping %d iterations and then doing %d\n", numberRead, force_iterations); for (int iteration = 0; iteration < numberRead; iteration++) { // read to next Clp0102 char temp[300]; while (fgets(temp, 250, fpFollow)) { if (!strncmp(temp, "dirOut", 6)) break; } sscanf(temp + 7, "%d%*s%d%*s%*s%d", &force_way_out, &force_way_in); fgets(temp, 250, fpFollow); char cin, cout; int whichIteration; sscanf(temp, "%d%*f%*s%*c%c%d%*s%*c%c%d", &whichIteration, &cin, &force_in, &cout, &force_out); assert(whichIterations == iteration + 1); if (cin == 'C') force_in += numberColumns_; if (cout == 'C') force_out += numberColumns_; printf("Iteration %d will force %d out (way %d) and %d in (way %d)\n", whichIteration, force_out, force_way_out, force_in, force_way_in); assert(getInternalStatus(force_out) == basic); assert(getInternalStatus(force_in) != basic); setInternalStatus(force_in, basic); if (force_way_out == -1) setInternalStatus(force_out, atUpperBound); else setInternalStatus(force_out, atLowerBound); } // clean up int numberBasic = 0; for (int i = 0; i < numberTotal_; i++) { if (getInternalStatus(i) == basic) { abcPivotVariable_[numberBasic++] = i; } else if (getInternalStatus(i) == atLowerBound) { if (abcLower_[i] < -1.0e30) { abcLower_[i] = abcUpper_[i] - dualBound_; setFakeLower(i); } } else if (getInternalStatus(i) == atUpperBound) { if (abcUpper_[i] > 1.0e30) { abcUpper_[i] = abcLower_[i] + dualBound_; setFakeUpper(i); } } } assert(numberBasic == numberRows_); } #endif // initialize - no values pass and algorithm_ is -1 // put in standard form (and make row copy) // create modifiable copies of model rim and do optional scaling // If problem looks okay // Do initial factorization numberFake_ = 0; // Number of variables at fake bounds numberChanged_ = 0; // Number of variables with changed costs /// Initial coding here statusOfProblemInDual(0); } void AbcSimplexDual::finishSolve() { assert(problemStatus_ || !sumPrimalInfeasibilities_); } void AbcSimplexDual::gutsOfDual() { double largestPrimalError = 0.0; double largestDualError = 0.0; lastCleaned_ = 0; // last time objective or bounds cleaned up numberDisasters_ = 0; // This says whether to restore things etc //startupSolve(); // startup will have factorized so can skip // Start check for cycles abcProgress_.startCheck(); // Say change made on first iteration changeMade_ = 1; // Say last objective infinite //lastObjectiveValue_=-COIN_DBL_MAX; if ((stateOfProblem_ & VALUES_PASS) != 0) { // modify costs so nonbasic dual feasible with fake for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { switch (getInternalStatus(iSequence)) { case basic: case AbcSimplex::isFixed: break; case isFree: case superBasic: //abcCost_[iSequence]-=djSaved_[iSequence]; //abcDj_[iSequence]-=djSaved_[iSequence]; djSaved_[iSequence] = 0.0; break; case atLowerBound: if (djSaved_[iSequence] < 0.0) { //abcCost_[iSequence]-=djSaved_[iSequence]; //abcDj_[iSequence]-=djSaved_[iSequence]; djSaved_[iSequence] = 0.0; } break; case atUpperBound: if (djSaved_[iSequence] > 0.0) { //abcCost_[iSequence]-=djSaved_[iSequence]; //abcDj_[iSequence]-=djSaved_[iSequence]; djSaved_[iSequence] = 0.0; } break; } } } progressFlag_ = 0; /* Status of problem: 0 - optimal 1 - infeasible 2 - unbounded -1 - iterating -2 - factorization wanted -3 - redo checking without factorization -4 - looks infeasible */ int factorizationType = 1; double pivotTolerance = abcFactorization_->pivotTolerance(); while (problemStatus_ < 0) { // If getting nowhere - return double limit = numberRows_ + numberColumns_; limit = 10000.0 + 2000.0 * limit; #if ABC_NORMAL_DEBUG > 0 if (numberIterations_ > limit) { printf("Set flag so ClpSimplexDual done\n"); abort(); } #endif bool disaster = false; if (disasterArea_ && disasterArea_->check()) { disasterArea_->saveInfo(); disaster = true; } if (numberIterations_ > baseIteration_) { // may factorize, checks if problem finished statusOfProblemInDual(factorizationType); factorizationType = 1; largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_); largestDualError = CoinMax(largestDualError, largestDualError_); if (disaster) problemStatus_ = 3; } int saveNumberIterations = numberIterations_; // exit if victory declared if (problemStatus_ >= 0) break; // test for maximum iterations if (hitMaximumIterations()) { problemStatus_ = 3; break; } // Check event int eventStatus = eventHandler_->event(ClpEventHandler::endOfFactorization); if (eventStatus >= 0) { problemStatus_ = 5; secondaryStatus_ = ClpEventHandler::endOfFactorization; break; } if (!numberPrimalInfeasibilities_ && problemStatus_ < 0) { #if ABC_NORMAL_DEBUG > 0 printf("Primal feasible - sum dual %g - go to primal\n", sumDualInfeasibilities_); #endif problemStatus_ = 10; break; } // Do iterations ? loop round #if PARTITION_ROW_COPY == 1 stateOfProblem_ |= NEED_BASIS_SORT; #endif while (true) { #if PARTITION_ROW_COPY == 1 if ((stateOfProblem_ & NEED_BASIS_SORT) != 0) { int iVector = getAvailableArray(); abcMatrix_->sortUseful(usefulArray_[iVector]); setAvailableArray(iVector); stateOfProblem_ &= ~NEED_BASIS_SORT; } #endif problemStatus_ = -1; if ((stateOfProblem_ & VALUES_PASS) != 0) { assert(djSaved_ - abcDj_ == numberTotal_); abcDj_ = djSaved_; djSaved_ = NULL; } #if ABC_PARALLEL if (parallelMode_ == 0) { #endif whileIteratingSerial(); #if ABC_PARALLEL } else { #if ABC_PARALLEL == 1 whileIteratingThread(); #else whileIteratingCilk(); #endif } #endif if ((stateOfProblem_ & VALUES_PASS) != 0) { djSaved_ = abcDj_; abcDj_ = djSaved_ - numberTotal_; } if (pivotRow_ != -100) { if (numberIterations_ > saveNumberIterations || problemStatus_ >= 0) { if (problemStatus_ == 10 && abcFactorization_->pivots() && numberTimesOptimal_ < 3) { factorizationType = 5; problemStatus_ = -4; if (!numberPrimalInfeasibilities_) factorizationType = 9; } else if (problemStatus_ == -2) { factorizationType = 5; } break; } else if (abcFactorization_->pivots() || problemStatus_ == -6 || pivotTolerance < abcFactorization_->pivotTolerance()) { factorizationType = 5; if (pivotTolerance < abcFactorization_->pivotTolerance()) { pivotTolerance = abcFactorization_->pivotTolerance(); if (!abcFactorization_->pivots()) abcFactorization_->minimumPivotTolerance(CoinMin(0.25, 1.05 * abcFactorization_->minimumPivotTolerance())); } break; } else { bool tryFake = numberAtFakeBound() > 0 && currentDualBound_ < 1.0e17; if (problemStatus_ == -5 || tryFake) { if (numberTimesOptimal_ > 10) problemStatus_ = 4; numberTimesOptimal_++; int numberChanged = -1; if (numberFlagged_) { numberFlagged_ = 0; for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; if (flagged(iPivot)) { clearFlagged(iPivot); } } numberChanged = 0; } else if (tryFake) { double dummyChange; numberChanged = changeBounds(0, dummyChange); if (numberChanged > 0) { handler_->message(CLP_DUAL_ORIGINAL, messages_) << CoinMessageEol; bounceTolerances(3 /*(numberChanged>0) ? 3 : 1*/); // temp #if ABC_NORMAL_DEBUG > 0 //printf("should not need to break out of loop\n"); #endif abcProgress_.reset(); break; } else if (numberChanged < -1) { // some flipped bounceTolerances(3); abcProgress_.reset(); } } else if (numberChanged < 0) { // what now printf("No iterations since last inverts A - nothing done - think\n"); abort(); } } else if (problemStatus_ != -4) { // what now printf("No iterations since last inverts B - nothing done - think\n"); abort(); } } } else { // end of values pass // we need to tell statusOf.. to restore costs etc // and take off flag bounceTolerances(2); assert(!solution_); delete[] solution_; solution_ = NULL; gutsOfSolution(3); makeNonFreeVariablesDualFeasible(); stateOfProblem_ &= ~VALUES_PASS; problemStatus_ = -2; factorizationType = 5; break; } } if (problemStatus_ == 1 && (progressFlag_ & 8) != 0 && fabs(rawObjectiveValue_) > 1.0e10) { #if ABC_NORMAL_DEBUG > 0 printf("fix looks infeasible when has been\n"); // goto primal #endif problemStatus_ = 10; // infeasible - but has looked feasible } if (!problemStatus_ && abcFactorization_->pivots()) { gutsOfSolution(1); // need to compute duals //computeInternalObjectiveValue(); checkCutoff(true); } } largestPrimalError_ = largestPrimalError; largestDualError_ = largestDualError; } /// The duals are updated by the given arrays. static void updateDualsInDualBit2(CoinPartitionedVector &row, double *COIN_RESTRICT djs, double theta, int iBlock) { int offset = row.startPartition(iBlock); double *COIN_RESTRICT work = row.denseVector() + offset; const int *COIN_RESTRICT which = row.getIndices() + offset; int number = row.getNumElements(iBlock); for (int i = 0; i < number; i++) { int iSequence = which[i]; double alphaI = work[i]; work[i] = 0.0; double value = djs[iSequence] - theta * alphaI; djs[iSequence] = value; } row.setNumElementsPartition(iBlock, 0); } static void updateDualsInDualBit(CoinPartitionedVector &row, double *djs, double theta, int numberBlocks) { for (int iBlock = 1; iBlock < numberBlocks; iBlock++) { cilk_spawn updateDualsInDualBit2(row, djs, theta, iBlock); } updateDualsInDualBit2(row, djs, theta, 0); cilk_sync; } /// The duals are updated by the given arrays. void AbcSimplexDual::updateDualsInDual() { //dual->usefulArray(3)->compact(); CoinPartitionedVector &array = usefulArray_[arrayForTableauRow_]; //array.compact(); int numberBlocks = array.getNumPartitions(); int number = array.getNumElements(); if (numberBlocks) { if (number > 100) { updateDualsInDualBit(array, abcDj_, theta_, numberBlocks); } else { for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { updateDualsInDualBit2(array, abcDj_, theta_, iBlock); } } array.compact(); } else { double *COIN_RESTRICT work = array.denseVector(); const int *COIN_RESTRICT which = array.getIndices(); double *COIN_RESTRICT reducedCost = abcDj_; #pragma cilk grainsize = 128 cilk_for(int i = 0; i < number; i++) { //for (int i = 0; i < number; i++) { int iSequence = which[i]; double alphaI = work[i]; work[i] = 0.0; double value = reducedCost[iSequence] - theta_ * alphaI; reducedCost[iSequence] = value; } array.setNumElements(0); } } int AbcSimplexDual::nextSuperBasic() { if (firstFree_ >= 0) { // make sure valid int returnValue = firstFree_; int iColumn = firstFree_ + 1; #ifdef PAN if (fakeSuperBasic(firstFree_) < 0) { // somehow cleaned up returnValue = -1; for (; iColumn < numberTotal_; iColumn++) { if (getInternalStatus(iColumn) == isFree || getInternalStatus(iColumn) == superBasic) { if (fakeSuperBasic(iColumn) >= 0) { if (fabs(abcDj_[iColumn]) > 1.0e2 * dualTolerance_) break; } } } if (iColumn < numberTotal_) { returnValue = iColumn; iColumn++; } } #endif for (; iColumn < numberTotal_; iColumn++) { if (getInternalStatus(iColumn) == isFree || getInternalStatus(iColumn) == superBasic) { #ifdef PAN if (fakeSuperBasic(iColumn) >= 0) { #endif if (fabs(abcDj_[iColumn]) > dualTolerance_) break; #ifdef PAN } #endif } } firstFree_ = iColumn; if (firstFree_ == numberTotal_) firstFree_ = -1; return returnValue; } else { return -1; } } /* Chooses dual pivot row For easy problems we can just choose one of the first rows we look at */ void AbcSimplexDual::dualPivotRow() { //int lastPivotRow=pivotRow_; pivotRow_ = -1; const double *COIN_RESTRICT lowerBasic = lowerBasic_; const double *COIN_RESTRICT upperBasic = upperBasic_; const double *COIN_RESTRICT solutionBasic = solutionBasic_; const int *COIN_RESTRICT pivotVariable = abcPivotVariable_; freeSequenceIn_ = -1; double worstValuesPass = 0.0; if ((stateOfProblem_ & VALUES_PASS) != 0) { for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { switch (getInternalStatus(iSequence)) { case basic: #if ABC_NORMAL_DEBUG > 0 if (fabs(abcDj_[iSequence]) > 1.0e-6) printf("valuespass basic dj %d status %d dj %g\n", iSequence, getInternalStatus(iSequence), abcDj_[iSequence]); #endif break; case AbcSimplex::isFixed: break; case isFree: case superBasic: #if ABC_NORMAL_DEBUG > 0 if (fabs(abcDj_[iSequence]) > 1.0e-6) printf("Bad valuespass dj %d status %d dj %g\n", iSequence, getInternalStatus(iSequence), abcDj_[iSequence]); #endif break; case atLowerBound: #if ABC_NORMAL_DEBUG > 0 if (abcDj_[iSequence] < -1.0e-6) printf("Bad valuespass dj %d status %d dj %g\n", iSequence, getInternalStatus(iSequence), abcDj_[iSequence]); #endif break; case atUpperBound: #if ABC_NORMAL_DEBUG > 0 if (abcDj_[iSequence] > 1.0e-6) printf("Bad valuespass dj %d status %d dj %g\n", iSequence, getInternalStatus(iSequence), abcDj_[iSequence]); #endif break; } } // get worst basic dual - later do faster using steepest infeasibility array // does not have to be worst - just bad int iWorst = -1; double worst = dualTolerance_; for (int i = 0; i < numberRows_; i++) { int iSequence = abcPivotVariable_[i]; if (fabs(abcDj_[iSequence]) > worst) { iWorst = i; worst = fabs(abcDj_[iSequence]); } } if (iWorst >= 0) { pivotRow_ = iWorst; worstValuesPass = abcDj_[abcPivotVariable_[iWorst]]; } else { // factorize and clean up costs pivotRow_ = -100; return; } } if (false && numberFreeNonBasic_ > abcFactorization_->maximumPivots()) { lastFirstFree_ = -1; firstFree_ = -1; } if (firstFree_ >= 0) { // free int nextFree = nextSuperBasic(); if (nextFree >= 0) { freeSequenceIn_ = nextFree; // **** should skip price (and maybe weights update) // unpack vector and find a good pivot unpack(usefulArray_[arrayForBtran_], nextFree); abcFactorization_->updateColumn(usefulArray_[arrayForBtran_]); double *COIN_RESTRICT work = usefulArray_[arrayForBtran_].denseVector(); int number = usefulArray_[arrayForBtran_].getNumElements(); int *COIN_RESTRICT which = usefulArray_[arrayForBtran_].getIndices(); double bestFeasibleAlpha = 0.0; int bestFeasibleRow = -1; double bestInfeasibleAlpha = 0.0; int bestInfeasibleRow = -1; // Go for big infeasibilities unless big changes double maxInfeas = (fabs(dualOut_) > 1.0e3) ? 0.001 : 1.0e6; for (int i = 0; i < number; i++) { int iRow = which[i]; double alpha = fabs(work[iRow]); if (alpha > 1.0e-3) { double value = solutionBasic[iRow]; double lower = lowerBasic[iRow]; double upper = upperBasic[iRow]; double infeasibility = 0.0; if (value > upper) infeasibility = value - upper; else if (value < lower) infeasibility = lower - value; infeasibility = CoinMin(maxInfeas, infeasibility); if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) { if (!flagged(pivotVariable[iRow])) { bestInfeasibleAlpha = infeasibility * alpha; bestInfeasibleRow = iRow; } } if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) { bestFeasibleAlpha = alpha; bestFeasibleRow = iRow; } } } if (bestInfeasibleRow >= 0) pivotRow_ = bestInfeasibleRow; else if (bestFeasibleAlpha > 1.0e-2) pivotRow_ = bestFeasibleRow; usefulArray_[arrayForBtran_].clear(); if (pivotRow_ < 0) { // no good freeSequenceIn_ = -1; if (!abcFactorization_->pivots()) { lastFirstFree_ = -1; firstFree_ = -1; } } } } if (pivotRow_ < 0) { // put back so can test //pivotRow_=lastPivotRow; // get pivot row using whichever method it is pivotRow_ = abcDualRowPivot_->pivotRow(); } if (pivotRow_ >= 0) { sequenceOut_ = pivotVariable[pivotRow_]; valueOut_ = solutionBasic[pivotRow_]; lowerOut_ = lowerBasic[pivotRow_]; upperOut_ = upperBasic[pivotRow_]; // if we have problems we could try other way and hope we get a // zero pivot? if (valueOut_ > upperOut_) { directionOut_ = -1; dualOut_ = valueOut_ - upperOut_; } else if (valueOut_ < lowerOut_) { directionOut_ = 1; dualOut_ = lowerOut_ - valueOut_; } else { // odd (could be free) - it's feasible - go to nearest if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) { directionOut_ = 1; dualOut_ = lowerOut_ - valueOut_; } else { directionOut_ = -1; dualOut_ = valueOut_ - upperOut_; } } if (worstValuesPass) { // in values pass so just use sign of dj // We don't want to go through any barriers so set dualOut low // free variables will never be here // need to think about tolerances with inexact solutions // could try crashing in variables away from bounds and with near zero djs // pivot to take out bad slack and increase objective double fakeValueOut = solution_[sequenceOut_]; dualOut_ = 1.0e-10; if (upperOut_ > lowerOut_ || fabs(fakeValueOut - upperOut_) < 100.0 * primalTolerance_) { if (abcDj_[sequenceOut_] > 0.0) { directionOut_ = 1; } else { directionOut_ = -1; } } else { // foolishly trust primal solution if (fakeValueOut > upperOut_) { directionOut_ = -1; } else { directionOut_ = 1; } } #if 1 // make sure plausible but only when fake primals stored if (directionOut_ < 0 && fabs(fakeValueOut - upperOut_) > dualBound_ + primalTolerance_) { if (fabs(fakeValueOut - upperOut_) > fabs(fakeValueOut - lowerOut_)) directionOut_ = 1; } else if (directionOut_ > 0 && fabs(fakeValueOut - lowerOut_) > dualBound_ + primalTolerance_) { if (fabs(fakeValueOut - upperOut_) < fabs(fakeValueOut - lowerOut_)) directionOut_ = -1; } #endif } #if ABC_NORMAL_DEBUG > 3 assert(dualOut_ >= 0.0); #endif } else { sequenceOut_ = -1; } alpha_ = 0.0; upperTheta_ = COIN_DBL_MAX; return; } // Checks if any fake bounds active - if so returns number and modifies // dualBound_ and everything. // Free variables will be left as free // Returns number of bounds changed if >=0 // Returns -1 if not initialize and no effect // Fills in changeVector which can be used to see if unbounded // initialize is 0 or 3 // and cost of change vector int AbcSimplexDual::changeBounds(int initialize, double &changeCost) { assert(initialize == 0 || initialize == 1 || initialize == 3 || initialize == 4); numberFake_ = 0; double *COIN_RESTRICT abcLower = abcLower_; double *COIN_RESTRICT abcUpper = abcUpper_; double *COIN_RESTRICT abcSolution = abcSolution_; double *COIN_RESTRICT abcCost = abcCost_; if (!initialize) { int numberInfeasibilities; double newBound; newBound = 5.0 * currentDualBound_; numberInfeasibilities = 0; changeCost = 0.0; // put back original bounds and then check CoinAbcMemcpy(abcLower, lowerSaved_, numberTotal_); CoinAbcMemcpy(abcUpper, upperSaved_, numberTotal_); const double *COIN_RESTRICT dj = abcDj_; int iSequence; // bounds will get bigger - just look at ones at bounds int numberFlipped = 0; double checkBound = 1.000000000001 * currentDualBound_; for (iSequence = 0; iSequence < numberTotal_; iSequence++) { double lowerValue = abcLower[iSequence]; double upperValue = abcUpper[iSequence]; double value = abcSolution[iSequence]; setFakeBound(iSequence, AbcSimplexDual::noFake); switch (getInternalStatus(iSequence)) { case basic: case AbcSimplex::isFixed: break; case isFree: case superBasic: break; case atUpperBound: if (fabs(value - upperValue) > primalTolerance_) { // can we flip if (value - lowerValue < checkBound && dj[iSequence] > -currentDualTolerance_) { abcSolution_[iSequence] = lowerValue; setInternalStatus(iSequence, atLowerBound); numberFlipped++; } else { numberInfeasibilities++; } } break; case atLowerBound: if (fabs(value - lowerValue) > primalTolerance_) { // can we flip if (upperValue - value < checkBound && dj[iSequence] < currentDualTolerance_) { abcSolution_[iSequence] = upperValue; setInternalStatus(iSequence, atUpperBound); numberFlipped++; } else { numberInfeasibilities++; } } break; } } // If dual infeasible then carry on if (numberInfeasibilities) { handler_->message(CLP_DUAL_CHECKB, messages_) << newBound << CoinMessageEol; int iSequence; for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { double lowerValue = abcLower[iSequence]; double upperValue = abcUpper[iSequence]; double newLowerValue; double newUpperValue; Status status = getInternalStatus(iSequence); if (status == atUpperBound || status == atLowerBound) { double value = abcSolution[iSequence]; if (value - lowerValue <= upperValue - value) { newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound); newUpperValue = CoinMin(upperValue, newLowerValue + newBound); } else { newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound); newLowerValue = CoinMax(lowerValue, newUpperValue - newBound); } abcLower[iSequence] = newLowerValue; abcUpper[iSequence] = newUpperValue; if (newLowerValue > lowerValue) { if (newUpperValue < upperValue) { setFakeBound(iSequence, AbcSimplexDual::bothFake); #ifdef CLP_INVESTIGATE abort(); // No idea what should happen here - I have never got here #endif numberFake_++; } else { setFakeBound(iSequence, AbcSimplexDual::lowerFake); assert(abcLower[iSequence] > lowerSaved_[iSequence]); assert(abcUpper[iSequence] == upperSaved_[iSequence]); numberFake_++; } } else { if (newUpperValue < upperValue) { setFakeBound(iSequence, AbcSimplexDual::upperFake); assert(abcLower[iSequence] == lowerSaved_[iSequence]); assert(abcUpper[iSequence] < upperSaved_[iSequence]); numberFake_++; } } if (status == atUpperBound) abcSolution[iSequence] = newUpperValue; else abcSolution[iSequence] = newLowerValue; double movement = abcSolution[iSequence] - value; changeCost += movement * abcCost[iSequence]; } } currentDualBound_ = newBound; #if ABC_NORMAL_DEBUG > 0 printf("new dual bound %g\n", currentDualBound_); #endif } else { numberInfeasibilities = -1 - numberFlipped; } return numberInfeasibilities; } else { int iSequence; if (initialize == 3) { for (iSequence = 0; iSequence < numberTotal_; iSequence++) { #if ABC_NORMAL_DEBUG > 1 int bad = -1; if (getFakeBound(iSequence) == noFake) { if (abcLower_[iSequence] != lowerSaved_[iSequence] || abcUpper_[iSequence] != upperSaved_[iSequence]) bad = 0; } else if (getFakeBound(iSequence) == lowerFake) { if (abcLower_[iSequence] == lowerSaved_[iSequence] || abcUpper_[iSequence] != upperSaved_[iSequence]) bad = 1; } else if (getFakeBound(iSequence) == upperFake) { if (abcLower_[iSequence] != lowerSaved_[iSequence] || abcUpper_[iSequence] == upperSaved_[iSequence]) bad = 2; } else { if (abcLower_[iSequence] == lowerSaved_[iSequence] || abcUpper_[iSequence] == upperSaved_[iSequence]) bad = 3; } if (bad >= 0) printf("%d status %d l %g,%g u %g,%g\n", iSequence, internalStatus_[iSequence], abcLower_[iSequence], lowerSaved_[iSequence], abcUpper_[iSequence], upperSaved_[iSequence]); #endif setFakeBound(iSequence, AbcSimplexDual::noFake); } CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_); CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_); } double testBound = /*0.999999 **/ currentDualBound_; for (iSequence = 0; iSequence < numberTotal_; iSequence++) { Status status = getInternalStatus(iSequence); if (status == atUpperBound || status == atLowerBound) { double lowerValue = abcLower[iSequence]; double upperValue = abcUpper[iSequence]; double value = abcSolution[iSequence]; if (lowerValue > -largeValue_ || upperValue < largeValue_) { if (true || lowerValue - value > -0.5 * currentDualBound_ || upperValue - value < 0.5 * currentDualBound_) { if (fabs(lowerValue - value) <= fabs(upperValue - value)) { if (upperValue > lowerValue + testBound) { if (getFakeBound(iSequence) == AbcSimplexDual::noFake) numberFake_++; abcUpper[iSequence] = lowerValue + currentDualBound_; double difference = upperSaved_[iSequence] - abcUpper[iSequence]; if (difference < 1.0e-10 * (1.0 + fabs(upperSaved_[iSequence]) + dualBound_)) abcUpper[iSequence] = upperSaved_[iSequence]; assert(abcUpper[iSequence] <= upperSaved_[iSequence] + 1.0e-3 * currentDualBound_); assert(abcLower[iSequence] >= lowerSaved_[iSequence] - 1.0e-3 * currentDualBound_); if (upperSaved_[iSequence] - abcUpper[iSequence] > abcLower[iSequence] - lowerSaved_[iSequence]) { setFakeBound(iSequence, AbcSimplexDual::upperFake); assert(abcLower[iSequence] == lowerSaved_[iSequence]); assert(abcUpper[iSequence] < upperSaved_[iSequence]); } else { setFakeBound(iSequence, AbcSimplexDual::lowerFake); assert(abcLower[iSequence] > lowerSaved_[iSequence]); assert(abcUpper[iSequence] == upperSaved_[iSequence]); } } } else { if (lowerValue < upperValue - testBound) { if (getFakeBound(iSequence) == AbcSimplexDual::noFake) numberFake_++; abcLower[iSequence] = upperValue - currentDualBound_; double difference = abcLower[iSequence] - lowerSaved_[iSequence]; if (difference < 1.0e-10 * (1.0 + fabs(lowerSaved_[iSequence]) + dualBound_)) abcLower[iSequence] = lowerSaved_[iSequence]; assert(abcUpper[iSequence] <= upperSaved_[iSequence] + 1.0e-3 * currentDualBound_); assert(abcLower[iSequence] >= lowerSaved_[iSequence] - 1.0e-3 * currentDualBound_); if (upperSaved_[iSequence] - abcUpper[iSequence] > abcLower[iSequence] - lowerSaved_[iSequence]) { setFakeBound(iSequence, AbcSimplexDual::upperFake); //assert (abcLower[iSequence]==lowerSaved_[iSequence]); abcLower[iSequence] = CoinMax(abcLower[iSequence], lowerSaved_[iSequence]); assert(abcUpper[iSequence] < upperSaved_[iSequence]); } else { setFakeBound(iSequence, AbcSimplexDual::lowerFake); assert(abcLower[iSequence] > lowerSaved_[iSequence]); //assert (abcUpper[iSequence]==upperSaved_[iSequence]); abcUpper[iSequence] = CoinMin(abcUpper[iSequence], upperSaved_[iSequence]); } } } } else { if (getFakeBound(iSequence) == AbcSimplexDual::noFake) numberFake_++; abcLower[iSequence] = -0.5 * currentDualBound_; abcUpper[iSequence] = 0.5 * currentDualBound_; setFakeBound(iSequence, AbcSimplexDual::bothFake); abort(); } if (status == atUpperBound) abcSolution[iSequence] = abcUpper[iSequence]; else abcSolution[iSequence] = abcLower[iSequence]; } else { // set non basic free variables to fake bounds // I don't think we should ever get here CoinAssert(!("should not be here")); abcLower[iSequence] = -0.5 * currentDualBound_; abcUpper[iSequence] = 0.5 * currentDualBound_; setFakeBound(iSequence, AbcSimplexDual::bothFake); numberFake_++; setInternalStatus(iSequence, atUpperBound); abcSolution[iSequence] = 0.5 * currentDualBound_; } } else if (status == basic) { // make sure not at fake bound and bounds correct setFakeBound(iSequence, AbcSimplexDual::noFake); double gap = abcUpper[iSequence] - abcLower[iSequence]; if (gap > 0.5 * currentDualBound_ && gap < 2.0 * currentDualBound_) { abcLower[iSequence] = lowerSaved_[iSequence]; ; abcUpper[iSequence] = upperSaved_[iSequence]; ; } } } return 0; } } /* order of variables is a) atLo/atUp b) free/superBasic pre pass done in AbcMatrix to get possibles dj's will be stored in a parallel array (tempArray) for first pass this pass to be done inside AbcMatrix then as now use perturbationArray as variable tolerance idea is to use large tolerances and then balance and reduce every so often this should be in conjunction with dualBound changes variables always stay at bounds on accuracy issues - only deliberate flips allowed variable tolerance should be 0.8-1.0 desired (at end) use going to clpMatrix_ as last resort */ /* Array has tableau row Puts candidates in list Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free */ void AbcSimplexDual::dualColumn1(bool doAll) { const CoinIndexedVector &update = usefulArray_[arrayForBtran_]; CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_]; CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; if ((stateOfProblem_ & VALUES_PASS) == 0) upperTheta_ = 1.0e31; else upperTheta_ = fabs(abcDj_[sequenceOut_]) - 0.5 * dualTolerance_; assert(upperTheta_ > 0.0); if (!doAll) upperTheta_ = abcMatrix_->dualColumn1(update, tableauRow, candidateList); else upperTheta_ = dualColumn1B(); //tableauRow.compact(); //candidateList.compact(); //tableauRow.setPackedMode(true); //candidateList.setPackedMode(true); } /* Array has tableau row Puts candidates in list Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free */ double AbcSimplexDual::dualColumn1A() { const CoinIndexedVector &update = usefulArray_[arrayForBtran_]; CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_]; CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; int numberRemaining = 0; double upperTheta = upperTheta_; // pivot elements double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); // indices int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); const double *COIN_RESTRICT abcDj = abcDj_; //const double * COIN_RESTRICT abcPerturbation = abcPerturbation_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; const double *COIN_RESTRICT pi = update.denseVector(); const int *COIN_RESTRICT piIndex = update.getIndices(); int number = update.getNumElements(); int *COIN_RESTRICT index = tableauRow.getIndices(); double *COIN_RESTRICT array = tableauRow.denseVector(); assert(!tableauRow.getNumElements()); int numberNonZero = 0; // do first pass to get possibles double bestPossible = 0.0; // We can also see if infeasible or pivoting on free double tentativeTheta = 1.0e25; // try with this much smaller as guess double acceptablePivot = currentAcceptablePivot_; double dualT = -currentDualTolerance_; const double multiplier[] = { 1.0, -1.0 }; if (ordinaryVariables_) { for (int i = 0; i < number; i++) { int iRow = piIndex[i]; unsigned char type = internalStatus[iRow]; if ((type & 4) == 0) { index[numberNonZero] = iRow; double piValue = pi[iRow]; array[numberNonZero++] = piValue; int iStatus = internalStatus[iRow] & 3; double mult = multiplier[iStatus]; double alpha = piValue * mult; double oldValue = abcDj[iRow] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iRow; } } } } else { double badFree = 0.0; double freePivot = currentAcceptablePivot_; alpha_ = 0.0; // We can also see if infeasible or pivoting on free for (int i = 0; i < number; i++) { int iRow = piIndex[i]; unsigned char type = internalStatus[iRow]; double piValue = pi[iRow]; if ((type & 4) == 0) { index[numberNonZero] = iRow; array[numberNonZero++] = piValue; int iStatus = internalStatus[iRow] & 3; double mult = multiplier[iStatus]; double alpha = piValue * mult; double oldValue = abcDj[iRow] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iRow; } } else if ((type & 7) < 6) { //} else if (getInternalStatus(iRow)==isFree||getInternalStatus(iRow)==superBasic) { bool keep; index[numberNonZero] = iRow; array[numberNonZero++] = piValue; bestPossible = CoinMax(bestPossible, fabs(piValue)); double oldValue = abcDj[iRow]; // If free has to be very large - should come in via dualRow //if (getInternalStatus(iRow+addSequence)==isFree&&fabs(piValue)<1.0e-3) //break; if (oldValue > currentDualTolerance_) { keep = true; } else if (oldValue < -currentDualTolerance_) { keep = true; } else { if (fabs(piValue) > CoinMax(10.0 * currentAcceptablePivot_, 1.0e-5)) { keep = true; } else { keep = false; badFree = CoinMax(badFree, fabs(piValue)); } } if (keep) { // free - choose largest if (fabs(piValue) > freePivot) { freePivot = fabs(piValue); sequenceIn_ = iRow; theta_ = oldValue / piValue; alpha_ = piValue; } } } } } //tableauRow.setNumElements(numberNonZero); //candidateList.setNumElements(numberRemaining); tableauRow.setNumElementsPartition(0, numberNonZero); candidateList.setNumElementsPartition(0, numberRemaining); if (!numberRemaining && sequenceIn_ < 0 && upperTheta_ > 1.0e29) { return COIN_DBL_MAX; // Looks infeasible } else { return upperTheta; } } /* Array has tableau row Puts candidates in list Returns upper theta (infinity if no pivot) and may set sequenceIn_ if free */ double AbcSimplexDual::dualColumn1B() { CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_]; CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; tableauRow.compact(); candidateList.clearAndReset(); tableauRow.setPackedMode(true); candidateList.setPackedMode(true); int numberRemaining = 0; double upperTheta = upperTheta_; // pivot elements double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); // indices int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); const double *COIN_RESTRICT abcDj = abcDj_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; int *COIN_RESTRICT index = tableauRow.getIndices(); double *COIN_RESTRICT array = tableauRow.denseVector(); int numberNonZero = tableauRow.getNumElements(); // do first pass to get possibles double bestPossible = 0.0; // We can also see if infeasible or pivoting on free double tentativeTheta = 1.0e25; // try with this much smaller as guess double acceptablePivot = currentAcceptablePivot_; double dualT = -currentDualTolerance_; const double multiplier[] = { 1.0, -1.0 }; if (ordinaryVariables_) { for (int i = 0; i < numberNonZero; i++) { int iSequence = index[i]; unsigned char iStatus = internalStatus[iSequence] & 7; if (iStatus < 4) { double tableauValue = array[i]; if (fabs(tableauValue) > zeroTolerance_) { double mult = multiplier[iStatus]; double alpha = tableauValue * mult; double oldValue = abcDj[iSequence] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iSequence; } } } } } else { double badFree = 0.0; double freePivot = currentAcceptablePivot_; double currentDualTolerance = currentDualTolerance_; for (int i = 0; i < numberNonZero; i++) { int iSequence = index[i]; unsigned char iStatus = internalStatus[iSequence] & 7; if (iStatus < 6) { double tableauValue = array[i]; if (fabs(tableauValue) > zeroTolerance_) { if ((iStatus & 4) == 0) { double mult = multiplier[iStatus]; double alpha = tableauValue * mult; double oldValue = abcDj[iSequence] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iSequence; } } else { bool keep; bestPossible = CoinMax(bestPossible, fabs(tableauValue)); double oldValue = abcDj[iSequence]; // If free has to be very large - should come in via dualRow //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3) //break; if (oldValue > currentDualTolerance) { keep = true; } else if (oldValue < -currentDualTolerance) { keep = true; } else { if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) { keep = true; } else { keep = false; badFree = CoinMax(badFree, fabs(tableauValue)); } } if (keep) { #ifdef PAN if (fakeSuperBasic(iSequence) >= 0) { #endif if (iSequence == freeSequenceIn_) tableauValue = COIN_DBL_MAX; // free - choose largest if (fabs(tableauValue) > fabs(freePivot)) { freePivot = tableauValue; sequenceIn_ = iSequence; theta_ = oldValue / tableauValue; alpha_ = tableauValue; } #ifdef PAN } #endif } } } } } } candidateList.setNumElements(numberRemaining); if (!numberRemaining && sequenceIn_ < 0 && upperTheta_ > 1.0e29) { return COIN_DBL_MAX; // Looks infeasible } else { return upperTheta; } } //#define MORE_CAREFUL //#define PRINT_RATIO_PROGRESS #define MINIMUMTHETA 1.0e-18 #define MAXTRY 100 #define TOL_TYPE -1 #if TOL_TYPE == -1 #define useTolerance1() (0.0) #define useTolerance2() (0.0) #define goToTolerance1() (0.0) #define goToTolerance2() (0.0) #elif TOL_TYPE == 0 #define useTolerance1() (currentDualTolerance_) #define useTolerance2() (currentDualTolerance_) #define goToTolerance1() (0.0) #define goToTolerance2() (0.0) #elif TOL_TYPE == 1 #define useTolerance1() (perturbedTolerance) #define useTolerance2() (0.0) #define goToTolerance1() (0.0) #define goToTolerance2() (0.0) #elif TOL_TYPE == 2 #define useTolerance1() (perturbedTolerance) #define useTolerance2() (perturbedTolerance) #define goToTolerance1() (0.0) #define goToTolerance2() (0.0) #elif TOL_TYPE == 3 #define useTolerance1() (perturbedTolerance) #define useTolerance2() (perturbedTolerance) #define goToTolerance1() (perturbedTolerance) #define goToTolerance2() (0.0) #elif TOL_TYPE == 4 #define useTolerance1() (perturbedTolerance) #define useTolerance2() (perturbedTolerance) #define goToTolerance1() (perturbedTolerance) #define goToTolerance2() (perturbedTolerance) #else XXXX #endif #if CLP_CAN_HAVE_ZERO_OBJ < 2 #define MODIFYCOST 2 #endif #define DONT_MOVE_OBJECTIVE static void dualColumn2Bit(AbcSimplexDual *dual, dualColumnResult *result, int numberBlocks) { for (int iBlock = 1; iBlock < numberBlocks; iBlock++) { cilk_spawn dual->dualColumn2First(result[iBlock]); } dual->dualColumn2First(result[0]); cilk_sync; } void AbcSimplexDual::dualColumn2First(dualColumnResult &result) { CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_]; int iBlock = result.block; int numberSwapped = result.numberSwapped; int numberRemaining = result.numberRemaining; double tentativeTheta = result.tentativeTheta; // pivot elements double *COIN_RESTRICT array = candidateList.denseVector(); // indices int *COIN_RESTRICT indices = candidateList.getIndices(); const double *COIN_RESTRICT abcLower = abcLower_; const double *COIN_RESTRICT abcUpper = abcUpper_; const double *COIN_RESTRICT abcDj = abcDj_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; const double multiplier[] = { 1.0, -1.0 }; // For list of flipped int *COIN_RESTRICT flippedIndex = flipList.getIndices(); double *COIN_RESTRICT flippedElement = flipList.denseVector(); // adjust //if (iBlock>=0) { int offset = candidateList.startPartition(iBlock); assert(offset == flipList.startPartition(iBlock)); array += offset; indices += offset; flippedIndex += offset; flippedElement += offset; //} double thruThis = 0.0; double bestPivot = currentAcceptablePivot_; int bestSequence = -1; int numberInList = numberRemaining; numberRemaining = 0; double upperTheta = 1.0e50; double increaseInThis = 0.0; //objective increase in this loop for (int i = 0; i < numberInList; i++) { int iSequence = indices[i]; assert(getInternalStatus(iSequence) != isFree && getInternalStatus(iSequence) != superBasic); int iStatus = internalStatus[iSequence] & 3; // treat as if at lower bound double mult = multiplier[iStatus]; double alpha = array[i]; double oldValue = abcDj[iSequence] * mult; double value = oldValue - tentativeTheta * alpha; assert(alpha > 0.0); //double perturbedTolerance = abcPerturbation[iSequence]; if (value < -currentDualTolerance_) { // possible here or could go through double range = abcUpper[iSequence] - abcLower[iSequence]; thruThis += range * alpha; increaseInThis += (oldValue - useTolerance1()) * range; // goes on swapped list (also means candidates if too many) flippedElement[numberSwapped] = alpha; flippedIndex[numberSwapped++] = iSequence; if (alpha > bestPivot) { bestPivot = alpha; bestSequence = iSequence; } } else { // won't go through this time value = oldValue - upperTheta * alpha; if (value < -currentDualTolerance_ && alpha >= currentAcceptablePivot_) upperTheta = (oldValue + currentDualTolerance_) / alpha; array[numberRemaining] = alpha; indices[numberRemaining++] = iSequence; } } result.bestEverPivot = bestPivot; result.sequence = bestSequence; result.numberSwapped = numberSwapped; result.numberRemaining = numberRemaining; result.thruThis = thruThis; result.increaseInThis = increaseInThis; result.theta = upperTheta; } void AbcSimplexDual::dualColumn2Most(dualColumnResult &result) { CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_]; int numberBlocksReal = candidateList.getNumPartitions(); flipList.setPartitions(numberBlocksReal, candidateList.startPartitions()); int numberBlocks; if (!numberBlocksReal) numberBlocks = 1; else numberBlocks = numberBlocksReal; //usefulArray_[arrayForTableauRow_].compact(); //candidateList.compact(); int numberSwapped = 0; int numberRemaining = candidateList.getNumElements(); double totalThru = 0.0; // for when variables flip int sequenceIn = -1; double bestEverPivot = currentAcceptablePivot_; // If we think we need to modify costs (not if something from broad sweep) bool modifyCosts = false; // Increase in objective due to swapping bounds (may be negative) double increaseInObjective = 0.0; // pivot elements double *COIN_RESTRICT array = candidateList.denseVector(); // indices int *COIN_RESTRICT indices = candidateList.getIndices(); const double *COIN_RESTRICT abcLower = abcLower_; const double *COIN_RESTRICT abcUpper = abcUpper_; //const double * COIN_RESTRICT abcSolution = abcSolution_; const double *COIN_RESTRICT abcDj = abcDj_; //const double * COIN_RESTRICT abcPerturbation = abcPerturbation_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; // We can also see if infeasible or pivoting on free double tentativeTheta = 1.0e25; // try with this much smaller as guess const double multiplier[] = { 1.0, -1.0 }; // For list of flipped int numberLastSwapped = 0; int *COIN_RESTRICT flippedIndex = flipList.getIndices(); double *COIN_RESTRICT flippedElement = flipList.denseVector(); // If sum of bad small pivots too much #ifdef MORE_CAREFUL bool badSumPivots = false; #endif // not free int lastSequence = -1; double lastTheta = 0.0; double lastPivotValue = 0.0; double thisPivotValue = 0.0; // if free then always choose , otherwise ... double theta = 1.0e50; // now flip flop between spare arrays until reasonable theta tentativeTheta = CoinMax(10.0 * upperTheta_, 1.0e-7); // loops increasing tentative theta until can't go through #ifdef PRINT_RATIO_PROGRESS int iPass = 0; #endif dualColumnResult result2[8]; for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { result2[iBlock].block = iBlock; if (numberBlocksReal) result2[iBlock].numberRemaining = candidateList.getNumElements(iBlock); else result2[iBlock].numberRemaining = candidateList.getNumElements(); result2[iBlock].numberSwapped = 0; result2[iBlock].numberLastSwapped = 0; } double useThru = 0.0; while (tentativeTheta < 1.0e22) { double thruThis = 0.0; double bestPivot = currentAcceptablePivot_; int bestSequence = -1; sequenceIn = -1; thisPivotValue = 0.0; int numberInList = numberRemaining; double upperTheta = 1.0e50; double increaseInThis = 0.0; //objective increase in this loop for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { result2[iBlock].tentativeTheta = tentativeTheta; //result2[iBlock].numberRemaining=numberRemaining; //result2[iBlock].numberSwapped=numberSwapped; //result2[iBlock].numberLastSwapped=numberLastSwapped; } for (int iBlock = 1; iBlock < numberBlocks; iBlock++) { cilk_spawn dualColumn2First(result2[iBlock]); } dualColumn2First(result2[0]); cilk_sync; //dualColumn2Bit(this,result2,numberBlocks); numberSwapped = 0; numberRemaining = 0; for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { if (result2[iBlock].bestEverPivot > bestPivot) { bestPivot = result2[iBlock].bestEverPivot; bestSequence = result2[iBlock].sequence; } numberSwapped += result2[iBlock].numberSwapped; numberRemaining += result2[iBlock].numberRemaining; thruThis += result2[iBlock].thruThis; increaseInThis += result2[iBlock].increaseInThis; if (upperTheta > result2[iBlock].theta) upperTheta = result2[iBlock].theta; } // end of pass #ifdef PRINT_RATIO_PROGRESS iPass++; printf("Iteration %d pass %d dualOut %g thru %g increase %g numberRemain %d\n", numberIterations_, iPass, fabs(dualOut_), totalThru + thruThis, increaseInObjective + increaseInThis, numberRemaining); #endif if (totalThru + thruThis > fabs(dualOut_) - fabs(1.0e-9 * dualOut_) - 1.0e-12 || (increaseInObjective + increaseInThis < 0.0 && bestSequence >= 0) || !numberRemaining) { // We should be pivoting in this batch // so compress down to this lot numberRemaining = 0; candidateList.clearAndReset(); for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { int offset = flipList.startPartition(iBlock); // set so can compact int numberKeep = result2[iBlock].numberLastSwapped; flipList.setNumElementsPartition(iBlock, numberKeep); for (int i = offset + numberKeep; i < offset + result2[iBlock].numberSwapped; i++) { array[numberRemaining] = flippedElement[i]; flippedElement[i] = 0; indices[numberRemaining++] = flippedIndex[i]; } } // make sure candidateList will be cleared correctly candidateList.setNumElements(numberRemaining); candidateList.setPackedMode(true); flipList.compact(); numberSwapped = numberLastSwapped; useThru = totalThru; assert(totalThru <= fabs(dualOut_)); int iTry; // first get ratio with tolerance for (iTry = 0; iTry < MAXTRY; iTry++) { upperTheta = 1.0e50; numberInList = numberRemaining; numberRemaining = 0; increaseInThis = 0.0; //objective increase in this loop thruThis = 0.0; for (int i = 0; i < numberInList; i++) { int iSequence = indices[i]; assert(getInternalStatus(iSequence) != isFree && getInternalStatus(iSequence) != superBasic); int iStatus = internalStatus[iSequence] & 3; // treat as if at lower bound double mult = multiplier[iStatus]; double alpha = array[i]; double oldValue = abcDj[iSequence] * mult; double value = oldValue - upperTheta * alpha; assert(alpha > 0.0); //double perturbedTolerance = abcPerturbation[iSequence]; if (value < -currentDualTolerance_) { if (alpha >= currentAcceptablePivot_) { upperTheta = (oldValue + currentDualTolerance_) / alpha; } } } bestPivot = currentAcceptablePivot_; sequenceIn = -1; double largestPivot = currentAcceptablePivot_; // Sum of bad small pivots #ifdef MORE_CAREFUL double sumBadPivots = 0.0; badSumPivots = false; #endif // Make sure upperTheta will work (-O2 and above gives problems) upperTheta *= 1.0000000001; for (int i = 0; i < numberInList; i++) { int iSequence = indices[i]; int iStatus = internalStatus[iSequence] & 3; // treat as if at lower bound double mult = multiplier[iStatus]; double alpha = array[i]; double oldValue = abcDj[iSequence] * mult; double value = oldValue - upperTheta * alpha; assert(alpha > 0.0); //double perturbedTolerance = abcPerturbation[iSequence]; double badDj = 0.0; if (value <= 0.0) { // add to list of swapped flippedElement[numberSwapped] = alpha; flippedIndex[numberSwapped++] = iSequence; badDj = oldValue - useTolerance2(); // select if largest pivot bool take = (alpha > bestPivot); #ifdef MORE_CAREFUL if (alpha < currentAcceptablePivot_ && upperTheta < 1.0e20) { if (value < -currentDualTolerance_) { double gap = abcUpper[iSequence] - abcLower[iSequence]; if (gap < 1.0e20) sumBadPivots -= value * gap; else sumBadPivots += 1.0e20; //printf("bad %d alpha %g dj at lower %g\n", // iSequence,alpha,value); } } #endif if (take && flagged(iSequence)) { if (bestPivot > currentAcceptablePivot_) take = false; } if (take) { sequenceIn = iSequence; thisPivotValue = alpha; bestPivot = alpha; if (flagged(iSequence)) bestPivot = currentAcceptablePivot_; bestSequence = iSequence; theta = (oldValue + goToTolerance1()) / alpha; //assert (oldValue==abcDj[iSequence]); largestPivot = CoinMax(largestPivot, 0.5 * bestPivot); } double range = abcUpper[iSequence] - abcLower[iSequence]; thruThis += range * fabs(alpha); increaseInThis += badDj * range; } else { // add to list of remaining array[numberRemaining] = alpha; indices[numberRemaining++] = iSequence; } } #ifdef MORE_CAREFUL // If we have done pivots and things look bad set alpha_ 0.0 to force factorization if (sumBadPivots > 1.0e4) { #if ABC_NORMAL_DEBUG > 0 if (handler_->logLevel() > 1) printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots, abcFactorization_->pivots()); #endif if (abcFactorization_->pivots() > 3) { badSumPivots = true; #ifdef PRINT_RATIO_PROGRESS printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots, abcFactorization_->pivots()); #endif break; } } #endif // If we stop now this will be increase in objective (I think) double increase = (fabs(dualOut_) - totalThru) * theta; increase += increaseInObjective; if (theta < 0.0) thruThis += fabs(dualOut_); // force using this one if (((increaseInObjective < 0.0 && increase < 0.0) || (theta > 1.0e9 && lastTheta < 1.0e-5)) && lastSequence >= 0) { // back // We may need to be more careful - we could do by // switch so we always do fine grained? bestPivot = 0.0; bestSequence = -1; } else { // add in totalThru += thruThis; increaseInObjective += increaseInThis; lastTheta = theta; } if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) { // back to previous one sequenceIn = lastSequence; thisPivotValue = lastPivotValue; #ifdef PRINT_RATIO_PROGRESS printf("Try %d back to previous bestPivot %g bestEver %g totalThru %g\n", iTry, bestPivot, bestEverPivot, totalThru); #endif break; } else if (sequenceIn == -1 && upperTheta > largeValue_) { if (lastPivotValue > currentAcceptablePivot_) { // back to previous one sequenceIn = lastSequence; thisPivotValue = lastPivotValue; #ifdef PRINT_RATIO_PROGRESS printf("Try %d no pivot this time large theta %g lastPivot %g\n", iTry, upperTheta, lastPivotValue); #endif } else { // can only get here if all pivots too small #ifdef PRINT_RATIO_PROGRESS printf("XXBAD Try %d no pivot this time large theta %g lastPivot %g\n", iTry, upperTheta, lastPivotValue); #endif } break; } else if (totalThru >= fabs(dualOut_)) { #ifdef PRINT_RATIO_PROGRESS printf("Try %d upperTheta %g totalThru %g - can modify costs\n", iTry, upperTheta, lastPivotValue); #endif modifyCosts = true; // fine grain - we can modify costs break; // no point trying another loop } else { lastSequence = sequenceIn; lastPivotValue = thisPivotValue; if (bestPivot > bestEverPivot) { bestEverPivot = bestPivot; //bestEverSequence=bestSequence; } modifyCosts = true; // fine grain - we can modify costs } #ifdef PRINT_RATIO_PROGRESS printf("Onto next pass totalthru %g bestPivot %g\n", totalThru, lastPivotValue); #endif numberLastSwapped = numberSwapped; } break; } else { // skip this lot if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) { #ifdef PRINT_RATIO_PROGRESS printf("Continuing bestPivot %g bestEver %g %d swapped - new tentative %g\n", bestPivot, bestEverPivot, numberSwapped, 2 * upperTheta); #endif bestEverPivot = bestPivot; //bestEverSequence=bestSequence; lastPivotValue = bestPivot; lastSequence = bestSequence; for (int iBlock = 0; iBlock < numberBlocks; iBlock++) result2[iBlock].numberLastSwapped = result2[iBlock].numberSwapped; numberLastSwapped = numberSwapped; } else { #ifdef PRINT_RATIO_PROGRESS printf("keep old swapped bestPivot %g bestEver %g %d swapped - new tentative %g\n", bestPivot, bestEverPivot, numberSwapped, 2 * upperTheta); #endif // keep old swapped numberRemaining = 0; for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { int offset = flipList.startPartition(iBlock); int numberRemaining2 = offset + result2[iBlock].numberRemaining; int numberLastSwapped = result2[iBlock].numberLastSwapped; for (int i = offset + numberLastSwapped; i < offset + result2[iBlock].numberSwapped; i++) { array[numberRemaining2] = flippedElement[i]; flippedElement[i] = 0; indices[numberRemaining2++] = flippedIndex[i]; } result2[iBlock].numberRemaining = numberRemaining2 - offset; numberRemaining += result2[iBlock].numberRemaining; result2[iBlock].numberSwapped = numberLastSwapped; } numberSwapped = numberLastSwapped; } increaseInObjective += increaseInThis; tentativeTheta = 2.0 * upperTheta; totalThru += thruThis; } } if (flipList.getNumPartitions()) { // need to compress numberRemaining = 0; candidateList.clearAndReset(); for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { int offset = flipList.startPartition(iBlock); // set so can compact int numberKeep = result2[iBlock].numberLastSwapped; flipList.setNumElementsPartition(iBlock, numberKeep); for (int i = offset + numberKeep; i < offset + result2[iBlock].numberSwapped; i++) { array[numberRemaining] = flippedElement[i]; flippedElement[i] = 0; indices[numberRemaining++] = flippedIndex[i]; } } // make sure candidateList will be cleared correctly candidateList.setNumElements(numberRemaining); candidateList.setPackedMode(true); flipList.compact(); numberSwapped = numberLastSwapped; } result.theta = theta; //? result.totalThru = totalThru; result.useThru = useThru; result.bestEverPivot = bestEverPivot; //? result.increaseInObjective = increaseInObjective; //? result.tentativeTheta = tentativeTheta; //? result.lastPivotValue = lastPivotValue; result.thisPivotValue = thisPivotValue; result.lastSequence = lastSequence; result.sequence = sequenceIn; //result.block; //result.pass; //result.phase; result.numberSwapped = numberSwapped; result.numberLastSwapped = numberLastSwapped; result.modifyCosts = modifyCosts; } /* Chooses incoming Puts flipped ones in list If necessary will modify costs */ void AbcSimplexDual::dualColumn2() { CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; CoinPartitionedVector &flipList = usefulArray_[arrayForFlipBounds_]; //usefulArray_[arrayForTableauRow_].compact(); usefulArray_[arrayForTableauRow_].computeNumberElements(); if (candidateList.getNumElements() < 100) candidateList.compact(); //usefulArray_[arrayForTableauRow_].compact(); // But factorizations complain if <1.0e-8 //acceptablePivot=CoinMax(acceptablePivot,1.0e-8); // do ratio test for normal iteration currentAcceptablePivot_ = 1.0e-1 * acceptablePivot_; if (numberIterations_ > 100) currentAcceptablePivot_ = acceptablePivot_; if (abcFactorization_->pivots() > 10 || (abcFactorization_->pivots() && sumDualInfeasibilities_)) currentAcceptablePivot_ = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict else if (abcFactorization_->pivots() > 5) currentAcceptablePivot_ = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict else if (abcFactorization_->pivots()) currentAcceptablePivot_ = acceptablePivot_; // relax // If very large infeasibility use larger acceptable pivot if (abcFactorization_->pivots()) { double value = CoinMax(currentAcceptablePivot_, 1.0e-12 * fabs(dualOut_)); currentAcceptablePivot_ = CoinMin(1.0e2 * currentAcceptablePivot_, value); } // return at once if no free but there were free if (lastFirstFree_ >= 0 && sequenceIn_ < 0) { candidateList.clearAndReset(); upperTheta_ = COIN_DBL_MAX; } if (upperTheta_ == COIN_DBL_MAX) { usefulArray_[arrayForTableauRow_].compact(); assert(!candidateList.getNumElements()); return; } int numberSwapped = 0; //int numberRemaining = candidateList.getNumElements(); double useThru = 0.0; // for when variables flip // If we think we need to modify costs (not if something from broad sweep) bool modifyCosts = false; // pivot elements //double * COIN_RESTRICT array=candidateList.denseVector(); // indices const double *COIN_RESTRICT abcLower = abcLower_; const double *COIN_RESTRICT abcUpper = abcUpper_; const double *COIN_RESTRICT abcSolution = abcSolution_; const double *COIN_RESTRICT abcDj = abcDj_; //const double * COIN_RESTRICT abcPerturbation = abcPerturbation_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; // We can also see if infeasible or pivoting on free const double multiplier[] = { 1.0, -1.0 }; // For list of flipped int numberLastSwapped = 0; int *COIN_RESTRICT flippedIndex = flipList.getIndices(); double *COIN_RESTRICT flippedElement = flipList.denseVector(); // If sum of bad small pivots too much #ifdef MORE_CAREFUL bool badSumPivots = false; #endif if (sequenceIn_ < 0) { int lastSequence = -1; double lastPivotValue = 0.0; double thisPivotValue = 0.0; dualColumnResult result; dualColumn2Most(result); sequenceIn_ = result.sequence; lastSequence = result.lastSequence; thisPivotValue = result.thisPivotValue; lastPivotValue = result.lastPivotValue; numberSwapped = result.numberSwapped; numberLastSwapped = result.numberLastSwapped; modifyCosts = result.modifyCosts; useThru = result.useThru; #ifdef PRINT_RATIO_PROGRESS printf("Iteration %d - out of loop - sequenceIn %d pivotValue %g\n", numberIterations_, sequenceIn_, thisPivotValue); #endif // can get here without sequenceIn_ set but with lastSequence if (sequenceIn_ < 0 && lastSequence >= 0) { // back to previous one sequenceIn_ = lastSequence; thisPivotValue = lastPivotValue; #ifdef PRINT_RATIO_PROGRESS printf("BUT back one - out of loop - sequenceIn %d pivotValue %g\n", sequenceIn_, thisPivotValue); #endif } // Movement should be minimum for anti-degeneracy - unless // fixed variable out // no need if just one candidate assert(numberSwapped - numberLastSwapped >= 0); double minimumTheta; if (upperOut_ > lowerOut_ && numberSwapped - numberLastSwapped > 1) #ifndef TRY_ABC_GUS minimumTheta = MINIMUMTHETA; #else minimumTheta = minimumThetaMovement_; #endif else minimumTheta = 0.0; if (sequenceIn_ >= 0) { double oldValue; alpha_ = thisPivotValue; // correct sign int iStatus = internalStatus[sequenceIn_] & 3; alpha_ *= multiplier[iStatus]; oldValue = abcDj[sequenceIn_]; theta_ = CoinMax(oldValue / alpha_, 0.0); if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) { // can't pivot to zero theta_ = minimumTheta; } // may need to adjust costs so all dual feasible AND pivoted is exactly 0 if (modifyCosts #ifdef MORE_CAREFUL && !badSumPivots #endif ) { for (int i = numberLastSwapped; i < numberSwapped; i++) { int iSequence = flippedIndex[i]; int iStatus = internalStatus[iSequence] & 3; // treat as if at lower bound double mult = multiplier[iStatus]; double alpha = flippedElement[i]; flippedElement[i] = 0.0; if (iSequence == sequenceIn_) continue; double oldValue = abcDj[iSequence] * mult; double value = oldValue - theta_ * alpha; //double perturbedTolerance = abcPerturbation[iSequence]; if (-value > currentDualTolerance_) { //thisIncrease = true; #if MODIFYCOST if (numberIterations_ >= normalDualColumnIteration_) { // modify cost to hit new tolerance double modification = alpha * theta_ - oldValue - currentDualTolerance_; if ((specialOptions_ & (2048 + 4096)) != 0) { if ((specialOptions_ & 2048) != 0) { if (fabs(modification) < 1.0e-10) modification = 0.0; } else { if (fabs(modification) < 1.0e-12) modification = 0.0; } } abcDj_[iSequence] += mult * modification; abcCost_[iSequence] += mult * modification; if (modification) numberChanged_++; // Say changed costs } #endif } } numberSwapped = numberLastSwapped; } } } #ifdef MORE_CAREFUL // If we have done pivots and things look bad set alpha_ 0.0 to force factorization if ((badSumPivots || fabs(theta_ * badFree) > 10.0 * currentDualTolerance_) && abcFactorization_->pivots()) { #if ABC_NORMAL_DEBUG > 0 if (handler_->logLevel() > 1) printf("forcing re-factorization\n"); #endif sequenceIn_ = -1; } #endif if (sequenceIn_ >= 0) { lowerIn_ = abcLower[sequenceIn_]; upperIn_ = abcUpper[sequenceIn_]; valueIn_ = abcSolution[sequenceIn_]; dualIn_ = abcDj_[sequenceIn_]; if (numberTimesOptimal_) { // can we adjust cost back closer to original //*** add coding } if (numberIterations_ >= normalDualColumnIteration_) { #if MODIFYCOST > 1 // modify cost to hit zero exactly // so (dualIn_+modification)==theta_*alpha_ //double perturbedTolerance = abcPerturbation[sequenceIn_]; double modification = theta_ * alpha_ - (dualIn_ + goToTolerance2()); // But should not move objective too much ?? #ifdef DONT_MOVE_OBJECTIVE double moveObjective = fabs(modification * abcSolution[sequenceIn_]); double smallMove = CoinMax(fabs(rawObjectiveValue_), 1.0e-3); if (moveObjective > smallMove) { #if ABC_NORMAL_DEBUG > 0 if (handler_->logLevel() > 1) printf("would move objective by %g - original mod %g sol value %g\n", moveObjective, modification, abcSolution[sequenceIn_]); #endif modification *= smallMove / moveObjective; } #endif #ifdef MORE_CAREFUL if (badSumPivots) modification = 0.0; #endif if (atFakeBound(sequenceIn_) && fabs(modification) < currentDualTolerance_) modification = 0.0; if ((specialOptions_ & (2048 + 4096)) != 0) { if ((specialOptions_ & 16384) != 0) { // in fast dual if (fabs(modification) < 1.0e-7) modification = 0.0; } else if ((specialOptions_ & 2048) != 0) { if (fabs(modification) < 1.0e-10) modification = 0.0; } else { if (fabs(modification) < 1.0e-12) modification = 0.0; } } dualIn_ += modification; abcDj_[sequenceIn_] = dualIn_; abcCost_[sequenceIn_] += modification; if (modification) numberChanged_++; // Say changed costs //int costOffset = numberRows_+numberColumns_; //abcCost[sequenceIn_+costOffset] += modification; // save change //assert (fabs(modification)<1.0e-6); #if ABC_NORMAL_DEBUG > 3 if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15) printf("exact %d new cost %g, change %g\n", sequenceIn_, abcCost_[sequenceIn_], modification); #endif #endif } if (alpha_ < 0.0) { // as if from upper bound directionIn_ = -1; upperIn_ = valueIn_; } else { // as if from lower bound directionIn_ = 1; lowerIn_ = valueIn_; } } else { // no pivot alpha_ = 0.0; } // clear array //printf("dualout %g theta %g alpha %g dj %g thru %g swapped %d product %g guess %g useThru %g\n", // dualOut_,theta_,alpha_,abcDj_[sequenceIn_],useThru,numberSwapped, // fabs(dualOut_)*theta_,objectiveValue()+(fabs(dualOut_)-useThru)*theta_,useThru); objectiveChange_ = (fabs(dualOut_) - useThru) * theta_; rawObjectiveValue_ += objectiveChange_; //CoinZeroN(array, maximumNumberInList); //candidateList.setNumElements(0); candidateList.clearAndReset(); flipList.setNumElements(numberSwapped); if (numberSwapped) flipList.setPackedMode(true); #ifdef PRINT_RATIO_PROGRESS printf("%d flipped\n", numberSwapped); #endif } /* Checks if tentative optimal actually means unbounded Returns -3 if not, 2 if is unbounded */ int AbcSimplexDual::checkUnbounded(CoinIndexedVector &ray, double changeCost) { int status = 2; // say unbounded abcFactorization_->updateColumn(ray); // get reduced cost int number = ray.getNumElements(); int *COIN_RESTRICT index = ray.getIndices(); double *COIN_RESTRICT array = ray.denseVector(); const int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_; for (int i = 0; i < number; i++) { int iRow = index[i]; int iPivot = abcPivotVariable[iRow]; changeCost -= cost(iPivot) * array[iRow]; } double way; if (changeCost > 0.0) { //try going down way = 1.0; } else if (changeCost < 0.0) { //try going up way = -1.0; } else { #if ABC_NORMAL_DEBUG > 3 printf("can't decide on up or down\n"); #endif way = 0.0; status = -3; } double movement = 1.0e10 * way; // some largish number double zeroTolerance = 1.0e-14 * dualBound_; for (int i = 0; i < number; i++) { int iRow = index[i]; int iPivot = abcPivotVariable[iRow]; double arrayValue = array[iRow]; if (fabs(arrayValue) < zeroTolerance) arrayValue = 0.0; double newValue = solution(iPivot) + movement * arrayValue; if (newValue > upper(iPivot) + primalTolerance_ || newValue < lower(iPivot) - primalTolerance_) status = -3; // not unbounded } if (status == 2) { // create ray delete[] ray_; ray_ = new double[numberColumns_]; CoinZeroN(ray_, numberColumns_); for (int i = 0; i < number; i++) { int iRow = index[i]; int iPivot = abcPivotVariable[iRow]; double arrayValue = array[iRow]; if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance) ray_[iPivot] = way * array[iRow]; } } ray.clear(); return status; } // Saves good status etc void AbcSimplex::saveGoodStatus() { // save arrays CoinMemcpyN(internalStatus_, numberTotal_, internalStatusSaved_); // save costs CoinMemcpyN(abcCost_, numberTotal_, abcCost_ + 2 * numberTotal_); CoinMemcpyN(abcSolution_, numberTotal_, solutionSaved_); } // Restores previous good status and says trouble void AbcSimplex::restoreGoodStatus(int type) { #if PARTITION_ROW_COPY == 1 stateOfProblem_ |= NEED_BASIS_SORT; #endif if (type) { if (alphaAccuracy_ != -1.0) alphaAccuracy_ = -2.0; // trouble - restore solution CoinMemcpyN(internalStatusSaved_, numberTotal_, internalStatus_); CoinMemcpyN(abcCost_ + 2 * numberTotal_, numberTotal_, abcCost_); CoinMemcpyN(solutionSaved_, numberTotal_, abcSolution_); int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_; // redo pivotVariable int numberBasic = 0; for (int i = 0; i < numberTotal_; i++) { if (getInternalStatus(i) == basic) abcPivotVariable[numberBasic++] = i; } assert(numberBasic == numberRows_); forceFactorization_ = 1; // a bit drastic but .. changeMade_++; // say something changed } } static int computePrimalsAndCheck(AbcSimplexDual *model, const int whichArray[2]) { CoinIndexedVector *array1 = model->usefulArray(whichArray[0]); CoinIndexedVector *array2 = model->usefulArray(whichArray[1]); int numberRefinements = model->computePrimals(array1, array2); model->checkPrimalSolution(true); return numberRefinements; } static int computeDualsAndCheck(AbcSimplexDual *model, const int whichArray[2]) { CoinIndexedVector *array1 = model->usefulArray(whichArray[0]); CoinIndexedVector *array2 = model->usefulArray(whichArray[1]); int numberRefinements = model->computeDuals(NULL, array1, array2); model->checkDualSolutionPlusFake(); return numberRefinements; } // Computes solutions - 1 do duals, 2 do primals, 3 both int AbcSimplex::gutsOfSolution(int type) { double oldLargestPrimalError = largestPrimalError_; double oldLargestDualError = largestDualError_; AbcSimplexDual *dual = reinterpret_cast< AbcSimplexDual * >(this); // do work and check int numberRefinements = 0; #if ABC_PARALLEL if (type != 3) { #endif if ((type & 2) != 0) { //work space int whichArray[2]; for (int i = 0; i < 2; i++) whichArray[i] = getAvailableArray(); numberRefinements = computePrimalsAndCheck(dual, whichArray); for (int i = 0; i < 2; i++) setAvailableArray(whichArray[i]); } if ((type & 1) != 0 #if CAN_HAVE_ZERO_OBJ > 1 && (specialOptions_ & 2097152) == 0 #endif ) { //work space int whichArray[2]; for (int i = 0; i < 2; i++) whichArray[i] = getAvailableArray(); numberRefinements += computeDualsAndCheck(dual, whichArray); for (int i = 0; i < 2; i++) setAvailableArray(whichArray[i]); } #if ABC_PARALLEL } else { #ifndef NDEBUG abcFactorization_->checkMarkArrays(); #endif // redo to do in parallel //work space int whichArray[5]; for (int i = 1; i < 5; i++) whichArray[i] = getAvailableArray(); if (parallelMode_ == 0) { numberRefinements += computePrimalsAndCheck(dual, whichArray + 3); numberRefinements += computeDualsAndCheck(dual, whichArray + 1); } else { #if ABC_PARALLEL == 1 threadInfo_[0].stuff[0] = whichArray[1]; threadInfo_[0].stuff[1] = whichArray[2]; stopStart_ = 1 + 32 * 1; startParallelStuff(1); #else whichArray[0] = 1; CoinThreadInfo info; info.status = 1; info.stuff[0] = whichArray[1]; info.stuff[1] = whichArray[2]; int n = cilk_spawn computeDualsAndCheck(dual, whichArray + 1); #endif numberRefinements = computePrimalsAndCheck(dual, whichArray + 3); #if ABC_PARALLEL == 1 numberRefinements += stopParallelStuff(1); #else cilk_sync; numberRefinements += n; #endif } for (int i = 1; i < 5; i++) setAvailableArray(whichArray[i]); } #endif rawObjectiveValue_ += sumNonBasicCosts_; setClpSimplexObjectiveValue(); if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2)) handler_->message(CLP_SIMPLEX_ACCURACY, messages_) << largestPrimalError_ << largestDualError_ << CoinMessageEol; if ((largestPrimalError_ > 1.0e1 || largestDualError_ > 1.0e1) && numberRows_ > 100 && numberIterations_) { #if ABC_NORMAL_DEBUG > 0 printf("Large errors - primal %g dual %g\n", largestPrimalError_, largestDualError_); #endif // Change factorization tolerance //if (abcFactorization_->zeroTolerance() > 1.0e-18) //abcFactorization_->zeroTolerance(1.0e-18); } bool notChanged = true; if (numberIterations_ && (forceFactorization_ > 2 || forceFactorization_ < 0 || abcFactorization_->pivotTolerance() < 0.9899999999) && (oldLargestDualError || oldLargestPrimalError)) { double useOldDualError = oldLargestDualError; double useDualError = largestDualError_; if (algorithm_ > 0 && abcNonLinearCost_ && abcNonLinearCost_->sumInfeasibilities()) { double factor = CoinMax(1.0, CoinMin(1.0e3, infeasibilityCost_ * 1.0e-6)); useOldDualError /= factor; useDualError /= factor; } if ((largestPrimalError_ > 1.0e3 && oldLargestPrimalError * 1.0e2 < largestPrimalError_) || (useDualError > 1.0e3 && useOldDualError * 1.0e2 < useDualError)) { double pivotTolerance = abcFactorization_->pivotTolerance(); double factor = (largestPrimalError_ > 1.0e10 || largestDualError_ > 1.0e10) ? 2.0 : 1.2; if (pivotTolerance < 0.1) abcFactorization_->pivotTolerance(0.1); else if (pivotTolerance < 0.98999999) abcFactorization_->pivotTolerance(CoinMin(0.99, pivotTolerance * factor)); notChanged = pivotTolerance == abcFactorization_->pivotTolerance(); #if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9) if (pivotTolerance < 0.9899999) { double newTolerance = abcFactorization_->pivotTolerance(); printf("Changing pivot tolerance from %g to %g and backtracking\n", pivotTolerance, newTolerance); } printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n", oldLargestPrimalError, largestPrimalError_, oldLargestDualError, largestDualError_, pivotTolerance); #endif if (pivotTolerance < 0.9899999) { //largestPrimalError_=0.0; //largestDualError_=0.0; //returnCode=1; } } } if (progress_.iterationNumber_[0] > 0 && progress_.iterationNumber_[CLP_PROGRESS - 1] - progress_.iterationNumber_[0] < CLP_PROGRESS * 3 && abcFactorization_->pivotTolerance() < 0.25 && notChanged) { double pivotTolerance = abcFactorization_->pivotTolerance(); abcFactorization_->pivotTolerance(pivotTolerance * 1.5); #if defined(CLP_USEFUL_PRINTOUT) && !defined(GCC_4_9) double newTolerance = abcFactorization_->pivotTolerance(); printf("Changing pivot tolerance from %g to %g - inverting too often\n", pivotTolerance, newTolerance); #endif } return numberRefinements; } // Make non free variables dual feasible by moving to a bound int AbcSimplexDual::makeNonFreeVariablesDualFeasible(bool changeCosts) { const double multiplier[] = { 1.0, -1.0 }; // sign????????????? double dualT = -100.0 * currentDualTolerance_; int numberMoved = 0; int numberChanged = 0; double *COIN_RESTRICT solution = abcSolution_; double *COIN_RESTRICT lower = abcLower_; double *COIN_RESTRICT upper = abcUpper_; double saveNonBasicCosts = sumNonBasicCosts_; double largestCost = dualTolerance_; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { unsigned char iStatus = internalStatus_[iSequence] & 7; if (iStatus < 4) { double mult = multiplier[iStatus]; double dj = abcDj_[iSequence] * mult; if (dj < dualT) { largestCost = CoinMax(fabs(abcPerturbation_[iSequence]), largestCost); } } else if (iStatus != 7) { largestCost = CoinMax(fabs(abcPerturbation_[iSequence]), largestCost); } } for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { unsigned char iStatus = internalStatus_[iSequence] & 7; if (iStatus < 4) { double mult = multiplier[iStatus]; double dj = abcDj_[iSequence] * mult; // ? should use perturbationArray if (dj < dualT && changeCosts) { double gap = upper[iSequence] - lower[iSequence]; if (gap > 1.1 * currentDualBound_) { assert(getFakeBound(iSequence) == AbcSimplexDual::noFake); if (!iStatus) { upper[iSequence] = lower[iSequence] + currentDualBound_; setFakeBound(iSequence, AbcSimplexDual::upperFake); } else { lower[iSequence] = upper[iSequence] - currentDualBound_; setFakeBound(iSequence, AbcSimplexDual::lowerFake); } } double costDifference = abcPerturbation_[iSequence] - abcCost_[iSequence]; if (costDifference * mult > 0.0) { #if ABC_NORMAL_DEBUG > 0 //printf("can improve situation on %d by %g\n",iSequence,costDifference*mult); #endif dj += costDifference * mult; abcDj_[iSequence] += costDifference; abcCost_[iSequence] = abcPerturbation_[iSequence]; costDifference = 0.0; } double changeInObjective = -gap * dj; // probably need to be more intelligent if ((changeInObjective > 1.0e3 || gap > CoinMax(0.5 * dualBound_, 1.0e3)) && dj > -1.0e-1) { // change cost costDifference = (-dj - 0.5 * dualTolerance_) * mult; if (fabs(costDifference) < 100000.0 * dualTolerance_) { dj += costDifference * mult; abcDj_[iSequence] += costDifference; abcCost_[iSequence] += costDifference; numberChanged++; } } } if (dj < dualT && !changeCosts) { numberMoved++; if (iStatus == 0) { // at lower bound // to upper bound setInternalStatus(iSequence, atUpperBound); solution[iSequence] = upper[iSequence]; sumNonBasicCosts_ += abcCost_[iSequence] * (upper[iSequence] - lower[iSequence]); } else { // at upper bound // to lower bound setInternalStatus(iSequence, atLowerBound); solution[iSequence] = lower[iSequence]; sumNonBasicCosts_ -= abcCost_[iSequence] * (upper[iSequence] - lower[iSequence]); } } } } #if ABC_NORMAL_DEBUG > 0 if (numberMoved + numberChanged) printf("makeDualfeasible %d moved, %d had costs changed - sum dual %g\n", numberMoved, numberChanged, sumDualInfeasibilities_); #endif rawObjectiveValue_ += (sumNonBasicCosts_ - saveNonBasicCosts); setClpSimplexObjectiveValue(); return numberMoved; } int AbcSimplexDual::whatNext() { int returnCode = 0; assert(!stateOfIteration_); // see if update stable #if ABC_NORMAL_DEBUG > 3 if ((handler_->logLevel() & 32)) printf("btran alpha %g, ftran alpha %g\n", btranAlpha_, alpha_); #endif int numberPivots = abcFactorization_->pivots(); //double checkValue = 1.0e-7; // numberPivots<5 ? 1.0e-7 : 1.0e-6; double checkValue = numberPivots ? 1.0e-7 : 1.0e-5; // relax if free variable in if ((internalStatus_[sequenceIn_] & 7) > 1) checkValue = 1.0e-5; // if can't trust much and long way from optimal then relax if (largestPrimalError_ > 10.0) checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_); if (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > checkValue * (1.0 + fabs(alpha_))) { handler_->message(CLP_DUAL_CHECK, messages_) << btranAlpha_ << alpha_ << CoinMessageEol; // see with more relaxed criterion double test; if (fabs(btranAlpha_) < 1.0e-8 || fabs(alpha_) < 1.0e-8) test = 1.0e-1 * fabs(alpha_); else test = 10.0 * checkValue; //1.0e-4 * (1.0 + fabs(alpha_)); bool needFlag = (fabs(btranAlpha_) < 1.0e-12 || fabs(alpha_) < 1.0e-12 || fabs(btranAlpha_ - alpha_) > test); #if ABC_NORMAL_DEBUG > 0 if (fabs(btranAlpha_ + alpha_) < 1.0e-5 * (1.0 + fabs(alpha_))) { printf("alpha diff (%g,%g) %g check %g\n", btranAlpha_, alpha_, fabs(btranAlpha_ - alpha_), checkValue * (1.0 + fabs(alpha_))); } #endif if (numberPivots) { if (needFlag && numberPivots < 10) { // need to reject something assert(sequenceOut_ >= 0); char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 0 printf("flag %d as alpha %g %g - after %d pivots\n", sequenceOut_, btranAlpha_, alpha_, numberPivots); #endif // Make safer? abcFactorization_->saferTolerances(1.0, -1.03); setFlagged(sequenceOut_); // so can't happen immediately again sequenceOut_ = -1; lastBadIteration_ = numberIterations_; // say be more cautious } // Make safer? //if (numberPivots<5) //abcFactorization_->saferTolerances (-0.99, -1.01); problemStatus_ = -2; // factorize now returnCode = -2; stateOfIteration_ = 2; } else { if (needFlag) { // need to reject something assert(sequenceOut_ >= 0); char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 3 printf("flag a %g %g\n", btranAlpha_, alpha_); #endif // Make safer? double tolerance = abcFactorization_->pivotTolerance(); abcFactorization_->saferTolerances(1.0, -1.03); setFlagged(sequenceOut_); // so can't happen immediately again sequenceOut_ = -1; //abcProgress_.clearBadTimes(); lastBadIteration_ = numberIterations_; // say be more cautious if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha_) < 1.0e-8 && numberIterations_ > 100) { //printf("I think should declare infeasible\n"); problemStatus_ = 1; returnCode = 1; stateOfIteration_ = 2; } else { stateOfIteration_ = 1; if (tolerance < abcFactorization_->pivotTolerance()) { stateOfIteration_ = 2; returnCode = -2; } } } } } if (!stateOfIteration_) { // check update returnCode = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_); } return returnCode; } /* Checks if finished. type 0 initial 1 normal 2 initial but we have solved problem before 3 can skip factorization Updates status */ void AbcSimplexDual::statusOfProblemInDual(int type) { // If lots of iterations then adjust costs if large ones int numberPivots = abcFactorization_->pivots(); #if 0 printf("statusOf %d pivots type %d pstatus %d forceFac %d dont %d\n", numberPivots,type,problemStatus_,forceFactorization_,dontFactorizePivots_); #endif int saveType = type; double lastSumDualInfeasibilities = sumDualInfeasibilities_; double lastSumPrimalInfeasibilities = sumPrimalInfeasibilities_; if (type >= 4) { // force factorization type &= 3; numberPivots = 9999999; } //double realDualInfeasibilities = 0.0; double changeCost; bool unflagVariables = numberFlagged_ > 0; int weightsSaved = 0; int numberBlackMarks = 0; if (!type) { // be optimistic stateOfProblem_ &= ~PESSIMISTIC; // but use 0.1 double newTolerance = CoinMax(0.1, saveData_.pivotTolerance_); abcFactorization_->pivotTolerance(newTolerance); } if (type != 3 && (numberPivots > dontFactorizePivots_ || numberIterations_ == baseIteration_)) { // factorize // save dual weights abcDualRowPivot_->saveWeights(this, 1); weightsSaved = 2; // is factorization okay? int factorizationStatus = 0; if ((largestPrimalError_ > 1.0e1 || largestDualError_ > 1.0e1) && numberRows_ > 100 && numberIterations_ > baseIteration_) { double newTolerance = CoinMin(1.1 * abcFactorization_->pivotTolerance(), 0.99); abcFactorization_->pivotTolerance(newTolerance); } #ifdef EARLY_FACTORIZE if ((numberEarly_ & 0xffff) <= 5 || !usefulArray_[ABC_NUMBER_USEFUL - 1].getNumElements()) { //#define KEEP_TRYING_EARLY #ifdef KEEP_TRYING_EARLY int savedEarly = numberEarly_ >> 16; int newEarly = (numberEarly_ & 0xffff); newEarly = CoinMin(savedEarly, 1 + 2 * newEarly); numberEarly_ = (savedEarly << 16) | newEarly; #endif factorizationStatus = internalFactorize(type ? 1 : 2); } else { CoinIndexedVector &vector = usefulArray_[ABC_NUMBER_USEFUL - 1]; bool bad = vector.getNumElements() < 0; const int *indices = vector.getIndices(); double *array = vector.denseVector(); int capacity = vector.capacity(); capacity--; int numberPivotsStored = indices[capacity]; #if ABC_NORMAL_DEBUG > 0 printf("early status %s nPivots %d\n", bad ? "bad" : "good", numberPivotsStored); #endif if (!bad) { // do remaining pivots int iterationsDone = abcEarlyFactorization_->pivots(); factorizationStatus = abcEarlyFactorization_->replaceColumns(this, vector, iterationsDone, numberPivotsStored, true); if (!factorizationStatus) { // should be able to compare e.g. basics 1.0 #define GOOD_EARLY #ifdef GOOD_EARLY AbcSimplexFactorization *temp = abcFactorization_; abcFactorization_ = abcEarlyFactorization_; abcEarlyFactorization_ = temp; #endif moveToBasic(); #ifndef GOOD_EARLY int *pivotSave = CoinCopyOfArray(abcPivotVariable_, numberRows_); factorizationStatus = internalFactorize(type ? 1 : 2); { double *array = usefulArray_[3].denseVector(); double *array2 = usefulArray_[4].denseVector(); for (int iPivot = 0; iPivot < numberRows_; iPivot++) { int iSequence = abcPivotVariable_[iPivot]; unpack(usefulArray_[3], iSequence); unpack(usefulArray_[4], iSequence); abcFactorization_->updateColumn(usefulArray_[3]); abcEarlyFactorization_->updateColumn(usefulArray_[4]); assert(fabs(array[iPivot] - 1.0) < 1.0e-4); array[iPivot] = 0.0; int iPivot2; for (iPivot2 = 0; iPivot2 < numberRows_; iPivot2++) { if (pivotSave[iPivot2] == iSequence) break; } assert(iPivot2 < numberRows_); assert(fabs(array2[iPivot2] - 1.0) < 1.0e-4); array2[iPivot2] = 0.0; for (int i = 0; i < numberRows_; i++) assert(fabs(array[i]) < 1.0e-4); for (int i = 0; i < numberRows_; i++) assert(fabs(array2[i]) < 1.0e-4); usefulArray_[3].clear(); usefulArray_[4].clear(); } delete[] pivotSave; } #endif } else { bad = true; } } // clean up vector.setNumElements(0); for (int i = 0; i < 2 * numberPivotsStored; i++) { array[--capacity] = 0.0; } if (bad) { // switch off #if ABC_NORMAL_DEBUG > 0 printf("switching off early factorization(in statusOf)\n"); #endif #ifndef KEEP_TRYING_EARLY numberEarly_ = 0; delete abcEarlyFactorization_; abcEarlyFactorization_ = NULL; #else numberEarly_ &= 0xffff0000; #endif factorizationStatus = internalFactorize(1); } } #else factorizationStatus = internalFactorize(type ? 1 : 2); #endif if (factorizationStatus) { if (type == 1 && lastFlaggedIteration_ != 1) { // use for real disaster lastFlaggedIteration_ = 1; // no - restore previous basis unflagVariables = false; // restore weights weightsSaved = 4; changeMade_++; // say something changed // Keep any flagged variables for (int i = 0; i < numberTotal_; i++) { if (flagged(i)) internalStatusSaved_[i] |= 64; //say flagged } restoreGoodStatus(type); // get correct bounds on all variables resetFakeBounds(1); numberBlackMarks = 10; if (sequenceOut_ >= 0) { // need to reject something char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 3 printf("flag d\n"); #endif setFlagged(sequenceOut_); lastBadIteration_ = numberIterations_; // say be more cautious // so can't happen immediately again sequenceOut_ = -1; } //abcProgress_.clearBadTimes(); // Go to safe abcFactorization_->pivotTolerance(0.99); } else { #if ABC_NORMAL_DEBUG > 0 printf("Bad initial basis\n"); #endif } if (internalFactorize(1)) { // try once more if (internalFactorize(0)) abort(); } #if PARTITION_ROW_COPY == 1 int iVector = getAvailableArray(); abcMatrix_->sortUseful(usefulArray_[iVector]); setAvailableArray(iVector); #endif } else { lastFlaggedIteration_ = 0; } } // get primal and dual solutions numberBlackMarks += gutsOfSolution(3) / 2; assert(type != 0 || numberIterations_ == baseIteration_); // may be wrong if (!type) { initialSumInfeasibilities_ = sumPrimalInfeasibilities_; initialNumberInfeasibilities_ = numberPrimalInfeasibilities_; CoinAbcMemcpy(solutionSaved_ + numberTotal_, abcSolution_, numberTotal_); CoinAbcMemcpy(internalStatusSaved_ + numberTotal_, internalStatus_, numberTotal_); } double savePrimalInfeasibilities = sumPrimalInfeasibilities_; if ((moreSpecialOptions_ & 16384) != 0 && (moreSpecialOptions_ & 8192) == 0) { if (numberIterations_ == baseIteration_) { // If primal feasible and looks suited to primal go to primal if (!sumPrimalInfeasibilities_ && sumDualInfeasibilities_ > 1.0 && numberRows_ * 4 > numberColumns_) { // do primal // mark so can perturb initialSumInfeasibilities_ = -1.0; problemStatus_ = 10; return; } } else if ((numberIterations_ > 1000 && numberIterations_ * 10 > numberRows_ && numberIterations_ * 4 < numberRows_) || (sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e10 * (initialSumInfeasibilities_ + 1000.0) && numberIterations_ > 10000 && sumDualInfeasibilities_ > 1.0)) { if (sumPrimalInfeasibilities_ + 1.0e2 * sumDualInfeasibilities_ > 2000.0 * (initialSumInfeasibilities_ + 100.0)) { #if ABC_NORMAL_DEBUG > 0 printf("sump %g sumd %g initial %g\n", sumPrimalInfeasibilities_, sumDualInfeasibilities_, initialSumInfeasibilities_); #endif // but only if average gone up as well if (sumPrimalInfeasibilities_ * initialNumberInfeasibilities_ > 100.0 * (initialSumInfeasibilities_ + 10.0) * numberPrimalInfeasibilities_) { #if ABC_NORMAL_DEBUG > 0 printf("returning to do primal\n"); #endif // do primal problemStatus_ = 10; CoinAbcMemcpy(abcSolution_, solutionSaved_ + numberTotal_, numberTotal_); CoinAbcMemcpy(internalStatus_, internalStatusSaved_ + numberTotal_, numberTotal_); // mark so can perturb initialSumInfeasibilities_ = -1.0; return; } } } } else if ((sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e6 * (initialSumInfeasibilities_ + 1000.0) && (moreSpecialOptions_ & 8192) == 0 && numberIterations_ > 1000 && numberIterations_ * 10 > numberRows_ && numberIterations_ * 4 < numberRows_) || (sumPrimalInfeasibilities_ + sumDualInfeasibilities_ > 1.0e10 * (initialSumInfeasibilities_ + 1000.0) && (moreSpecialOptions_ & 8192) == 0 && numberIterations_ > 10000 && sumDualInfeasibilities_ > 1.0)) { // do primal problemStatus_ = 10; CoinAbcMemcpy(abcSolution_, solutionSaved_ + numberTotal_, numberTotal_); CoinAbcMemcpy(internalStatus_, internalStatusSaved_ + numberTotal_, numberTotal_); // mark so can perturb initialSumInfeasibilities_ = -1.0; return; } //double saveDualInfeasibilities=sumDualInfeasibilities_; // 0 don't ignore, 1 ignore errors, 2 values pass int ignoreStuff = (firstFree_ < 0 && lastFirstFree_ < 0) ? 0 : 1; if ((stateOfProblem_ & VALUES_PASS) != 0) ignoreStuff = 2; { double thisObj = rawObjectiveValue(); double lastObj = abcProgress_.lastObjective(0); if ((lastObj > thisObj + 1.0e-3 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4 || (sumPrimalInfeasibilities_ > 10.0 * lastSumPrimalInfeasibilities + 1.0e3 && sumDualInfeasibilities_ > 10.0 * CoinMax(lastSumDualInfeasibilities, 1.0))) && !ignoreStuff) { #if ABC_NORMAL_DEBUG > 0 printf("bad - objective backwards %g %g errors %g %g suminf %g %g\n", lastObj, thisObj, largestDualError_, largestPrimalError_, sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif numberBlackMarks++; if (largestDualError_ > 1.0e-3 || largestPrimalError_ > 1.0e-3 || (sumDualInfeasibilities_ > 10.0 * CoinMax(lastSumDualInfeasibilities, 1.0) && firstFree_ < 0)) { #if ABC_NORMAL_DEBUG > 0 printf("backtracking %d iterations - dual error %g primal error %g - sumDualInf %g\n", numberPivots, largestDualError_, largestPrimalError_, sumDualInfeasibilities_); #endif // restore previous basis unflagVariables = false; // restore weights weightsSaved = 4; changeMade_++; // say something changed // Keep any flagged variables for (int i = 0; i < numberTotal_; i++) { if (flagged(i)) internalStatusSaved_[i] |= 64; //say flagged } restoreGoodStatus(type); // get correct bounds on all variables resetFakeBounds(1); numberBlackMarks += 10; if (sequenceOut_ >= 0) { // need to reject something #if ABC_NORMAL_DEBUG > 0 if (flagged(sequenceOut_)) printf("BAD %x already flagged\n", sequenceOut_); #endif char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 3 printf("flag d\n"); #endif setFlagged(sequenceOut_); lastBadIteration_ = numberIterations_; // say be more cautious #if ABC_NORMAL_DEBUG > 0 } else { printf("backtracking - no valid sequence out\n"); #endif } // so can't happen immediately again sequenceOut_ = -1; //abcProgress_.clearBadTimes(); // Go to safe //abcFactorization_->pivotTolerance(0.99); if (internalFactorize(1)) { abort(); } moveToBasic(); #if PARTITION_ROW_COPY == 1 int iVector = getAvailableArray(); abcMatrix_->sortUseful(usefulArray_[iVector]); setAvailableArray(iVector); #endif // get primal and dual solutions numberBlackMarks += gutsOfSolution(3); //gutsOfSolution(1); makeNonFreeVariablesDualFeasible(); } } } //#define PAN #ifdef PAN if (!type && (specialOptions_ & COIN_CBC_USING_CLP) == 0 && firstFree_ < 0) { // first time - set large bad ones superbasic int whichArray = getAvailableArray(); double *COIN_RESTRICT work = usefulArray_[whichArray].denseVector(); int number = 0; int *COIN_RESTRICT which = usefulArray_[whichArray].getIndices(); const double *COIN_RESTRICT reducedCost = abcDj_; const double *COIN_RESTRICT lower = abcLower_; const double *COIN_RESTRICT upper = abcUpper_; #if PRINT_PAN int nAtUb = 0; #endif for (int i = 0; i < numberTotal_; i++) { Status status = getInternalStatus(i); double djValue = reducedCost[i]; double badValue = 0.0; if (status == atLowerBound && djValue < -currentDualTolerance_) { double gap = upper[i] - lower[i]; if (gap > 100000000.0) badValue = -djValue * CoinMin(gap, 1.0e10); } else if (status == atUpperBound && djValue > currentDualTolerance_) { double gap = upper[i] - lower[i]; if (gap > 100000000.0) { badValue = djValue * CoinMin(gap, 1.0e10); #if PRINT_PAN nAtUb++; #endif } } if (badValue) { work[number] = badValue; which[number++] = i; } } // for now don't do if primal feasible (probably should go to primal) if (!numberPrimalInfeasibilities_) { number = 0; #if ABC_NORMAL_DEBUG > 1 printf("XXXX doesn't work if primal feasible - also fix switch to Clp primal\n"); #endif } if (number && ignoreStuff != 2) { CoinSort_2(work, work + number, which); int numberToDo = CoinMin(number, numberRows_ / 10); // ?? #if PRINT_PAN > 1 CoinSort_2(which, which + numberToDo, work); int n = 0; #endif #if PRINT_PAN printf("Panning %d variables (out of %d at lb and %d at ub)\n", numberToDo, number - nAtUb, nAtUb); #endif firstFree_ = numberTotal_; for (int i = 0; i < numberToDo; i++) { int iSequence = which[i]; #if PRINT_PAN > 1 if (!n) printf("Pan "); printf("%d ", iSequence); n++; if (n == 10) { n = 0; printf("\n"); } #endif setInternalStatus(iSequence, superBasic); firstFree_ = CoinMin(firstFree_, iSequence); } #if PRINT_PAN > 1 if (n) printf("\n"); #endif ordinaryVariables_ = 0; CoinAbcMemset0(work, number); stateOfProblem_ |= FAKE_SUPERBASIC; } setAvailableArray(whichArray); } #endif if (!sumOfRelaxedPrimalInfeasibilities_ && sumDualInfeasibilities_ && type && saveType < 9 && ignoreStuff != 2) { // say been feasible #ifdef ABC_INHERIT abcState_ |= CLP_ABC_BEEN_FEASIBLE; #endif int numberFake = numberAtFakeBound(); if (!numberFake) { makeNonFreeVariablesDualFeasible(); numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; // just primal gutsOfSolution(2); } } bool computeNewSumInfeasibilities = numberIterations_ == baseIteration_; if ((!type || (firstFree_ < 0 && lastFirstFree_ >= 0)) && ignoreStuff != 2) { #ifdef PAN if ((stateOfProblem_ & FAKE_SUPERBASIC) != 0 && firstFree_ < 0 && lastFirstFree_ >= 0 && type) { // clean up #if PRINT_PAN printf("Pan switch off after %d iterations\n", numberIterations_); #endif computeNewSumInfeasibilities = true; ordinaryVariables_ = 1; const double *COIN_RESTRICT reducedCost = abcDj_; const double *COIN_RESTRICT lower = abcLower_; const double *COIN_RESTRICT upper = abcUpper_; for (int i = 0; i < numberTotal_; i++) { Status status = getInternalStatus(i); if (status == superBasic) { double djValue = reducedCost[i]; double value = abcSolution_[i]; if (value < lower[i] + primalTolerance_) { setInternalStatus(i, atLowerBound); #if ABC_NORMAL_DEBUG > 0 if (djValue < -dualTolerance_) printf("Odd at lb %d dj %g\n", i, djValue); #endif } else if (value > upper[i] - primalTolerance_) { setInternalStatus(i, atUpperBound); #if ABC_NORMAL_DEBUG > 0 if (djValue > dualTolerance_) printf("Odd at ub %d dj %g\n", i, djValue); #endif } else { assert(status != superBasic); } } if (status == isFree) ordinaryVariables_ = 0; } stateOfProblem_ &= ~FAKE_SUPERBASIC; } #endif // set up perturbation and dualBound_ // make dual feasible if (firstFree_ < 0) { if (numberFlagged_) { numberFlagged_ = 0; for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; if (flagged(iPivot)) { clearFlagged(iPivot); } } } int numberChanged = resetFakeBounds(0); #if ABC_NORMAL_DEBUG > 0 if (numberChanged > 0) printf("Re-setting %d fake bounds %d gone dual infeasible - can expect backward objective\n", numberChanged, numberDualInfeasibilities_); #endif if (numberDualInfeasibilities_ || numberChanged != -1) { makeNonFreeVariablesDualFeasible(); numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; sumOfRelaxedDualInfeasibilities_ = 0.0; // just primal gutsOfSolution(2); if (!numberIterations_ && sumPrimalInfeasibilities_ > 1.0e5 * (savePrimalInfeasibilities + 1.0e3) && (moreSpecialOptions_ & (256 | 8192)) == 0) { // Use primal problemStatus_ = 10; // mark so can perturb initialSumInfeasibilities_ = -1.0; CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_); CoinAbcMemcpy(internalStatus_, internalStatusSaved_, numberTotal_); return; } // redo if (computeNewSumInfeasibilities) initialSumInfeasibilities_ = sumPrimalInfeasibilities_; //computeObjective(); } } //bounceTolerances(type); } else if (numberIterations_ > numberRows_ && stateDualColumn_ == -2) { bounceTolerances(4); } // If bad accuracy treat as singular if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_ > 1000) { problemStatus_ = 10; //abort(); #if ABC_NORMAL_DEBUG > 0 printf("Big problems - errors %g %g try primal\n", largestPrimalError_, largestDualError_); #endif } else if (goodAccuracy() && numberIterations_ > lastBadIteration_ + 200) { // Can reduce tolerance double newTolerance = CoinMax(0.995 * abcFactorization_->pivotTolerance(), saveData_.pivotTolerance_); if (!type) newTolerance = saveData_.pivotTolerance_; if (newTolerance > abcFactorization_->minimumPivotTolerance()) abcFactorization_->pivotTolerance(newTolerance); } bestObjectiveValue_ = CoinMax(bestObjectiveValue_, rawObjectiveValue_ - bestPossibleImprovement_); // Double check infeasibility if no action if (abcProgress_.lastIterationNumber(0) == numberIterations_) { if (abcDualRowPivot_->looksOptimal()) { // say been feasible #ifdef ABC_INHERIT abcState_ |= CLP_ABC_BEEN_FEASIBLE; #endif numberPrimalInfeasibilities_ = 0; sumPrimalInfeasibilities_ = 0.0; } } else { double thisObj = rawObjectiveValue_ - bestPossibleImprovement_; #ifdef CLP_INVESTIGATE assert(bestPossibleImprovement_ > -1000.0 && rawObjectiveValue_ > -1.0e100); if (bestPossibleImprovement_) printf("obj %g add in %g -> %g\n", rawObjectiveValue_, bestPossibleImprovement_, thisObj); #endif double lastObj = abcProgress_.lastObjective(0); #ifndef NDEBUG #if ABC_NORMAL_DEBUG > 3 resetFakeBounds(-1); #endif #endif if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3 && saveType < 9) { numberTimesOptimal_ = 0; } } // Up tolerance if looks a bit odd if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) { if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) { int backIteration = abcProgress_.lastIterationNumber(CLP_PROGRESS - 1); if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) { if (abcFactorization_->pivotTolerance() < 0.9) { // up tolerance abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91)); //printf("tol now %g\n",abcFactorization_->pivotTolerance()); abcProgress_.clearIterationNumbers(); } } } } // Check if looping int loop; if (type != 2) loop = abcProgress_.looping(); else loop = -1; if (abcProgress_.reallyBadTimes() > 10) { problemStatus_ = 10; // instead - try other algorithm #if ABC_NORMAL_DEBUG > 3 printf("returning at %d\n", __LINE__); #endif } if (loop >= 0) { problemStatus_ = loop; //exit if in loop if (!problemStatus_) { // declaring victory // say been feasible #ifdef ABC_INHERIT abcState_ |= CLP_ABC_BEEN_FEASIBLE; #endif numberPrimalInfeasibilities_ = 0; sumPrimalInfeasibilities_ = 0.0; } else { problemStatus_ = 10; // instead - try other algorithm #if ABC_NORMAL_DEBUG > 3 printf("returning at %d\n", __LINE__); #endif } return; } // really for free variables in if ((progressFlag_ & 2) != 0) { //situationChanged = 2; } progressFlag_ &= (~3); //reset progress flag // mark as having gone optimal if looks like it if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilities_) progressFlag_ |= 8; if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) { handler_->message(CLP_SIMPLEX_STATUS, messages_) << numberIterations_ << objectiveValue(); handler_->printing(sumPrimalInfeasibilities_ > 0.0) << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_; handler_->printing(sumDualInfeasibilities_ > 0.0) << sumDualInfeasibilities_ << numberDualInfeasibilities_; handler_->printing(numberDualInfeasibilitiesWithoutFree_ < numberDualInfeasibilities_) << numberDualInfeasibilitiesWithoutFree_; handler_->message() << CoinMessageEol; } double approximateObjective = rawObjectiveValue_; //realDualInfeasibilities = sumDualInfeasibilities_; // If we need to carry on cleaning variables if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) { for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; if (!flagged(iPivot) && pivoted(iPivot)) { // carry on numberPrimalInfeasibilities_ = -1; sumOfRelaxedPrimalInfeasibilities_ = 1.0; sumPrimalInfeasibilities_ = 1.0; break; } } } /* If we are primal feasible and any dual infeasibilities are on free variables then it is better to go to primal */ if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ && numberDualInfeasibilities_) { // say been feasible #ifdef ABC_INHERIT abcState_ |= CLP_ABC_BEEN_FEASIBLE; #endif if (type && ignoreStuff != 2) problemStatus_ = 10; else numberPrimalInfeasibilities_ = 1; } // check optimal // give code benefit of doubt if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0 && ignoreStuff != 2) { // say been feasible #ifdef ABC_INHERIT abcState_ |= CLP_ABC_BEEN_FEASIBLE; #endif // say optimal (with these bounds etc) numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; numberPrimalInfeasibilities_ = 0; sumPrimalInfeasibilities_ = 0.0; } else if (sumOfRelaxedDualInfeasibilities_ > 1.0e-3 && firstFree_ < 0 && ignoreStuff != 2 && numberIterations_ > lastCleaned_ + 10) { #if ABC_NORMAL_DEBUG > 0 printf("Making dual feasible\n"); #endif makeNonFreeVariablesDualFeasible(true); lastCleaned_ = numberIterations_; numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; // just primal gutsOfSolution(2); } // see if cutoff reached double limit = 0.0; getDblParam(ClpDualObjectiveLimit, limit); int numberChangedBounds = 0; if (primalFeasible() && ignoreStuff != 2) { // may be optimal - or may be bounds are wrong numberChangedBounds = changeBounds(0, changeCost); if (numberChangedBounds) gutsOfSolution(2); } if (primalFeasible() && ignoreStuff != 2) { if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) { //looks optimal - do we need to reset tolerance if (lastCleaned_ < numberIterations_ && /*numberTimesOptimal_*/ stateDualColumn_ < 4 && (numberChanged_ || (specialOptions_ & 4096) == 0)) { #if CLP_CAN_HAVE_ZERO_OBJ if ((specialOptions_ & 2097152) == 0) { #endif numberTimesOptimal_++; if (stateDualColumn_ == -2) { #if ABC_NORMAL_DEBUG > 0 printf("Going straight from state -2 to 0\n"); #endif stateDualColumn_ = -1; } changeMade_++; // say something changed if (numberTimesOptimal_ == 1) { if (fabs(currentDualTolerance_ - dualTolerance_) > 1.0e-12) { #if ABC_NORMAL_DEBUG > 0 printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n", currentDualTolerance_, dualTolerance_, numberTimesOptimal_); #endif currentDualTolerance_ = dualTolerance_; } } else { if (numberTimesOptimal_ == 2) { // better to have small tolerance even if slower abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15)); } #if ABC_NORMAL_DEBUG > 0 printf("changing dual tolerance from %g to %g (numberTimesOptimal_==%d)\n", currentDualTolerance_, dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1), numberTimesOptimal_); #endif currentDualTolerance_ = dualTolerance_ * pow(2.0, numberTimesOptimal_ - 1); } if (numberChanged_ > 0) { handler_->message(CLP_DUAL_ORIGINAL, messages_) << CoinMessageEol; //realDualInfeasibilities = sumDualInfeasibilities_; } else if (numberChangedBounds < -1) { // bounds were flipped //gutsOfSolution(2); } if (stateDualColumn_ >= 0) stateDualColumn_++; #ifndef TRY_ABC_GUS int returnCode = bounceTolerances(1); if (returnCode) #endif // always go to primal problemStatus_ = 10; lastCleaned_ = numberIterations_; #if CLP_CAN_HAVE_ZERO_OBJ } else { // no cost - skip checks problemStatus_ = 0; } #endif } else { problemStatus_ = 0; // optimal if (lastCleaned_ < numberIterations_ && numberChanged_) { handler_->message(CLP_SIMPLEX_GIVINGUP, messages_) << CoinMessageEol; } } } else if (numberChangedBounds < -1) { // some flipped bounceTolerances(2); #if ABC_NORMAL_DEBUG > 0 printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n", numberChangedBounds, numberAtFakeBound(), __LINE__, __FILE__); #endif } else if (numberChangedBounds > 0 || numberAtFakeBound()) { handler_->message(CLP_DUAL_BOUNDS, messages_) << currentDualBound_ << CoinMessageEol; #if ABC_NORMAL_DEBUG > 0 printf("changing bounds\n"); printf("numberChangedBounds %d numberAtFakeBound %d at line %d in file %s\n", numberChangedBounds, numberAtFakeBound(), __LINE__, __FILE__); #endif if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_) { // check unbounded // find a variable with bad dj int iChosen = -1; if ((specialOptions_ & 0x03000000) == 0) { double largest = 100.0 * primalTolerance_; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { double djValue = abcDj_[iSequence]; double originalLo = lowerSaved_[iSequence]; double originalUp = upperSaved_[iSequence]; if (fabs(djValue) > fabs(largest)) { if (getInternalStatus(iSequence) != basic) { if (djValue > 0 && originalLo < -1.0e20) { if (djValue > fabs(largest)) { largest = djValue; iChosen = iSequence; } } else if (djValue < 0 && originalUp > 1.0e20) { if (-djValue > fabs(largest)) { largest = djValue; iChosen = iSequence; } } } } } } if (iChosen >= 0) { int iVector = getAvailableArray(); unpack(usefulArray_[iVector], iChosen); //double changeCost; problemStatus_ = checkUnbounded(usefulArray_[iVector], 0.0 /*changeCost*/); usefulArray_[iVector].clear(); setAvailableArray(iVector); } } gutsOfSolution(3); //bounceTolerances(2); //realDualInfeasibilities = sumDualInfeasibilities_; } // unflag all variables (we may want to wait a bit?) if (unflagVariables) { int numberFlagged = numberFlagged_; numberFlagged_ = 0; for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; if (flagged(iPivot)) { clearFlagged(iPivot); } } #if ABC_NORMAL_DEBUG > 3 if (numberFlagged) { printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus, problemStatus_, numberPrimalInfeasibilities_, numberTimesOptimal_); } #endif unflagVariables = numberFlagged > 0; if (numberFlagged && !numberPivots) { /* looks like trouble as we have not done any iterations. Try changing pivot tolerance then give it a few goes and give up */ if (abcFactorization_->pivotTolerance() < 0.9) { abcFactorization_->pivotTolerance(0.99); problemStatus_ = -1; } else if (numberTimesOptimal_ < 3) { numberTimesOptimal_++; problemStatus_ = -1; } else { unflagVariables = false; //secondaryStatus_ = 1; // and say probably infeasible #ifdef ABC_INHERIT if ((moreSpecialOptions_ & 256) == 0 || (abcState_ & CLP_ABC_BEEN_FEASIBLE) != 0) #else if ((moreSpecialOptions_ & 256) == 0) #endif { // try primal problemStatus_ = 10; } else { // almost certainly infeasible problemStatus_ = 1; } #if ABC_NORMAL_DEBUG > 3 printf("returning at %d\n", __LINE__); #endif } } } } if (problemStatus_ < 0) { saveGoodStatus(); if (weightsSaved) { // restore weights (if saved) - also recompute infeasibility list abcDualRowPivot_->saveWeights(this, weightsSaved); } else { abcDualRowPivot_->recomputeInfeasibilities(); } } // see if cutoff reached checkCutoff(false); // If we are in trouble and in branch and bound give up if ((specialOptions_ & 1024) != 0) { int looksBad = 0; if (largestPrimalError_ * largestDualError_ > 1.0e2) { looksBad = 1; } else if (largestPrimalError_ > 1.0e-2 && rawObjectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) { looksBad = 2; } if (looksBad) { if (abcFactorization_->pivotTolerance() < 0.9) { // up tolerance abcFactorization_->pivotTolerance(CoinMin(abcFactorization_->pivotTolerance() * 1.05 + 0.02, 0.91)); } else if (numberIterations_ > 10000) { #if ABC_NORMAL_DEBUG > 0 if (handler_->logLevel() > 2) printf("bad dual - saying infeasible %d\n", looksBad); #endif problemStatus_ = 1; secondaryStatus_ = 1; // and say was on cutoff } else if (largestPrimalError_ > 1.0e5) { #if ABC_NORMAL_DEBUG > 0 if (handler_->logLevel() > 2) printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_); #endif allSlackBasis(); problemStatus_ = 10; } } } //if (problemStatus_ < 0 && !changeMade_) { //problemStatus_ = 4; // unknown //} lastGoodIteration_ = numberIterations_; if (numberIterations_ > lastBadIteration_ + 200) { moreSpecialOptions_ &= ~16; // clear check accuracy flag if (numberFlagged_) { numberFlagged_ = 0; for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; if (flagged(iPivot)) { clearFlagged(iPivot); } } } } if (problemStatus_ < 0 && ignoreStuff != 2) { //sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful if (sumDualInfeasibilities_) { numberDualInfeasibilities_ = 1; } else if (!numberPrimalInfeasibilities_) { // Looks good problemStatus_ = 0; } } double thisObj = abcProgress_.lastObjective(0); double lastObj = abcProgress_.lastObjective(1); if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4 && lastFirstFree_ < 0) { //printf("BAD - objective backwards\n"); //bounceTolerances(3); numberBlackMarks += 3; } if (numberBlackMarks > 0 && numberIterations_ > baseIteration_) { // be very careful int maxFactor = abcFactorization_->maximumPivots(); if (maxFactor > 10) { if ((stateOfProblem_ & PESSIMISTIC) == 0 || true) { if (largestDualError_ > 10.0 * CoinMax(lastDualError_, 1.0e-6) || largestPrimalError_ > 10.0 * CoinMax(lastPrimalError_, 1.0e-6)) maxFactor = CoinMin(maxFactor, 20); forceFactorization_ = CoinMax(1, (maxFactor >> 1)); } else if (numberBlackMarks > 2) { forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1)); } } stateOfProblem_ |= PESSIMISTIC; if (numberBlackMarks >= 10) forceFactorization_ = 1; } lastPrimalError_ = largestPrimalError_; lastDualError_ = largestDualError_; lastFirstFree_ = firstFree_; if (ignoreStuff == 2) { // in values pass lastFirstFree_ = -1; // firstFree_=-1; } if (alphaAccuracy_ > 0.0) alphaAccuracy_ = 1.0; // If we are stopping - use plausible objective // Maybe only in fast dual if (problemStatus_ > 2) rawObjectiveValue_ = approximateObjective; if (problemStatus_ == 1 && (progressFlag_ & 8) != 0 && fabs(rawObjectiveValue_) > 1.0e10) problemStatus_ = 10; // infeasible - but has looked feasible checkMoveBack(true); #if 0 { double * array = usefulArray_[arrayForTableauRow_].denseVector(); for (int iPivot = 0; iPivot < numberRows_; iPivot++) { int iSequence = abcPivotVariable_[iPivot]; unpack(usefulArray_[arrayForTableauRow_],iSequence); abcFactorization_->updateColumn(usefulArray_[arrayForTableauRow_]); assert (fabs(array[iPivot] - 1.0) < 1.0e-4); array[iPivot] = 0.0; for (int i = 0; i < numberRows_; i++) assert (fabs(array[i]) < 1.0e-4); usefulArray_[arrayForTableauRow_].clear(); } } #endif #if 0 { // save status, cost and solution unsigned char * saveStatus = CoinCopyOfArray(internalStatus_,numberTotal_); double * saveSolution = CoinCopyOfArray(abcSolution_,numberTotal_); double * saveCost = CoinCopyOfArray(abcCost_,numberTotal_); printf("TestA sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_); if (perturbation_>100) CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_); else CoinAbcMemcpy(abcCost_,abcPerturbation_,numberTotal_); moveToBasic(); gutsOfSolution(3); for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { bool changed=fabs(saveCost[iSequence]-abcCost_[iSequence])>1.0e-10; bool bad=false; int iStatus=internalStatus_[iSequence]&7; if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) { // not basic double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence]; double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence]; double value = abcDj_[iSequence]; if (distanceUp > primalTolerance_) { // Check if "free" if (distanceDown <= primalTolerance_) { // should not be negative if (value < 0.0) { value = - value; if (value > currentDualTolerance_) { bad=true; } } } else { // free value=fabs(value); if (value > currentDualTolerance_) { bad=true; } } } else if (distanceDown > primalTolerance_) { // should not be positive if (value > 0.0) { if (value > currentDualTolerance_) { bad=true; } } } } if (changed||bad) printf("seq %d status %d current cost %g original %g dj %g %s\n", iSequence,iStatus,saveCost[iSequence],abcCost_[iSequence], abcDj_[iSequence],bad ? "(bad)" : ""); } printf("TestB sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_); if (numberDualInfeasibilities_) { makeNonFreeVariablesDualFeasible(); moveToBasic(); gutsOfSolution(3); printf("TestC sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_); } CoinAbcMemcpy(internalStatus_,saveStatus,numberTotal_); CoinAbcMemcpy(abcSolution_,saveSolution,numberTotal_); CoinAbcMemcpy(abcCost_,saveCost,numberTotal_); delete [] saveStatus; delete [] saveSolution; delete [] saveCost; moveToBasic(); gutsOfSolution(3); printf("TestD sumDual %g sumPrimal %g\n",sumDualInfeasibilities_,sumPrimalInfeasibilities_); } #endif #if 0 { const double multiplier[] = { 1.0, -1.0}; // sign????????????? double dualT = - currentDualTolerance_; double largestCostDifference=0.1*currentDualTolerance_; int numberCostDifference=0; double largestBadDj=0.0; int numberBadDj=0; double sumBadDj=0.0; double sumCanMoveBad=0.0; int numberCanMoveBad=0; double sumCanMoveOthers=0.0; int numberCanMoveOthers=0; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { double costDifference=abcPerturbation_[iSequence]-abcCost_[iSequence]; if (fabs(costDifference)>0.1*currentDualTolerance_) { if (fabs(costDifference)>largestCostDifference) largestCostDifference=fabs(costDifference); numberCostDifference++; } unsigned char iStatus = internalStatus_[iSequence] & 7; if (iStatus<4) { double mult = multiplier[iStatus]; double dj = abcDj_[iSequence] * mult; if (dj < dualT) { sumBadDj+=dualT-dj; numberBadDj++; if (dualT-dj>largestBadDj) largestBadDj=dualT-dj; if (costDifference*mult>0.0) { sumCanMoveBad+=costDifference*mult; numberCanMoveBad++; //dj+=costDifference*mult; abcDj_[iSequence]+=costDifference; abcCost_[iSequence]=abcPerturbation_[iSequence]; } } else if (costDifference*mult>0.0) { sumCanMoveOthers+=costDifference*mult; numberCanMoveOthers++; //dj+=costDifference*mult; abcDj_[iSequence]+=costDifference; abcCost_[iSequence]=abcPerturbation_[iSequence]; } } } if (numberCostDifference+numberBadDj+numberCanMoveBad+numberCanMoveOthers) { printf("largest diff %g (%d), largest bad dj %g (%d - sum %g), can move bad %g (%d), good %g (%d) - %d at fake bound\n",largestCostDifference,numberCostDifference,largestBadDj, numberBadDj, sumBadDj,sumCanMoveBad=0.0,numberCanMoveBad, sumCanMoveOthers,numberCanMoveOthers,numberAtFakeBound()); } } #endif } // see if cutoff reached bool AbcSimplexDual::checkCutoff(bool computeObjective) { double objValue = COIN_DBL_MAX; if (computeObjective) objValue = computeInternalObjectiveValue(); else objValue = minimizationObjectiveValue(); double limit = 0.0; getDblParam(ClpDualObjectiveLimit, limit); // Looks as if numberAtFakeBound() being computed too often if (fabs(limit) < 1.0e30 && objValue > limit && !numberAtFakeBound()) { bool looksInfeasible = !numberDualInfeasibilities_; if (objValue > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 && sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000) looksInfeasible = true; if (looksInfeasible && !computeObjective) objValue = computeInternalObjectiveValue(); if (looksInfeasible) { // Even if not perturbed internal costs may have changed // be careful if (objValue > limit) { problemStatus_ = 1; secondaryStatus_ = 1; // and say was on cutoff return true; } } } return false; } #define CLEANUP_DJS 0 /* While updateDualsInDual sees what effect is of flip this does actual flipping. If change >0.0 then value in array >0.0 => from lower to upper returns 3 if skip this iteration and re-factorize */ int AbcSimplexDual::flipBounds() { CoinIndexedVector &array = usefulArray_[arrayForFlipBounds_]; CoinIndexedVector &output = usefulArray_[arrayForFlipRhs_]; output.clear(); int numberFlipped = array.getNumElements(); if (numberFlipped) { #if CLP_CAN_HAVE_ZERO_OBJ > 1 if ((specialOptions_ & 2097152) == 0) { #endif // do actual flips // do in parallel once got spare space in factorization #if ABC_PARALLEL == 2 #endif int number = array.getNumElements(); const int *COIN_RESTRICT which = array.getIndices(); double *COIN_RESTRICT work = array.denseVector(); double *COIN_RESTRICT solution = abcSolution_; const double *COIN_RESTRICT lower = abcLower_; const double *COIN_RESTRICT upper = abcUpper_; #if CLEANUP_DJS double *COIN_RESTRICT cost = abcCost_; double *COIN_RESTRICT dj = abcDj_; // see if we can clean up djs // may have perturbation const double *COIN_RESTRICT originalCost = abcPerturbation_; #else const double *COIN_RESTRICT cost = abcCost_; const double *COIN_RESTRICT dj = abcDj_; #endif const double multiplier[] = { 1.0, -1.0 }; for (int i = 0; i < number; i++) { double value = work[i] * theta_; work[i] = 0.0; int iSequence = which[i]; unsigned char iStatus = internalStatus_[iSequence] & 3; assert((internalStatus_[iSequence] & 7) == 0 || (internalStatus_[iSequence] & 7) == 1); double mult = multiplier[iStatus]; double djValue = dj[iSequence] * mult - value; //double perturbedTolerance = abcPerturbation_[iSequence]; if (djValue < -currentDualTolerance_) { // might just happen - if so just skip assert(iSequence != sequenceIn_); double movement = (upper[iSequence] - lower[iSequence]) * mult; if (iStatus == 0) { // at lower bound // to upper bound setInternalStatus(iSequence, atUpperBound); solution[iSequence] = upper[iSequence]; #if CLEANUP_DJS if (cost[iSequence] < originalCost[iSequence]) { double difference = CoinMax(cost[iSequence] - originalCost[iSequence], djValue); dj[iSequence] -= difference; cost[iSequence] -= difference; } #endif } else { // at upper bound // to lower bound setInternalStatus(iSequence, atLowerBound); solution[iSequence] = lower[iSequence]; #if CLEANUP_DJS if (cost[iSequence] > originalCost[iSequence]) { double difference = CoinMin(cost[iSequence] - originalCost[iSequence], -djValue); dj[iSequence] -= difference; cost[iSequence] -= difference; } #endif } abcMatrix_->add(output, iSequence, movement); objectiveChange_ += cost[iSequence] * movement; } } array.setNumElements(0); // double check something flipped if (output.getNumElements()) { #if ABC_PARALLEL == 0 abcFactorization_->updateColumn(output); #else assert(FACTOR_CPU > 2); abcFactorization_->updateColumnCpu(output, 1); #endif abcDualRowPivot_->updatePrimalSolution(output, 1.0); } #if CLP_CAN_HAVE_ZERO_OBJ > 1 } else { clearArrays(3); } #endif } // recompute dualOut_ valueOut_ = solutionBasic_[pivotRow_]; if (directionOut_ < 0) { dualOut_ = valueOut_ - upperOut_; } else { dualOut_ = lowerOut_ - valueOut_; } #if 0 // amount primal will move movement_ = -dualOut_ * directionOut_ / alpha_; double movementOld = oldDualOut * directionOut_ / alpha_; // so objective should increase by fabs(dj)*movement_ // but we already have objective change - so check will be good if (objectiveChange_ + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(rawObjectiveValue_))) { #if ABC_NORMAL_DEBUG > 3 if (handler_->logLevel() & 32) printf("movement %g, movementOld %g swap change %g, rest %g * %g\n", objectiveChange_ + fabs(movement_ * dualIn_), movementOld,objectiveChange_, movement_, dualIn_); #endif } if (0) { if(abcFactorization_->pivots()) { // going backwards - factorize problemStatus_ = -2; // factorize now stateOfIteration_=2; } } #endif return numberFlipped; } /* Undo a flip */ void AbcSimplexDual::flipBack(int number) { CoinIndexedVector &array = usefulArray_[arrayForFlipBounds_]; const int *COIN_RESTRICT which = array.getIndices(); double *COIN_RESTRICT solution = abcSolution_; const double *COIN_RESTRICT lower = abcLower_; const double *COIN_RESTRICT upper = abcUpper_; for (int i = 0; i < number; i++) { int iSequence = which[i]; unsigned char iStatus = internalStatus_[iSequence] & 3; assert((internalStatus_[iSequence] & 7) == 0 || (internalStatus_[iSequence] & 7) == 1); if (iStatus == 0) { // at lower bound // to upper bound setInternalStatus(iSequence, atUpperBound); solution[iSequence] = upper[iSequence]; } else { // at upper bound // to lower bound setInternalStatus(iSequence, atLowerBound); solution[iSequence] = lower[iSequence]; } } } // Restores bound to original bound void AbcSimplexDual::originalBound(int iSequence) { if (getFakeBound(iSequence) != noFake) { numberFake_--; ; setFakeBound(iSequence, noFake); abcLower_[iSequence] = lowerSaved_[iSequence]; abcUpper_[iSequence] = upperSaved_[iSequence]; } } /* As changeBounds but just changes new bounds for a single variable. Returns true if change */ bool AbcSimplexDual::changeBound(int iSequence) { // old values double oldLower = abcLower_[iSequence]; double oldUpper = abcUpper_[iSequence]; double value = abcSolution_[iSequence]; bool modified = false; // original values double lowerValue = lowerSaved_[iSequence]; double upperValue = upperSaved_[iSequence]; setFakeBound(iSequence, noFake); if (value == oldLower) { if (upperValue > oldLower + currentDualBound_) { abcUpper_[iSequence] = oldLower + currentDualBound_; setFakeBound(iSequence, upperFake); assert(abcLower_[iSequence] == lowerSaved_[iSequence]); assert(abcUpper_[iSequence] < upperSaved_[iSequence]); modified = true; numberFake_++; } } else if (value == oldUpper) { if (lowerValue < oldUpper - currentDualBound_) { abcLower_[iSequence] = oldUpper - currentDualBound_; setFakeBound(iSequence, lowerFake); assert(abcLower_[iSequence] > lowerSaved_[iSequence]); assert(abcUpper_[iSequence] == upperSaved_[iSequence]); modified = true; numberFake_++; } } else { assert(value == oldLower || value == oldUpper); } return modified; } /* Checks number of variables at fake bounds. This is used by fastDual so can exit gracefully before end */ int AbcSimplexDual::numberAtFakeBound() { int numberFake = 0; if (numberFake_) { for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { FakeBound bound = getFakeBound(iSequence); switch (getInternalStatus(iSequence)) { case basic: break; case isFree: case superBasic: case AbcSimplex::isFixed: //setFakeBound (iSequence, noFake); break; case atUpperBound: if (bound == upperFake || bound == bothFake) numberFake++; break; case atLowerBound: if (bound == lowerFake || bound == bothFake) numberFake++; break; } } } return numberFake; } int AbcSimplexDual::resetFakeBounds(int type) { if (type == 0 || type == 1) { // put back original bounds and then check //CoinAbcMemcpy(abcLower_,lowerSaved_,numberTotal_); //CoinAbcMemcpy(abcUpper_,upperSaved_,numberTotal_); double dummyChangeCost = 0.0; return changeBounds(3 + type, dummyChangeCost); } else if (type < 0) { #ifndef NDEBUG // just check int iSequence; int nFake = 0; int nErrors = 0; int nSuperBasic = 0; int nWarnings = 0; for (iSequence = 0; iSequence < numberTotal_; iSequence++) { FakeBound fakeStatus = getFakeBound(iSequence); Status status = getInternalStatus(iSequence); bool isFake = false; #ifdef CLP_INVESTIGATE char RC = 'C'; int jSequence = iSequence; if (jSequence >= numberColumns_) { RC = 'R'; jSequence -= numberColumns_; } #endif double lowerValue = lowerSaved_[iSequence]; double upperValue = upperSaved_[iSequence]; double value = abcSolution_[iSequence]; CoinRelFltEq equal; if (status == atUpperBound || status == atLowerBound) { if (fakeStatus == AbcSimplexDual::upperFake) { if (!equal(abcUpper_[iSequence], (lowerValue + currentDualBound_)) || !(equal(abcUpper_[iSequence], value) || equal(abcLower_[iSequence], value))) { nErrors++; #ifdef CLP_INVESTIGATE printf("** upperFake %c%d %g <= %g <= %g true %g, %g\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue); #endif } isFake = true; ; } else if (fakeStatus == AbcSimplexDual::lowerFake) { if (!equal(abcLower_[iSequence], (upperValue - currentDualBound_)) || !(equal(abcUpper_[iSequence], value) || equal(abcLower_[iSequence], value))) { nErrors++; #ifdef CLP_INVESTIGATE printf("** lowerFake %c%d %g <= %g <= %g true %g, %g\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue); #endif } isFake = true; ; } else if (fakeStatus == AbcSimplexDual::bothFake) { nWarnings++; #ifdef CLP_INVESTIGATE printf("** %d at bothFake?\n", iSequence); #endif } else if (abcUpper_[iSequence] - abcLower_[iSequence] > 2.0 * currentDualBound_) { nErrors++; #ifdef CLP_INVESTIGATE printf("** noFake! %c%d %g <= %g <= %g true %g, %g\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue); #endif } } else if (status == superBasic || status == isFree) { nSuperBasic++; //printf("** free or superbasic %c%d %g <= %g <= %g true %g, %g - status %d\n", // RC,jSequence,abcLower_[iSequence],abcSolution_[iSequence], // abcUpper_[iSequence],lowerValue,upperValue,status); } else if (status == basic) { bool odd = false; if (!equal(abcLower_[iSequence], lowerValue)) odd = true; if (!equal(abcUpper_[iSequence], upperValue)) odd = true; if (odd) { #ifdef CLP_INVESTIGATE printf("** basic %c%d %g <= %g <= %g true %g, %g\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue); #endif nWarnings++; } } else if (status == isFixed) { if (!equal(abcUpper_[iSequence], abcLower_[iSequence])) { nErrors++; #ifdef CLP_INVESTIGATE printf("** fixed! %c%d %g <= %g <= %g true %g, %g\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue); #endif } } if (isFake) { nFake++; } else { if (fakeStatus != AbcSimplexDual::noFake) { nErrors++; #ifdef CLP_INVESTIGATE printf("** bad fake status %c%d %d\n", RC, jSequence, fakeStatus); #endif } } } if (nFake != numberFake_) { #ifdef CLP_INVESTIGATE printf("nfake %d numberFake %d\n", nFake, numberFake_); #endif nErrors++; } if (nErrors || type <= -1000) { #ifdef CLP_INVESTIGATE printf("%d errors, %d warnings, %d free/superbasic, %d fake\n", nErrors, nWarnings, nSuperBasic, numberFake_); printf("currentDualBound %g\n", currentDualBound_); #endif if (type <= -1000) { iSequence = -type; iSequence -= 1000; #ifdef CLP_INVESTIGATE char RC = 'C'; int jSequence = iSequence; if (jSequence >= numberColumns_) { RC = 'R'; jSequence -= numberColumns_; } double lowerValue = lowerSaved_[iSequence]; double upperValue = upperSaved_[iSequence]; printf("*** movement>1.0e30 for %c%d %g <= %g <= %g true %g, %g - status %d\n", RC, jSequence, abcLower_[iSequence], abcSolution_[iSequence], abcUpper_[iSequence], lowerValue, upperValue, internalStatus_[iSequence]); #endif assert(nErrors); // should have been picked up } assert(!nErrors); } #endif } else { CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_); CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_); // reset using status numberFake_ = 0; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { FakeBound fakeStatus = getFakeBound(iSequence); if (fakeStatus != AbcSimplexDual::noFake) { Status status = getInternalStatus(iSequence); if (status == basic) { setFakeBound(iSequence, AbcSimplexDual::noFake); continue; } double lowerValue = abcLower_[iSequence]; double upperValue = abcUpper_[iSequence]; double value = abcSolution_[iSequence]; numberFake_++; if (fakeStatus == AbcSimplexDual::upperFake) { abcUpper_[iSequence] = lowerValue + currentDualBound_; if (status == AbcSimplex::atLowerBound) { abcSolution_[iSequence] = lowerValue; } else if (status == AbcSimplex::atUpperBound) { abcSolution_[iSequence] = abcUpper_[iSequence]; } else { abort(); } } else if (fakeStatus == AbcSimplexDual::lowerFake) { abcLower_[iSequence] = upperValue - currentDualBound_; if (status == AbcSimplex::atLowerBound) { abcSolution_[iSequence] = abcLower_[iSequence]; } else if (status == AbcSimplex::atUpperBound) { abcSolution_[iSequence] = upperValue; } else { abort(); } } else { assert(fakeStatus == AbcSimplexDual::bothFake); if (status == AbcSimplex::atLowerBound) { abcLower_[iSequence] = value; abcUpper_[iSequence] = value + currentDualBound_; } else if (status == AbcSimplex::atUpperBound) { abcUpper_[iSequence] = value; abcLower_[iSequence] = value - currentDualBound_; } else if (status == AbcSimplex::isFree || status == AbcSimplex::superBasic) { abcLower_[iSequence] = value - 0.5 * currentDualBound_; abcUpper_[iSequence] = value + 0.5 * currentDualBound_; } else { abort(); } } } } } return -1; } // maybe could JUST check against perturbationArray void AbcSimplex::checkDualSolutionPlusFake() { sumDualInfeasibilities_ = 0.0; numberDualInfeasibilities_ = 0; sumFakeInfeasibilities_ = 0.0; int firstFreePrimal = -1; int firstFreeDual = -1; int numberSuperBasicWithDj = 0; bestPossibleImprovement_ = 0.0; // we can't really trust infeasibilities if there is dual error double error = CoinMin(1.0e-2, largestDualError_); // allow tolerance at least slightly bigger than standard double relaxedTolerance = currentDualTolerance_ + error; // allow bigger tolerance for possible improvement double possTolerance = 5.0 * relaxedTolerance; sumOfRelaxedDualInfeasibilities_ = 0.0; //rawObjectiveValue_ -=sumNonBasicCosts_; sumNonBasicCosts_ = 0.0; const double *COIN_RESTRICT abcLower = abcLower_; const double *COIN_RESTRICT abcUpper = abcUpper_; const double *COIN_RESTRICT abcSolution = abcSolution_; #ifdef CLEANUP_DJS double *COIN_RESTRICT abcDj = abcDj_; double *COIN_RESTRICT abcCost = abcCost_; // see if we can clean up djs // may have perturbation const double *COIN_RESTRICT originalCost = abcPerturbation_; #else const double *COIN_RESTRICT abcDj = abcDj_; const double *COIN_RESTRICT abcCost = abcCost_; #endif //const double * COIN_RESTRICT abcPerturbation = abcPerturbation_; int badDjSequence = -1; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { if (getInternalStatus(iSequence) != basic) { sumNonBasicCosts_ += abcSolution[iSequence] * abcCost[iSequence]; // not basic double distanceUp = abcUpper[iSequence] - abcSolution[iSequence]; double distanceDown = abcSolution[iSequence] - abcLower[iSequence]; double value = abcDj[iSequence]; if (distanceUp > primalTolerance_) { // Check if "free" if (distanceDown <= primalTolerance_) { // should not be negative if (value < 0.0) { value = -value; if (value > currentDualTolerance_) { bool treatAsSuperBasic = false; if (getInternalStatus(iSequence) == superBasic) { #ifdef PAN if (fakeSuperBasic(iSequence) >= 0) { #endif treatAsSuperBasic = true; numberSuperBasicWithDj++; if (firstFreeDual < 0) firstFreeDual = iSequence; #ifdef PAN } #endif } if (!treatAsSuperBasic) { sumDualInfeasibilities_ += value - currentDualTolerance_; badDjSequence = iSequence; if (value > possTolerance) bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value; if (value > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; numberDualInfeasibilities_++; } } #ifdef CLEANUP_DJS } else { // adjust cost etc if (abcCost[iSequence] > originalCost[iSequence]) { double difference = CoinMin(abcCost[iSequence] - originalCost[iSequence], value); //printf("Lb %d changing cost from %g to %g, dj from %g to %g\n", // iSequence,abcCost[iSequence],abcCost[iSequence]-difference, // abcDj[iSequence],abcDj[iSequence]-difference); abcDj[iSequence] = value - difference; abcCost[iSequence] -= difference; } #endif } } else { // free value = fabs(value); if (value > 1.0 * relaxedTolerance) { numberSuperBasicWithDj++; if (firstFreeDual < 0) firstFreeDual = iSequence; } if (value > currentDualTolerance_) { sumDualInfeasibilities_ += value - currentDualTolerance_; if (value > possTolerance) bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value; if (value > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; numberDualInfeasibilities_++; } if (firstFreePrimal < 0) firstFreePrimal = iSequence; } } else if (distanceDown > primalTolerance_) { // should not be positive if (value > 0.0) { if (value > currentDualTolerance_) { bool treatAsSuperBasic = false; if (getInternalStatus(iSequence) == superBasic) { #ifdef PAN if (fakeSuperBasic(iSequence) >= 0) { #endif treatAsSuperBasic = true; numberSuperBasicWithDj++; if (firstFreeDual < 0) firstFreeDual = iSequence; #ifdef PAN } #endif } if (!treatAsSuperBasic) { sumDualInfeasibilities_ += value - currentDualTolerance_; badDjSequence = iSequence; if (value > possTolerance) bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10); if (value > relaxedTolerance) sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; numberDualInfeasibilities_++; } } #ifdef CLEANUP_DJS } else { // adjust cost etc if (abcCost[iSequence] < originalCost[iSequence]) { double difference = CoinMax(abcCost[iSequence] - originalCost[iSequence], value); //printf("Ub %d changing cost from %g to %g, dj from %g to %g\n", // iSequence,abcCost[iSequence],abcCost[iSequence]-difference, // abcDj[iSequence],abcDj[iSequence]-difference); abcDj[iSequence] = value - difference; abcCost[iSequence] -= difference; } #endif } } } } #ifdef CLEANUP_DJS if (algorithm_ < 0 && numberDualInfeasibilities_ == 1 && firstFreeDual < 0 && sumDualInfeasibilities_ < 1.0e-1) { // clean bad one assert(badDjSequence >= 0); abcCost[badDjSequence] -= abcDj_[badDjSequence]; abcDj_[badDjSequence] = 0.0; } #endif numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj; numberFreeNonBasic_ = numberSuperBasicWithDj; firstFree_ = -1; if (algorithm_ < 0 && firstFreeDual >= 0) { // dual firstFree_ = firstFreeDual; } else if ((numberSuperBasicWithDj || (abcProgress_.lastIterationNumber(0) <= 0)) && algorithm_ > 0) { firstFree_ = firstFreePrimal; } } // Perturbs problem void AbcSimplexDual::perturb(double factor) { const double *COIN_RESTRICT lower = lowerSaved_; const double *COIN_RESTRICT upper = upperSaved_; const double *COIN_RESTRICT perturbation = abcPerturbation_; double *COIN_RESTRICT cost = abcCost_; for (int i = 0; i < numberTotal_; i++) { if (getInternalStatus(i) != basic) { if (fabs(lower[i]) < fabs(upper[i])) { cost[i] += factor * perturbation[i]; } else if (upper[i] != COIN_DBL_MAX && upper[i] != lower[i]) { cost[i] -= factor * perturbation[i]; } else { // free //cost[i]=costSaved[i]; } } } #ifdef CLEANUP_DJS // copy cost to perturbation CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_); #endif } #if ABC_NORMAL_DEBUG > 0 #define PERT_STATISTICS #endif #ifdef PERT_STATISTICS static void breakdown(const char *name, int numberLook, const double *region) { double range[] = { -COIN_DBL_MAX, -1.0e15, -1.0e11, -1.0e8, -1.0e5, -3.0e4, -1.0e4, -3.0e3, -1.0e3, -3.0e2, -1.0e2, -3.0e1, -1.0e1, -3.0, -1.0, -3.0e-1, -1.0e-1, -3.0e-1, -1.0e-2, -3.0e-2, -1.0e-3, -3.0e-3, -1.0e-4, -3.0e-4, -1.0e-5, -1.0e-8, -1.0e-11, -1.0e-15, 0.0, 1.0e-15, 1.0e-11, 1.0e-8, 1.0e-5, 3.0e-5, 1.0e-4, 3.0e-4, 1.0e-3, 3.0e-3, 1.0e-2, 3.0e-2, 1.0e-1, 3.0e-1, 1.0, 3.0, 1.0e1, 3.0e1, 1.0e2, 3.0e2, 1.0e3, 3.0e3, 1.0e4, 3.0e4, 1.0e5, 1.0e8, 1.0e11, 1.0e15, COIN_DBL_MAX }; int nRanges = static_cast< int >(sizeof(range) / sizeof(double)); int *number = new int[nRanges]; memset(number, 0, nRanges * sizeof(int)); int *numberExact = new int[nRanges]; memset(numberExact, 0, nRanges * sizeof(int)); int i; for (i = 0; i < numberLook; i++) { double value = region[i]; for (int j = 0; j < nRanges; j++) { if (value == range[j]) { numberExact[j]++; break; } else if (value < range[j]) { number[j]++; break; } } } printf("\n%s has %d entries\n", name, numberLook); for (i = 0; i < nRanges; i++) { if (number[i]) printf("%d between %g and %g", number[i], range[i - 1], range[i]); if (numberExact[i]) { if (number[i]) printf(", "); printf("%d exactly at %g", numberExact[i], range[i]); } if (number[i] + numberExact[i]) printf("\n"); } delete[] number; delete[] numberExact; } #endif // Perturbs problem void AbcSimplexDual::perturbB(double /*factorIn*/, int /*type*/) { double overallMultiplier = (perturbation_ == 50 || perturbation_ > 54) ? 1.0 : 0.1; // dual perturbation double perturbation = 1.0e-20; // maximum fraction of cost to perturb double maximumFraction = 1.0e-5; int maxLength = 0; int minLength = numberRows_; double averageCost = 0.0; int numberNonZero = 0; // See if we need to perturb double *COIN_RESTRICT sort = new double[numberColumns_]; // Use objective BEFORE scaling ?? const double *COIN_RESTRICT obj = objective(); for (int i = 0; i < numberColumns_; i++) { double value = fabs(obj[i]); sort[i] = value; averageCost += fabs(abcCost_[i + maximumAbcNumberRows_]); if (value) numberNonZero++; } if (numberNonZero) averageCost /= static_cast< double >(numberNonZero); else averageCost = 1.0; // allow for cost scaling averageCost *= objectiveScale_; std::sort(sort, sort + numberColumns_); int number = 1; double last = sort[0]; for (int i = 1; i < numberColumns_; i++) { if (last != sort[i]) number++; last = sort[i]; } delete[] sort; #ifdef PERT_STATISTICS printf("%d different values\n", number); #endif #if 0 printf("nnz %d percent %d", number, (number * 100) / numberColumns_); if (number * 4 > numberColumns_) printf(" - Would not perturb\n"); else printf(" - Would perturb\n"); #endif // If small costs then more dangerous (better to scale) double numberMultiplier = 1.25; if (averageCost < 1.0e-3) numberMultiplier = 20.0; else if (averageCost < 1.0e-2) numberMultiplier = 10.0; else if (averageCost < 1.0e-1) numberMultiplier = 3.0; else if (averageCost < 1.0) numberMultiplier = 2.0; numberMultiplier *= 2.0; // try if (number * numberMultiplier > numberColumns_ || perturbation_ >= 100) { perturbation_ = 101; #if ABC_NORMAL_DEBUG > 0 printf("Not perturbing - %d different costs, average %g\n", number, averageCost); #endif #ifdef CLEANUP_DJS // copy cost to perturbation CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_); #endif return; // good enough } const int *COIN_RESTRICT columnLength = matrix_->getVectorLengths() - maximumAbcNumberRows_; const double *COIN_RESTRICT lower = lowerSaved_; const double *COIN_RESTRICT upper = upperSaved_; for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) { if (lower[iSequence] < upper[iSequence]) { int length = columnLength[iSequence]; if (length > 2) { maxLength = CoinMax(maxLength, length); minLength = CoinMin(minLength, length); } } } bool uniformChange = false; bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0; double constantPerturbation = 10.0 * dualTolerance_; #ifdef PERT_STATISTICS breakdown("Objective before perturbation", numberColumns_, abcCost_ + numberRows_); #endif //#define PERTURB_OLD_WAY #ifndef PERTURB_OLD_WAY #if 1 #ifdef HEAVY_PERTURBATION if (perturbation_ == 50) perturbation_ = HEAVY_PERTURBATION; #endif if (perturbation_ >= 50) { // Experiment // maximumFraction could be 1.0e-10 to 1.0 double m[] = { 1.0e-10, 1.0e-9, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0 }; int whichOne = CoinMin(perturbation_ - 51, 10); if (whichOne < 0) { if (numberRows_ < 5000) whichOne = 2; // treat 50 as if 53 else whichOne = 2; // treat 50 as if 52 } if (inCbcOrOther && whichOne > 0) whichOne--; maximumFraction = m[whichOne]; if (whichOne > 7) constantPerturbation *= 10.0; //if (whichOne<2) //constantPerturbation = maximumFraction*1.0e-10 * dualTolerance_; } else if (inCbcOrOther) { maximumFraction = 1.0e-6; } #endif double smallestNonZero = 1.0e100; bool modifyRowCosts = false; numberNonZero = 0; //perturbation = 1.0e-10; perturbation = CoinMax(1.0e-10, maximumFraction); double *COIN_RESTRICT cost = abcCost_; bool allSame = true; double lastValue = 0.0; for (int iRow = 0; iRow < numberRows_; iRow++) { double lo = lower[iRow]; double up = upper[iRow]; if (lo < up) { double value = fabs(cost[iRow]); perturbation = CoinMax(perturbation, value); if (value) { modifyRowCosts = true; smallestNonZero = CoinMin(smallestNonZero, value); } } if (lo && lo > -1.0e10) { numberNonZero++; lo = fabs(lo); if (!lastValue) lastValue = lo; else if (fabs(lo - lastValue) > 1.0e-7) allSame = false; } if (up && up < 1.0e10) { numberNonZero++; up = fabs(up); if (!lastValue) lastValue = up; else if (fabs(up - lastValue) > 1.0e-7) allSame = false; } } double lastValue2 = 0.0; for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) { double lo = lower[iSequence]; double up = upper[iSequence]; if (lo < up) { double value = fabs(cost[iSequence]); perturbation = CoinMax(perturbation, value); if (value) { smallestNonZero = CoinMin(smallestNonZero, value); } } if (lo && lo > -1.0e10) { //numberNonZero++; lo = fabs(lo); if (!lastValue2) lastValue2 = lo; else if (fabs(lo - lastValue2) > 1.0e-7) allSame = false; } if (up && up < 1.0e10) { //numberNonZero++; up = fabs(up); if (!lastValue2) lastValue2 = up; else if (fabs(up - lastValue2) > 1.0e-7) allSame = false; } } if (allSame) { // Check elements double smallestNegative; double largestNegative; double smallestPositive; double largestPositive; matrix_->rangeOfElements(smallestNegative, largestNegative, smallestPositive, largestPositive); if (smallestNegative == largestNegative && smallestPositive == largestPositive) { // Really hit perturbation double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2)); maximumFraction = CoinMax(adjust, maximumFraction); } } perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction); double largestZero = 0.0; double largest = 0.0; double largestPerCent = 0.0; // modify costs bool printOut = (handler_->logLevel() == 63); printOut = false; //assert (!modifyRowCosts); modifyRowCosts = false; const double *COIN_RESTRICT perturbationArray = perturbationSaved_; if (modifyRowCosts) { for (int iRow = 0; iRow < numberRows_; iRow++) { if (lower[iRow] < upper[iRow]) { double value = perturbation; double currentValue = cost[iRow]; value = CoinMin(value, maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-3)); double perturbationValue = 2.0 * (perturbationArray[iRow] - 0.5); // for now back to normal random if (lower[iRow] > -largeValue_) { if (fabs(lower[iRow]) < fabs(upper[iRow])) value *= perturbationValue; else value *= -perturbationValue; } else if (upper[iRow] < largeValue_) { value *= -perturbationValue; } else { value = 0.0; } if (currentValue) { largest = CoinMax(largest, fabs(value)); if (fabs(value) > fabs(currentValue) * largestPerCent) largestPerCent = fabs(value / currentValue); } else { largestZero = CoinMax(largestZero, fabs(value)); } if (printOut) printf("row %d cost %g change %g\n", iRow, cost[iRow], value); cost[iRow] += value; } } } //double extraWeight = 10.0; // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0}; // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0}; double weight[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0 }; //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0}; // Scale back if wanted double weight2[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 }; if (constantPerturbation < 99.0 * dualTolerance_ && false) { perturbation *= 0.1; //extraWeight = 0.5; memcpy(weight, weight2, sizeof(weight2)); } // adjust weights if all columns long double factor = 1.0; if (maxLength) { factor = 3.0 / static_cast< double >(minLength); } // Make variables with more elements more expensive const double m1 = 0.5; double smallestAllowed = overallMultiplier * CoinMin(1.0e-2 * dualTolerance_, maximumFraction); double largestAllowed = overallMultiplier * CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost); // smaller if in BAB //if (inCbcOrOther) //largestAllowed=CoinMin(largestAllowed,1.0e-5); //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed); randomNumberGenerator_.setSeed(1000000000); #if ABC_NORMAL_DEBUG > 0 printf("perturbation %g constantPerturbation %g smallestNonZero %g maximumFraction %g smallestAllowed %g\n", perturbation, constantPerturbation, smallestNonZero, maximumFraction, smallestAllowed); #endif for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) { if (lower[iSequence] < upper[iSequence] && getInternalStatus(iSequence) != basic) { double value = perturbation; double currentValue = cost[iSequence]; double currentLargestAllowed = fabs(currentValue) * 1.0e-3; if (!currentValue) currentLargestAllowed = 1.0e-3; if (currentLargestAllowed <= 10.0 * smallestAllowed) continue; // don't perturb tiny value value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8)); //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue)); double value2 = constantPerturbation + 1.0e-1 * smallestNonZero; if (uniformChange) { value = maximumFraction; value2 = maximumFraction; } //double perturbationValue=2.0*(perturbationArray[iSequence]-0.5); // for now back to normal random double perturbationValue = randomNumberGenerator_.randomDouble(); perturbationValue = 1.0 - m1 + m1 * perturbationValue; // clean up if (lower[iSequence] > -largeValue_) { if (fabs(lower[iSequence]) < fabs(upper[iSequence])) { value *= perturbationValue; value2 *= perturbationValue; } else { //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble()); //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble()); value = 0.0; } } else if (upper[iSequence] < largeValue_) { value *= -perturbationValue; value2 *= -perturbationValue; } else { value = 0.0; } if (value) { int length = columnLength[iSequence]; if (length > 3) { length = static_cast< int >(static_cast< double >(length) * factor); length = CoinMax(3, length); } double multiplier; if (length < 10) multiplier = weight[length]; else multiplier = weight[10]; value *= multiplier; value = CoinMin(value, value2); if (value) { // get in range if (fabs(value) <= smallestAllowed) { value *= 10.0; while (fabs(value) <= smallestAllowed) value *= 10.0; } else if (fabs(value) > currentLargestAllowed) { value *= 0.1; while (fabs(value) > currentLargestAllowed) value *= 0.1; } } if (currentValue) { largest = CoinMax(largest, fabs(value)); if (fabs(value) > fabs(currentValue) * largestPerCent) largestPerCent = fabs(value / currentValue); } else { largestZero = CoinMax(largestZero, fabs(value)); } // but negative if at ub if (getInternalStatus(iSequence) == atUpperBound) value = -value; if (printOut) printf("col %d cost %g change %g\n", iSequence, cost[iSequence], value); cost[iSequence] += value; } } } #else // old way if (perturbation_ > 50) { // Experiment // maximumFraction could be 1.0e-10 to 1.0 double m[] = { 1.0e-10, 1.0e-9, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0 }; int whichOne = perturbation_ - 51; //if (inCbcOrOther&&whichOne>0) //whichOne--; maximumFraction = m[CoinMin(whichOne, 10)]; } else if (inCbcOrOther) { //maximumFraction = 1.0e-6; } int iRow; double smallestNonZero = 1.0e100; numberNonZero = 0; if (perturbation_ >= 50) { perturbation = 1.0e-8; if (perturbation_ > 50 && perturbation_ < 60) perturbation = CoinMax(1.0e-8, maximumFraction); bool allSame = true; double lastValue = 0.0; for (iRow = 0; iRow < numberRows_; iRow++) { double lo = abcLower_[iRow]; double up = abcUpper_[iRow]; if (lo && lo > -1.0e10) { numberNonZero++; lo = fabs(lo); if (!lastValue) lastValue = lo; else if (fabs(lo - lastValue) > 1.0e-7) allSame = false; } if (up && up < 1.0e10) { numberNonZero++; up = fabs(up); if (!lastValue) lastValue = up; else if (fabs(up - lastValue) > 1.0e-7) allSame = false; } } double lastValue2 = 0.0; for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) { double lo = abcLower_[iColumn]; double up = abcUpper_[iColumn]; if (lo < up) { double value = fabs(abcCost_[iColumn]); perturbation = CoinMax(perturbation, value); if (value) { smallestNonZero = CoinMin(smallestNonZero, value); } } if (lo && lo > -1.0e10) { //numberNonZero++; lo = fabs(lo); if (!lastValue2) lastValue2 = lo; else if (fabs(lo - lastValue2) > 1.0e-7) allSame = false; } if (up && up < 1.0e10) { //numberNonZero++; up = fabs(up); if (!lastValue2) lastValue2 = up; else if (fabs(up - lastValue2) > 1.0e-7) allSame = false; } } if (allSame) { // Check elements double smallestNegative; double largestNegative; double smallestPositive; double largestPositive; matrix_->rangeOfElements(smallestNegative, largestNegative, smallestPositive, largestPositive); if (smallestNegative == largestNegative && smallestPositive == largestPositive) { // Really hit perturbation double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2)); maximumFraction = CoinMax(adjust, maximumFraction); } } perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction); } double largestZero = 0.0; double largest = 0.0; double largestPerCent = 0.0; // modify costs bool printOut = (handler_->logLevel() == 63); printOut = false; // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0}; // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0}; double weight[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0 }; //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0}; //double extraWeight = 10.0; // Scale back if wanted double weight2[] = { 1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 }; if (constantPerturbation < 99.0 * dualTolerance_) { perturbation *= 0.1; //extraWeight = 0.5; memcpy(weight, weight2, sizeof(weight2)); } // adjust weights if all columns long double factor = 1.0; if (maxLength) { factor = 3.0 / static_cast< double >(minLength); } // Make variables with more elements more expensive const double m1 = 0.5; double smallestAllowed = CoinMin(1.0e-2 * dualTolerance_, maximumFraction); double largestAllowed = CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost); // smaller if in BAB //if (inCbcOrOther) //largestAllowed=CoinMin(largestAllowed,1.0e-5); //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed); for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) { if (abcLower_[iColumn] < abcUpper_[iColumn] && getInternalStatus(iColumn) != basic) { double value = perturbation; double currentValue = abcCost_[iColumn]; value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8)); //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue)); double value2 = constantPerturbation + 1.0e-1 * smallestNonZero; if (uniformChange) { value = maximumFraction; value2 = maximumFraction; } if (abcLower_[iColumn] > -largeValue_) { if (fabs(abcLower_[iColumn]) < fabs(abcUpper_[iColumn])) { value *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble()); value2 *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble()); } else { //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble()); //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble()); value = 0.0; } } else if (abcUpper_[iColumn] < largeValue_) { value *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble()); value2 *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble()); } else { value = 0.0; } if (value) { int length = columnLength[iColumn]; if (length > 3) { length = static_cast< int >(static_cast< double >(length) * factor); length = CoinMax(3, length); } double multiplier; #if 1 if (length < 10) multiplier = weight[length]; else multiplier = weight[10]; #else if (length < 10) multiplier = weight[length]; else multiplier = weight[10] + extraWeight * (length - 10); multiplier *= 0.5; #endif value *= multiplier; value = CoinMin(value, value2); if (value) { // get in range if (fabs(value) <= smallestAllowed) { value *= 10.0; while (fabs(value) <= smallestAllowed) value *= 10.0; } else if (fabs(value) > largestAllowed) { value *= 0.1; while (fabs(value) > largestAllowed) value *= 0.1; } } if (currentValue) { largest = CoinMax(largest, fabs(value)); if (fabs(value) > fabs(currentValue) * largestPerCent) largestPerCent = fabs(value / currentValue); } else { largestZero = CoinMax(largestZero, fabs(value)); } // but negative if at ub if (getInternalStatus(iColumn) == atUpperBound) value = -value; if (printOut) printf("col %d cost %g change %g\n", iColumn, objectiveWork_[iColumn], value); abcCost_[iColumn] += value; } } } #endif #ifdef PERT_STATISTICS { double averageCost = 0.0; int numberNonZero = 0; double *COIN_RESTRICT sort = new double[numberColumns_]; for (int i = 0; i < numberColumns_; i++) { double value = fabs(abcCost_[i + maximumAbcNumberRows_]); sort[i] = value; averageCost += value; if (value) numberNonZero++; } if (numberNonZero) averageCost /= static_cast< double >(numberNonZero); else averageCost = 1.0; std::sort(sort, sort + numberColumns_); int number = 1; double last = sort[0]; for (int i = 1; i < numberColumns_; i++) { if (last != sort[i]) number++; last = sort[i]; } delete[] sort; #if ABC_NORMAL_DEBUG > 0 printf("nnz %d percent %d", number, (number * 100) / numberColumns_); breakdown("Objective after perturbation", numberColumns_, abcCost_ + numberRows_); #endif } #endif perturbation_ = 100; // copy cost to perturbation #ifdef CLEANUP_DJS CoinAbcMemcpy(abcPerturbation_, abcCost_, numberTotal_); #endif handler_->message(CLP_SIMPLEX_PERTURB, messages_) << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero << CoinMessageEol; } /* Does something about fake tolerances -1 initializes 1 might be optimal (back to original costs) 2 after change currentDualBound 3 might be optimal BUT change currentDualBound */ int AbcSimplexDual::bounceTolerances(int type) { // for now turn off if (type == -1) { if (!perturbationFactor_) perturbationFactor_ = 0.001; CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); if (!perturbationFactor_) { //currentDualBound_=CoinMin(dualBound_,1.0e11); //CoinAbcMultiplyAdd(perturbationSaved_,numberTotal_,currentDualTolerance_, // abcPerturbation_,0.0); } else { //currentDualBound_=CoinMin(dualBound_,1.0e11); double factor = currentDualTolerance_ * perturbationFactor_; perturbB(factor, 0); } #if ABC_NORMAL_DEBUG > 0 //printf("new dual bound %g (bounceTolerances)\n",currentDualBound_); #endif } else { if (stateDualColumn_ == -2) { double newTolerance = CoinMin(1.0e-1, dualTolerance_ * 100.0); #if ABC_NORMAL_DEBUG > 0 printf("After %d iterations reducing dual tolerance from %g to %g\n", numberIterations_, currentDualTolerance_, newTolerance); #endif currentDualTolerance_ = newTolerance; stateDualColumn_ = -1; } else if (stateDualColumn_ == -1) { #if ABC_NORMAL_DEBUG > 0 printf("changing dual tolerance from %g to %g (vaguely optimal)\n", currentDualTolerance_, dualTolerance_); #endif currentDualTolerance_ = dualTolerance_; stateDualColumn_ = 0; } else { if (stateDualColumn_ > 1) { normalDualColumnIteration_ = numberIterations_ + CoinMax(10, numberRows_ / 4); #if ABC_NORMAL_DEBUG > 0 printf("no cost modification until iteration %d\n", normalDualColumnIteration_); #endif } } resetFakeBounds(0); CoinAbcMemcpy(abcCost_, abcPerturbation_, numberTotal_); numberChanged_ = 0; // Number of variables with changed costs moveToBasic(); if (type != 1) return 0; // save status and solution CoinAbcMemcpy(status_, internalStatus_, numberTotal_); CoinAbcMemcpy(solutionSaved_, abcSolution_, numberTotal_); #if ABC_NORMAL_DEBUG > 0 printf("A sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif gutsOfSolution(3); #if ABC_NORMAL_DEBUG > 0 printf("B sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif double saveSum1 = sumPrimalInfeasibilities_; double saveSumD1 = sumDualInfeasibilities_; bool goToPrimal = false; if (!sumOfRelaxedDualInfeasibilities_) { if (perturbation_ < 101) { // try going back to original costs CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); perturbation_ = 101; numberChanged_ = 0; // Number of variables with changed costs moveToBasic(); gutsOfSolution(3); #if ABC_NORMAL_DEBUG > 0 printf("B2 sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif saveSumD1 = sumDualInfeasibilities_; } if (sumOfRelaxedDualInfeasibilities_) { makeNonFreeVariablesDualFeasible(); moveToBasic(); abcProgress_.reset(); gutsOfSolution(3); #if ABC_NORMAL_DEBUG > 0 printf("C sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif if (sumPrimalInfeasibilities_ > 1.0e2 * saveSum1 + currentDualBound_ * 0.2 /*|| sumOfRelaxedDualInfeasibilities_<1.0e-4*sumPrimalInfeasibilities_*/) { numberDisasters_++; if (numberDisasters_ > 2) goToPrimal = true; } } else { numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; } } else { // carry on makeNonFreeVariablesDualFeasible(); moveToBasic(); abcProgress_.reset(); gutsOfSolution(3); #if ABC_NORMAL_DEBUG > 0 printf("C2 sumDual %g sumPrimal %g\n", sumDualInfeasibilities_, sumPrimalInfeasibilities_); #endif numberDisasters_++; if (sumPrimalInfeasibilities_ > 1.0e2 * saveSum1 + currentDualBound_ * 0.2 || (saveSumD1 < 1.0e-4 * sumPrimalInfeasibilities_ && saveSumD1 < 1.0)) { numberDisasters_++; if (numberDisasters_ > 2 || saveSumD1 < 1.0e-10 * sumPrimalInfeasibilities_) goToPrimal = true; } } if (goToPrimal) { // go to primal CoinAbcMemcpy(internalStatus_, status_, numberTotal_); CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_); int *COIN_RESTRICT abcPivotVariable = abcPivotVariable_; // redo pivotVariable int numberBasic = 0; for (int i = 0; i < numberTotal_; i++) { if (getInternalStatus(i) == basic) abcPivotVariable[numberBasic++] = i; } assert(numberBasic == numberRows_); moveToBasic(); checkPrimalSolution(false); #if ABC_NORMAL_DEBUG > 0 printf("D sumPrimal %g\n", sumPrimalInfeasibilities_); #endif return 1; } } CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_); return 0; } #if ABC_NORMAL_DEBUG > 0 extern int cilk_conflict; #endif int AbcSimplexDual::noPivotRow() { #if ABC_NORMAL_DEBUG > 3 if (handler_->logLevel() & 32) printf("** no row pivot\n"); #endif int numberPivots = abcFactorization_->pivots(); bool specialCase; int useNumberFake; int returnCode = 0; if (numberPivots <= CoinMax(dontFactorizePivots_, 20) && (specialOptions_ & 2048) != 0 && currentDualBound_ >= 1.0e8) { specialCase = true; // as dual bound high - should be okay useNumberFake = 0; } else { specialCase = false; useNumberFake = numberAtFakeBound(); } if (!numberPivots || specialCase) { // may have crept through - so may be optimal // check any flagged variables if (numberFlagged_ && numberPivots) { // try factorization returnCode = -2; } if (useNumberFake || numberDualInfeasibilities_) { // may be dual infeasible if ((specialOptions_ & 1024) == 0) problemStatus_ = -5; else if (!useNumberFake && numberPrimalInfeasibilities_ && !numberPivots) problemStatus_ = 1; } else { if (numberFlagged_) { #if ABC_NORMAL_DEBUG > 3 std::cout << "Flagged variables at end - infeasible?" << std::endl; printf("Probably infeasible - pivot was %g\n", alpha_); #endif #if ABC_NORMAL_DEBUG > 3 abort(); #endif problemStatus_ = -5; } else { problemStatus_ = 0; #ifndef CLP_CHECK_NUMBER_PIVOTS #define CLP_CHECK_NUMBER_PIVOTS 10 #endif #if CLP_CHECK_NUMBER_PIVOTS < 20 if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) { gutsOfSolution(2); rawObjectiveValue_ += sumNonBasicCosts_; setClpSimplexObjectiveValue(); if (numberPrimalInfeasibilities_) { problemStatus_ = -1; returnCode = -2; } } #endif if (!problemStatus_) { // make it look OK numberPrimalInfeasibilities_ = 0; sumPrimalInfeasibilities_ = 0.0; numberDualInfeasibilities_ = 0; sumDualInfeasibilities_ = 0.0; if (numberChanged_) { numberChanged_ = 0; // Number of variables with changed costs CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_); CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_); // make sure duals are current gutsOfSolution(1); setClpSimplexObjectiveValue(); abcProgress_.modifyObjective(-COIN_DBL_MAX); if (numberDualInfeasibilities_) { problemStatus_ = 10; // was -3; numberTimesOptimal_++; } else { computeObjectiveValue(true); } } else if (numberPivots) { computeObjectiveValue(true); } sumPrimalInfeasibilities_ = 0.0; } if ((specialOptions_ & (1024 + 16384)) != 0 && !problemStatus_ && abcFactorization_->pivots()) { CoinIndexedVector *arrayVector = &usefulArray_[arrayForFtran_]; double *rhs = arrayVector->denseVector(); arrayVector->clear(); CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_); abcMatrix_->timesModifyExcludingSlacks(-1.0, abcSolution_, rhs); //#define CHECK_ACCURACY #ifdef CHECK_ACCURACY bool bad = false; #endif bool bad2 = false; int i; for (i = 0; i < numberRows_; i++) { if (rhs[i] < abcLower_[i] - primalTolerance_ || rhs[i] > abcUpper_[i] + primalTolerance_) { bad2 = true; #ifdef CHECK_ACCURACY printf("row %d out of bounds %g, %g correct %g bad %g\n", i, abcLower_[i], abcUpper_[i], rhs[i], abcSolution_[i]); #endif } else if (fabs(rhs[i] - abcSolution_[i]) > 1.0e-3) { #ifdef CHECK_ACCURACY bad = true; printf("row %d correct %g bad %g\n", i, rhs[i], abcSolution_[i]); #endif } rhs[i] = 0.0; } for (int i = 0; i < numberTotal_; i++) { if (abcSolution_[i] < abcLower_[i] - primalTolerance_ || abcSolution_[i] > abcUpper_[i] + primalTolerance_) { bad2 = true; #ifdef CHECK_ACCURACY printf("column %d out of bounds %g, %g value %g\n", i, abcLower_[i], abcUpper_[i], abcSolution_[i]); #endif } } if (bad2) { problemStatus_ = -3; returnCode = -2; // Force to re-factorize early next time int numberPivots = abcFactorization_->pivots(); if (forceFactorization_ < 0) forceFactorization_ = 100000; forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); } } } } } else { problemStatus_ = -3; returnCode = -2; // Force to re-factorize early next time int numberPivots = abcFactorization_->pivots(); if (forceFactorization_ < 0) forceFactorization_ = 100000; forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); } return returnCode; } void AbcSimplexDual::dualPivotColumn() { double saveAcceptable = acceptablePivot_; if (largestPrimalError_ > 1.0e-5 || largestDualError_ > 1.0e-5) { //if ((largestPrimalError_>1.0e-5||largestDualError_>1.0e-5)&&false) { if (!abcFactorization_->pivots()) acceptablePivot_ *= 1.0e2; else if (abcFactorization_->pivots() < 5) acceptablePivot_ *= 1.0e1; } dualColumn2(); if (sequenceIn_ < 0 && acceptablePivot_ > saveAcceptable) { acceptablePivot_ = saveAcceptable; #if ABC_NORMAL_DEBUG > 0 printf("No pivot!!!\n"); #endif // try again dualColumn1(true); dualColumn2(); } else { acceptablePivot_ = saveAcceptable; } if ((stateOfProblem_ & VALUES_PASS) != 0) { // check no flips assert(!usefulArray_[arrayForFlipBounds_].getNumElements()); // If no pivot - then try other way if plausible if (sequenceIn_ < 0) { // could use fake primal basic values if ((directionOut_ < 0 && lowerOut_ > -1.0e30) || (directionOut_ > 0 && upperOut_ < 1.30)) { directionOut_ = -directionOut_; if (directionOut_ < 0 && abcDj_[sequenceOut_] < 0.0) upperTheta_ = -abcDj_[sequenceOut_]; else if (directionOut_ > 0 && abcDj_[sequenceOut_] > 0.0) upperTheta_ = abcDj_[sequenceOut_]; else upperTheta_ = 1.0e30; CoinPartitionedVector &tableauRow = usefulArray_[arrayForTableauRow_]; CoinPartitionedVector &candidateList = usefulArray_[arrayForDualColumn_]; assert(!candidateList.getNumElements()); tableauRow.compact(); candidateList.compact(); // redo candidate list int *COIN_RESTRICT index = tableauRow.getIndices(); double *COIN_RESTRICT array = tableauRow.denseVector(); double *COIN_RESTRICT arrayCandidate = candidateList.denseVector(); int *COIN_RESTRICT indexCandidate = candidateList.getIndices(); const double *COIN_RESTRICT abcDj = abcDj_; const unsigned char *COIN_RESTRICT internalStatus = internalStatus_; // do first pass to get possibles double bestPossible = 0.0; // We can also see if infeasible or pivoting on free double tentativeTheta = 1.0e25; // try with this much smaller as guess double acceptablePivot = currentAcceptablePivot_; double dualT = -currentDualTolerance_; // fixed will have been taken out by now const double multiplier[] = { 1.0, -1.0 }; int freeSequence = -1; double freeAlpha = 0.0; int numberNonZero = tableauRow.getNumElements(); int numberRemaining = 0; if (ordinaryVariables()) { for (int i = 0; i < numberNonZero; i++) { int iSequence = index[i]; double tableauValue = array[i]; unsigned char iStatus = internalStatus[iSequence] & 7; double mult = multiplier[iStatus]; double alpha = tableauValue * mult; double oldValue = abcDj[iSequence] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta_ * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta_ = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iSequence; } } } else { double badFree = 0.0; freeAlpha = currentAcceptablePivot_; double currentDualTolerance = currentDualTolerance_; for (int i = 0; i < numberNonZero; i++) { int iSequence = index[i]; double tableauValue = array[i]; unsigned char iStatus = internalStatus[iSequence] & 7; if ((iStatus & 4) == 0) { double mult = multiplier[iStatus]; double alpha = tableauValue * mult; double oldValue = abcDj[iSequence] * mult; double value = oldValue - tentativeTheta * alpha; if (value < dualT) { bestPossible = CoinMax(bestPossible, alpha); value = oldValue - upperTheta_ * alpha; if (value < dualT && alpha >= acceptablePivot) { upperTheta_ = (oldValue - dualT) / alpha; } // add to list arrayCandidate[numberRemaining] = alpha; indexCandidate[numberRemaining++] = iSequence; } } else { bool keep; bestPossible = CoinMax(bestPossible, fabs(tableauValue)); double oldValue = abcDj[iSequence]; // If free has to be very large - should come in via dualRow //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3) //break; if (oldValue > currentDualTolerance) { keep = true; } else if (oldValue < -currentDualTolerance) { keep = true; } else { if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) { keep = true; } else { keep = false; badFree = CoinMax(badFree, fabs(tableauValue)); } } if (keep) { // free - choose largest if (fabs(tableauValue) > fabs(freeAlpha)) { freeAlpha = tableauValue; freeSequence = iSequence; } } } } } candidateList.setNumElements(numberRemaining); if (freeSequence >= 0) { sequenceIn_ = freeSequence; } else if (fabs(upperTheta_ - fabs(abcDj_[sequenceOut_])) < dualTolerance_) { // can just move stateOfIteration_ = -1; return; } dualColumn2(); } } // redo dualOut_ if (directionOut_ < 0) { dualOut_ = valueOut_ - upperOut_; } else { dualOut_ = lowerOut_ - valueOut_; } } usefulArray_[arrayForDualColumn_].clear(); //clearArrays(arrayForDualColumn_); if (sequenceIn_ >= 0 && upperTheta_ == COIN_DBL_MAX) { lowerIn_ = abcLower_[sequenceIn_]; upperIn_ = abcUpper_[sequenceIn_]; valueIn_ = abcSolution_[sequenceIn_]; dualIn_ = abcDj_[sequenceIn_]; if (valueIn_ < lowerIn_ + primalTolerance_) { directionIn_ = 1; } else if (valueIn_ > upperIn_ - primalTolerance_) { directionIn_ = -1; } else { if (alpha_ < 0.0) { // as if from upper bound directionIn_ = -1; //assert(upperIn_ == valueIn_); } else { // as if from lower bound directionIn_ = 1; //assert(lowerIn_ == valueIn_); } } } #if ABC_DEBUG checkArrays(4); #endif #if CAN_HAVE_ZERO_OBJ > 1 if ((specialOptions_ & 2097152) != 0) theta_ = 0.0; #endif movement_ = 0.0; objectiveChange_ = 0.0; /* 0 - take iteration 1 - don't take but continue 2 - don't take and break */ btranAlpha_ = -alpha_ * directionOut_; // for check } void AbcSimplexDual::getTableauColumnPart2() { #if ABC_PARALLEL if (parallelMode_ != 0) { abcFactorization_->updateColumnFTPart2(usefulArray_[arrayForFtran_]); // pivot element alpha_ = usefulArray_[arrayForFtran_].denseVector()[pivotRow_]; } #endif } void AbcSimplexDual::replaceColumnPart3() { abcFactorization_->replaceColumnPart3(this, &usefulArray_[arrayForReplaceColumn_], &usefulArray_[arrayForFtran_], lastPivotRow_, ftAlpha_ ? ftAlpha_ : alpha_); } void AbcSimplexDual::checkReplacePart1() { //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_); //usefulArray_[arrayForReplaceColumn_].print(); ftAlpha_ = abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_], pivotRow_); } #if 0 void AbcSimplexDual::checkReplacePart1a() { abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_); } void AbcSimplexDual::checkReplacePart1b() { //usefulArray_[arrayForReplaceColumn_].print(); ftAlpha_=abcFactorization_->checkReplacePart1b(&usefulArray_[arrayForReplaceColumn_],pivotRow_); } #endif /* returns 0 - OK 1 - take iteration then re-factorize 2 - flag something and skip 3 - break and re-factorize 5 - take iteration then re-factorize because of memory */ int AbcSimplexDual::checkReplace() { int returnCode = 0; if (!stateOfIteration_) { int saveStatus = problemStatus_; // check update #if MOVE_REPLACE_PART1A < 0 checkReplacePart1(); #endif int updateStatus = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_); #if ABC_DEBUG checkArrays(); #endif // If looks like bad pivot - refactorize if (fabs(dualOut_) > 1.0e50) updateStatus = 2; // if no pivots, bad update but reasonable alpha - take and invert if (updateStatus == 2 && !abcFactorization_->pivots() && fabs(alpha_) > 1.0e-5) updateStatus = 4; if (updateStatus == 1 || updateStatus == 4) { // slight error if (abcFactorization_->pivots() > 5 || updateStatus == 4) { problemStatus_ = -2; // factorize now returnCode = -3; } } else if (updateStatus == 2) { // major error // later we may need to unwind more e.g. fake bounds if (abcFactorization_->pivots() && ((moreSpecialOptions_ & 16) == 0 || abcFactorization_->pivots() > 4)) { problemStatus_ = -2; // factorize now returnCode = -2; moreSpecialOptions_ |= 16; stateOfIteration_ = 2; } else { assert(sequenceOut_ >= 0); // need to reject something char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 0 printf("BAD flag b %g pivotRow %d (seqout %d) sequenceIn %d\n", alpha_, pivotRow_, sequenceOut_, sequenceIn_); #endif setFlagged(sequenceOut_); // so can't happen immediately again sequenceOut_ = -1; //abcProgress_.clearBadTimes(); lastBadIteration_ = numberIterations_; // say be more cautious //clearArrays(-1); //if(primalSolutionUpdated) //gutsOfSolution(2); stateOfIteration_ = 1; } } else if (updateStatus == 3) { // out of memory // increase space if not many iterations if (abcFactorization_->pivots() < 0.5 * abcFactorization_->maximumPivots() && abcFactorization_->pivots() < 200) abcFactorization_->areaFactor( abcFactorization_->areaFactor() * 1.1); problemStatus_ = -2; // factorize now } else if (updateStatus == 5) { problemStatus_ = -2; // factorize now } if (theta_ < 0.0) { #if ABC_NORMAL_DEBUG > 3 if (handler_->logLevel() & 32) printf("negative theta %g\n", theta_); #endif theta_ = 0.0; } //printf("check pstatus %d ustatus %d returncode %d nels %d\n", // problemStatus_,updateStatus,returnCode,usefulArray_[arrayForReplaceColumn_].getNumElements()); if (saveStatus != -1) problemStatus_ = saveStatus; } return returnCode; } int AbcSimplexDual::noPivotColumn() { // no incoming column is valid int returnCode = 0; if (abcFactorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8) acceptablePivot_ = 1.0e-8; returnCode = 1; #if ABC_NORMAL_DEBUG > 3 if (handler_->logLevel() & 32) printf("** no column pivot\n"); #endif if (abcFactorization_->pivots() <= dontFactorizePivots_ && acceptablePivot_ <= 1.0e-8 && lastFirstFree_ < 0) { // If not in branch and bound etc save ray delete[] ray_; if ((specialOptions_ & (1024 | 4096)) == 0 || (specialOptions_ & 32) != 0) { // create ray anyway ray_ = new double[numberRows_]; setUsedArray(0); usefulArray_[arrayForBtran_].createOneUnpackedElement(pivotRow_, -directionOut_); abcFactorization_->updateColumnTranspose(usefulArray_[arrayForBtran_]); CoinMemcpyN(usefulArray_[arrayForBtran_].denseVector(), numberRows_, ray_); clearArrays(0); } else { ray_ = NULL; } // If we have just factorized and infeasibility reasonable say infeas double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13; if (largestGap_) dualTest = CoinMin(largestGap_, dualTest); // but ignore if none at fake bound if (currentDualBound_ < dualTest) { if (!numberAtFakeBound()) dualTest = 1.0e-12; } if (((specialOptions_ & 4096) != 0 || fabs(alpha_) < 1.0e-11) && currentDualBound_ >= dualTest) { double testValue = 1.0e-4; if (!abcFactorization_->pivots() && numberPrimalInfeasibilities_ == 1) testValue = 1.0e-6; if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue || (specialOptions_ & 64) == 0) { // say infeasible problemStatus_ = 1; // unless primal feasible!!!! //#define TEST_CLP_NODE #ifndef TEST_CLP_NODE // Should be correct - but ... int numberFake = numberAtFakeBound(); double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_; if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 || (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && abcFactorization_->pivots())) { if (sumPrimal > 50.0 && abcFactorization_->pivots() > 2) { problemStatus_ = -4; } else { problemStatus_ = 10; // Get rid of objective if ((specialOptions_ & 16384) == 0) objective_ = new ClpLinearObjective(NULL, numberColumns_); } } #else if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) { if ((specialOptions_ & 1024) != 0 && abcFactorization_->pivots()) { problemStatus_ = 10; // Get rid of objective if ((specialOptions_ & 16384) == 0) objective_ = new ClpLinearObjective(NULL, numberColumns_); } } #endif returnCode = 1; stateOfIteration_ = 2; } } if (abcFactorization_->pivots() != 0) { // If special option set - put off as long as possible if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) { if (numberTimesOptimal_ < 3 && abcFactorization_->pivots() > dontFactorizePivots_) problemStatus_ = -4; //say looks infeasible else problemStatus_ = 1; //say looks really infeasible } else { assert(sequenceOut_ >= 0); // flag char x = isColumn(sequenceOut_) ? 'C' : 'R'; handler_->message(CLP_SIMPLEX_FLAG, messages_) << x << sequenceWithin(sequenceOut_) << CoinMessageEol; #if ABC_NORMAL_DEBUG > 3 printf("flag c\n"); #endif setFlagged(sequenceOut_); lastBadIteration_ = numberIterations_; // say be more cautious // so can't happen immediately again sequenceOut_ = -1; if (!abcFactorization_->pivots()) { stateOfIteration_ = 1; } } } } pivotRow_ = -1; lastPivotRow_ = -1; return returnCode; } #if 0 #include "CoinTime.hpp" // Timed wait in nanoseconds - if negative then seconds static int timedWait(AbcSimplexDual * dual,int time, int which) { pthread_mutex_lock(dual->mutexPointer(which)); struct timespec absTime; clock_gettime(CLOCK_REALTIME, &absTime); if (time > 0) { absTime.tv_nsec += time; if (absTime.tv_nsec >= 1000000000) { absTime.tv_nsec -= 1000000000; absTime.tv_sec++; } } else { absTime.tv_sec -= time; } int rc = 0; while (! dual->threadInfoPointer()->status && rc == 0) { printf("waiting for signal on %d\n",which); rc = pthread_cond_timedwait(dual->conditionPointer(which), dual->mutexPointer(which), &absTime); } if (rc == 0) rc=dual->threadInfoPointer()->status; printf("got signal %d on %d\n",rc,which); //pthread_mutex_unlock(dual->mutexPointer(which)); return rc; } #endif #if ABC_PARALLEL == 1 void *abc_parallelManager(void *simplex) { AbcSimplexDual *dual = reinterpret_cast< AbcSimplexDual * >(simplex); int whichThread = dual->whichThread(); CoinThreadInfo *threadInfo = dual->threadInfoPointer(whichThread); pthread_mutex_lock(dual->mutexPointer(2, whichThread)); pthread_barrier_wait(dual->barrierPointer()); #if 0 int status=-1; while (status!=100) status=timedWait(dual,1000,2); pthread_cond_signal(dual->conditionPointer(1)); pthread_mutex_unlock(dual->mutexPointer(1,whichThread)); #endif // so now mutex_ is locked int whichLocked = 0; while (true) { pthread_mutex_t *mutexPointer = dual->mutexPointer(whichLocked, whichThread); // wait //printf("Child waiting for %d - status %d %d %d\n", // whichLocked,lockedX[0],lockedX[1],lockedX[2]); pthread_mutex_lock(mutexPointer); whichLocked++; if (whichLocked == 3) whichLocked = 0; int unLock = whichLocked + 1; if (unLock == 3) unLock = 0; //printf("child pointer %x status %d\n",threadInfo,threadInfo->status); assert(threadInfo->status >= 0); if (threadInfo->status == 1000) pthread_exit(NULL); int type = threadInfo->status; int *which = threadInfo->stuff; int &returnCode = which[0]; CoinIndexedVector *array; double dummy; switch (type) { // dummy case 0: break; case 1: returnCode = computeDualsAndCheck(dual, which); break; case 2: array = dual->usefulArray(dual->arrayForBtran()); // assuming not big overlap dual->dualRowPivot()->updateWeightsOnly(*array); //array->checkClean(); break; case 3: //dual->replaceColumnPart3(); dual->updateDualsInDual(); break; case 4: which[2] = -1; dual->abcMatrix()->dualColumn1Part(which[1], which[2], threadInfo->result, *dual->usefulArray(dual->arrayForBtran()), *dual->usefulArray(dual->arrayForTableauRow()), *dual->usefulArray(dual->arrayForDualColumn())); break; case 5: which[2] = -1; threadInfo->result = dual->abcMatrix()->dualColumn1Row(which[1], COIN_DBL_MAX, which[2], *dual->usefulArray(dual->arrayForBtran()), *dual->usefulArray(dual->arrayForTableauRow()), *dual->usefulArray(dual->arrayForDualColumn())); break; case 6: // not used dual->getTableauColumnPart2(); break; case 7: dual->checkReplacePart1(); //dual->checkReplacePart1a(); //dual->checkReplacePart1b(); break; case 8: which[2] = dual->flipBounds(); break; case 9: dual->replaceColumnPart3(); break; case 100: // initialization break; } threadInfo->status = -1; pthread_mutex_unlock(dual->mutexPointer(unLock, whichThread)); } } void AbcSimplex::startParallelStuff(int type) { /* first time 0,1 owned by main 2 by child at end of cycle should be 1,2 by main 0 by child then 2,0 by main 1 by child */ int iStart = 1; for (int i = 0; i < NUMBER_THREADS; i++) { if ((stopStart_ & iStart) != 0) { threadInfo_[i].status = type; pthread_mutex_unlock(&mutex_[locked_[i] + 3 * i]); } iStart = iStart << 1; } } #define NUMBER_PARALLEL_TASKS 10 //static int count_x[NUMBER_PARALLEL_TASKS]={0,0,0,0,0,0,0,0,0,0}; //static int count_y[NUMBER_PARALLEL_TASKS]={0,0,0,0,0,0,0,0,0,0}; int AbcSimplex::stopParallelStuff(int type) { //if (threadInfo_[0].status==-1) //count_y[type-1]++; //count_x[type-1]++; //printf("XXStop type %d iteration %d\n",type,numberIterations_); //if (threadInfo_[0].status!=-1) { //printf("Main waiting for %d - status %d\n", // locked,type); //pthread_mutex_lock(&mutex_[locked]); int iStop = 32; for (int i = 0; i < NUMBER_THREADS; i++) { if ((stopStart_ & iStop) != 0) { int locked = locked_[i] + 2; if (locked >= 3) locked -= 3; pthread_mutex_lock(&mutex_[locked + i * 3]); locked_[i]++; if (locked_[i] == 3) locked_[i] = 0; } iStop = iStop << 1; } return threadInfo_[0].stuff[0]; } #endif //#define PRINT_WEIGHTS_FREQUENCY #ifdef PRINT_WEIGHTS_FREQUENCY extern int xweights1; extern int xweights2; extern int xweights3; #endif //static int whichModel=0; int AbcSimplexDual::dual() { //handler_->setLogLevel(63); #if 0 whichModel++; printf("STARTing dual %d rows - model %d\n",numberRows_,whichModel); char modelName[20]; sprintf(modelName,"model%d.mps",whichModel); writeMps(modelName,2); sprintf(modelName,"model%d.bas",whichModel); writeBasis(modelName,true,2); #endif bestObjectiveValue_ = -COIN_DBL_MAX; algorithm_ = -1; moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy normalDualColumnIteration_ = 0; problemStatus_ = -1; numberIterations_ = baseIteration_; //baseIteration_=0; if (!lastDualBound_) { dualTolerance_ = dblParam_[ClpDualTolerance]; primalTolerance_ = dblParam_[ClpPrimalTolerance]; int tryType = moreSpecialOptions_ / 65536; if (tryType < 5) { currentDualBound_ = 1.0e8; } else { currentDualBound_ = 1.0e4; tryType -= 5; } currentDualBound_ = CoinMin(currentDualBound_, dualBound_); #if ABC_NORMAL_DEBUG > 0 printf("dual bound %g (initialize)\n", currentDualBound_); #endif switch (tryType) { default: currentDualTolerance_ = dualTolerance_; stateDualColumn_ = 0; break; case 1: currentDualTolerance_ = CoinMin(1.0e-1, dualTolerance_ * 10.0); stateDualColumn_ = -1; break; case 2: currentDualTolerance_ = CoinMin(1.0e-1, dualTolerance_ * 100.0); stateDualColumn_ = -1; break; case 3: currentDualTolerance_ = CoinMin(1.0e-1, dualTolerance_ * 10000.0); stateDualColumn_ = -1; // was -2 but that doesn't seem to work break; case 4: currentDualTolerance_ = CoinMin(1.0e-1, dualTolerance_ * 1000000.0); stateDualColumn_ = -1; // was -2 but that doesn't seem to work break; } } else { // done before currentDualBound_ = 1.1 * lastDualBound_; currentDualTolerance_ = dualTolerance_; #if ABC_NORMAL_DEBUG > 0 printf("new dual bound %g (resolve)\n", currentDualBound_); #endif } // save data saveData_ = saveData(); // initialize copyFromSaved(); // temp if (fabs(primalTolerance() - 1.001e-6) < 1.0e-12 && false) { printf("trying idiot idea - need scaled matrix\n"); primalTolerance_ = 1.0e-6; //for (int i=0;itransposeTimesAll(cost_, abcCost_); delete[] cost_; cost_ = NULL; CoinAbcMemcpy(costSaved_, abcCost_, numberTotal_); CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_); static_cast< AbcSimplexDual * >(this)->makeNonFreeVariablesDualFeasible(); const double multiplier[] = { 1.0, -1.0 }; double dualT = -currentDualTolerance_; double *COIN_RESTRICT solution = abcSolution_; double *COIN_RESTRICT lower = abcLower_; double *COIN_RESTRICT upper = abcUpper_; for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { int iStatus = internalStatus_[iSequence] & 7; if (iStatus < 4) { double lowerValue = lower[iSequence]; double upperValue = upper[iSequence]; if (lowerValue > -1.0e10 || upperValue < 1.0e10) { double mult = multiplier[iStatus]; double dj = abcDj_[iSequence] * mult; if (dj < dualT) iStatus = 1 - iStatus; if (lowerValue <= -1.0e10) iStatus = 1; else if (upperValue >= 1.0e10) iStatus = 0; if (iStatus == 0) { setInternalStatus(iSequence, atLowerBound); solution[iSequence] = lower[iSequence]; } else { setInternalStatus(iSequence, atUpperBound); solution[iSequence] = upper[iSequence]; } } assert(fabs(solution[iSequence]) < 1.0e10); } } CoinAbcMemcpy(solutionSaved_, abcSolution_, numberTotal_); } numberFlagged_ = 0; #ifndef NDEBUG for (int iRow = 0; iRow < numberRows_; iRow++) { int iPivot = abcPivotVariable_[iRow]; assert(!flagged(iPivot)); } #endif #if ABC_PARALLEL == 1 // redo to test on parallelMode_ if (parallelMode_ != 0) { // For waking up thread memset(mutex_, 0, sizeof(mutex_)); for (int iThread = 0; iThread < NUMBER_THREADS; iThread++) { for (int i = 0; i < 3; i++) { pthread_mutex_init(&mutex_[i + 3 * iThread], NULL); if (i < 2) pthread_mutex_lock(&mutex_[i + 3 * iThread]); } threadInfo_[iThread].status = 100; } //pthread_barrierattr_t attr; pthread_barrier_init(&barrier_, /*&attr*/ NULL, NUMBER_THREADS + 1); for (int iThread = 0; iThread < NUMBER_THREADS; iThread++) { pthread_create(&abcThread_[iThread], NULL, abc_parallelManager, reinterpret_cast< void * >(this)); } pthread_barrier_wait(&barrier_); pthread_barrier_destroy(&barrier_); for (int iThread = 0; iThread < NUMBER_THREADS; iThread++) { threadInfo_[iThread].status = 0; } for (int i = 0; i < NUMBER_THREADS; i++) { locked_[i] = 0; } } #endif #if ABC_PARALLEL == 2 abcFactorization_->setParallelMode(parallelMode_); #endif static_cast< AbcSimplexDual * >(this)->bounceTolerances(-1); int saveDont = dontFactorizePivots_; if ((specialOptions_ & 2048) == 0) dontFactorizePivots_ = 0; else if (!dontFactorizePivots_) dontFactorizePivots_ = 20; if (alphaAccuracy_ != -1.0) alphaAccuracy_ = 1.0; startupSolve(); lastPivotRow_ = -1; pivotRow_ = -1; // Save so can see if doing after primal swappedAlgorithm_ = (problemStatus_ == 10) ? 10 : 0; if (problemStatus_) { #ifdef PRINT_WEIGHTS_FREQUENCY xweights1 = 0; xweights2 = 0; xweights3 = 0; #endif numberFreeNonBasic_ = 0; gutsOfDual(); #ifdef PRINT_WEIGHTS_FREQUENCY if (dynamic_cast< AbcDualRowSteepest * >(abcDualRowPivot_)) { printf("ZZZZ %d iterations weights updated %d times\n", numberIterations_, xweights1 + xweights2 + xweights3); assert(numberIterations_ == xweights1 + xweights2 + xweights3); } #endif } if (!problemStatus_) { // see if cutoff reached double limit = 0.0; getDblParam(ClpDualObjectiveLimit, limit); if (fabs(limit) < 1.0e30 && minimizationObjectiveValue() > limit + 1.0e-7 + 1.0e-8 * fabs(limit)) { // actually infeasible on objective problemStatus_ = 1; secondaryStatus_ = 1; } } // If infeasible but primal errors - try dual if (problemStatus_ == 1 && numberPrimalInfeasibilities_) { bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0; double factor = (!inCbcOrOther) ? 1.0 : 0.3; double averageInfeasibility = sumPrimalInfeasibilities_ / static_cast< double >(numberPrimalInfeasibilities_); if (averageInfeasibility < factor * largestPrimalError_) problemStatus_ = 10; } finishSolve(); // Restore any saved stuff restoreData(saveData_); dontFactorizePivots_ = saveDont; if (problemStatus_ == 3) { rawObjectiveValue_ = CoinMax(bestObjectiveValue_, rawObjectiveValue_ - bestPossibleImprovement_); setClpSimplexObjectiveValue(); } #if ABC_PARALLEL == 1 if (parallelMode_ != 0) { //for (int i=0;i