[1910] | 1 | /* $Id: AbcSimplex.cpp 2429 2019-03-11 16:34:25Z stefan $ */ |
---|
[1878] | 2 | // Copyright (C) 2000, International Business Machines |
---|
| 3 | // Corporation and others, Copyright (C) 2012, FasterCoin. All Rights Reserved. |
---|
| 4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
| 5 | |
---|
| 6 | #include "ClpConfig.h" |
---|
| 7 | |
---|
| 8 | #include "CoinPragma.hpp" |
---|
| 9 | #include <math.h> |
---|
| 10 | //#define ABC_DEBUG 2 |
---|
| 11 | |
---|
[2385] | 12 | #if SLIM_CLP == 2 |
---|
[1878] | 13 | #define SLIM_NOIO |
---|
| 14 | #endif |
---|
| 15 | #include "CoinHelperFunctions.hpp" |
---|
| 16 | #include "CoinFloatEqual.hpp" |
---|
| 17 | #include "ClpSimplex.hpp" |
---|
| 18 | #include "AbcSimplex.hpp" |
---|
| 19 | #include "AbcSimplexDual.hpp" |
---|
| 20 | #include "AbcSimplexFactorization.hpp" |
---|
| 21 | #include "AbcNonLinearCost.hpp" |
---|
| 22 | #include "CoinAbcCommon.hpp" |
---|
| 23 | #include "AbcMatrix.hpp" |
---|
| 24 | #include "CoinIndexedVector.hpp" |
---|
| 25 | #include "AbcDualRowDantzig.hpp" |
---|
| 26 | #include "AbcDualRowSteepest.hpp" |
---|
| 27 | #include "AbcPrimalColumnDantzig.hpp" |
---|
| 28 | #include "AbcPrimalColumnSteepest.hpp" |
---|
| 29 | #include "ClpMessage.hpp" |
---|
| 30 | #include "ClpEventHandler.hpp" |
---|
| 31 | #include "ClpLinearObjective.hpp" |
---|
| 32 | #include "CoinAbcHelperFunctions.hpp" |
---|
| 33 | #include "CoinModel.hpp" |
---|
| 34 | #include "CoinLpIO.hpp" |
---|
| 35 | #include <cfloat> |
---|
| 36 | |
---|
| 37 | #include <string> |
---|
| 38 | #include <stdio.h> |
---|
| 39 | #include <iostream> |
---|
| 40 | //############################################################################# |
---|
[2385] | 41 | AbcSimplex::AbcSimplex(bool emptyMessages) |
---|
| 42 | : |
---|
| 43 | |
---|
[1878] | 44 | ClpSimplex(emptyMessages) |
---|
| 45 | { |
---|
| 46 | gutsOfDelete(0); |
---|
[2385] | 47 | gutsOfInitialize(0, 0, true); |
---|
| 48 | assert(maximumAbcNumberRows_ >= 0); |
---|
[1878] | 49 | //printf("XX %x AbcSimplex constructor\n",this); |
---|
| 50 | } |
---|
| 51 | |
---|
| 52 | //----------------------------------------------------------------------------- |
---|
| 53 | |
---|
[2385] | 54 | AbcSimplex::~AbcSimplex() |
---|
[1878] | 55 | { |
---|
| 56 | //printf("XX %x AbcSimplex destructor\n",this); |
---|
| 57 | gutsOfDelete(1); |
---|
| 58 | } |
---|
| 59 | // Copy constructor. |
---|
[2385] | 60 | AbcSimplex::AbcSimplex(const AbcSimplex &rhs) |
---|
| 61 | : ClpSimplex(rhs) |
---|
[1878] | 62 | { |
---|
| 63 | //printf("XX %x AbcSimplex constructor from %x\n",this,&rhs); |
---|
| 64 | gutsOfDelete(0); |
---|
[2385] | 65 | gutsOfInitialize(numberRows_, numberColumns_, false); |
---|
[1878] | 66 | gutsOfCopy(rhs); |
---|
[2385] | 67 | assert(maximumAbcNumberRows_ >= 0); |
---|
[1878] | 68 | } |
---|
| 69 | #include "ClpDualRowSteepest.hpp" |
---|
| 70 | #include "ClpPrimalColumnSteepest.hpp" |
---|
| 71 | #include "ClpFactorization.hpp" |
---|
| 72 | // Copy constructor from model |
---|
[2385] | 73 | AbcSimplex::AbcSimplex(const ClpSimplex &rhs) |
---|
| 74 | : ClpSimplex(rhs) |
---|
[1878] | 75 | { |
---|
| 76 | //printf("XX %x AbcSimplex constructor from ClpSimplex\n",this); |
---|
| 77 | gutsOfDelete(0); |
---|
[2385] | 78 | gutsOfInitialize(numberRows_, numberColumns_, true); |
---|
| 79 | gutsOfResize(numberRows_, numberColumns_); |
---|
[1878] | 80 | // Set up row/column selection and factorization type |
---|
[2385] | 81 | ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(rhs.dualRowPivot()); |
---|
[1878] | 82 | if (pivot) { |
---|
| 83 | AbcDualRowSteepest steep(pivot->mode()); |
---|
| 84 | setDualRowPivotAlgorithm(steep); |
---|
| 85 | } else { |
---|
| 86 | AbcDualRowDantzig dantzig; |
---|
| 87 | setDualRowPivotAlgorithm(dantzig); |
---|
| 88 | } |
---|
[2385] | 89 | ClpPrimalColumnSteepest *pivotColumn = dynamic_cast< ClpPrimalColumnSteepest * >(rhs.primalColumnPivot()); |
---|
[1878] | 90 | if (pivotColumn) { |
---|
| 91 | AbcPrimalColumnSteepest steep(pivotColumn->mode()); |
---|
| 92 | setPrimalColumnPivotAlgorithm(steep); |
---|
| 93 | } else { |
---|
| 94 | AbcPrimalColumnDantzig dantzig; |
---|
| 95 | setPrimalColumnPivotAlgorithm(dantzig); |
---|
| 96 | } |
---|
| 97 | //if (rhs.factorization()->) |
---|
| 98 | //factorization_->forceOtherFactorization(); |
---|
[2385] | 99 | factorization()->synchronize(rhs.factorization(), this); |
---|
[1878] | 100 | //factorization_->setGoDenseThreshold(rhs.factorization()->goDenseThreshold()); |
---|
| 101 | //factorization_->setGoSmallThreshold(rhs.factorization()->goSmallThreshold()); |
---|
[2385] | 102 | translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION); |
---|
| 103 | assert(maximumAbcNumberRows_ >= 0); |
---|
[1878] | 104 | } |
---|
| 105 | // Assignment operator. This copies the data |
---|
| 106 | AbcSimplex & |
---|
[2385] | 107 | AbcSimplex::operator=(const AbcSimplex &rhs) |
---|
[1878] | 108 | { |
---|
| 109 | if (this != &rhs) { |
---|
| 110 | gutsOfDelete(1); |
---|
| 111 | ClpSimplex::operator=(rhs); |
---|
| 112 | gutsOfCopy(rhs); |
---|
[2385] | 113 | assert(maximumAbcNumberRows_ >= 0); |
---|
[1878] | 114 | } |
---|
| 115 | return *this; |
---|
| 116 | } |
---|
| 117 | // fills in perturbationSaved_ from start with 0.5+random |
---|
[2385] | 118 | void AbcSimplex::fillPerturbation(int start, int number) |
---|
[1878] | 119 | { |
---|
[2385] | 120 | double *array = perturbationSaved_ + start; |
---|
| 121 | for (int i = 0; i < number; i++) |
---|
| 122 | array[i] = 0.5 + 0.5 * randomNumberGenerator_.randomDouble(); |
---|
[1878] | 123 | } |
---|
| 124 | // Sets up all extra pointers |
---|
[2385] | 125 | void AbcSimplex::setupPointers(int maxRows, int maxColumns) |
---|
[1878] | 126 | { |
---|
[2385] | 127 | int numberTotal = maxRows + maxColumns; |
---|
| 128 | scaleToExternal_ = scaleFromExternal_ + numberTotal; |
---|
| 129 | tempArray_ = offset_ + numberTotal; |
---|
| 130 | lowerSaved_ = abcLower_ + numberTotal; |
---|
| 131 | upperSaved_ = abcUpper_ + numberTotal; |
---|
| 132 | costSaved_ = abcCost_ + numberTotal; |
---|
| 133 | solutionSaved_ = abcSolution_ + numberTotal; |
---|
| 134 | djSaved_ = abcDj_ + numberTotal; |
---|
| 135 | internalStatusSaved_ = internalStatus_ + numberTotal; |
---|
| 136 | perturbationSaved_ = abcPerturbation_ + numberTotal; |
---|
| 137 | offsetRhs_ = tempArray_ + numberTotal; |
---|
| 138 | lowerBasic_ = lowerSaved_ + numberTotal; |
---|
| 139 | upperBasic_ = upperSaved_ + numberTotal; |
---|
| 140 | costBasic_ = costSaved_ + 2 * numberTotal; |
---|
| 141 | solutionBasic_ = solutionSaved_ + numberTotal; |
---|
| 142 | djBasic_ = djSaved_ + numberTotal; |
---|
| 143 | perturbationBasic_ = perturbationSaved_ + numberTotal; |
---|
| 144 | columnUseScale_ = scaleToExternal_ + maxRows; |
---|
| 145 | inverseColumnUseScale_ = scaleFromExternal_ + maxRows; |
---|
[1878] | 146 | } |
---|
| 147 | // Copies all saved versions to working versions and may do something for perturbation |
---|
[2385] | 148 | void AbcSimplex::copyFromSaved(int which) |
---|
[1878] | 149 | { |
---|
[2385] | 150 | if ((which & 1) != 0) |
---|
| 151 | CoinAbcMemcpy(abcSolution_, solutionSaved_, maximumNumberTotal_); |
---|
| 152 | if ((which & 2) != 0) |
---|
| 153 | CoinAbcMemcpy(abcCost_, costSaved_, maximumNumberTotal_); |
---|
| 154 | if ((which & 4) != 0) |
---|
| 155 | CoinAbcMemcpy(abcLower_, lowerSaved_, maximumNumberTotal_); |
---|
| 156 | if ((which & 8) != 0) |
---|
| 157 | CoinAbcMemcpy(abcUpper_, upperSaved_, maximumNumberTotal_); |
---|
| 158 | if ((which & 16) != 0) |
---|
| 159 | CoinAbcMemcpy(abcDj_, djSaved_, maximumNumberTotal_); |
---|
| 160 | if ((which & 32) != 0) { |
---|
| 161 | CoinAbcMemcpy(abcLower_, abcPerturbation_, numberTotal_); |
---|
| 162 | CoinAbcMemcpy(abcUpper_, abcPerturbation_ + numberTotal_, numberTotal_); |
---|
[1878] | 163 | } |
---|
| 164 | } |
---|
[2385] | 165 | void AbcSimplex::gutsOfCopy(const AbcSimplex &rhs) |
---|
[1878] | 166 | { |
---|
[1910] | 167 | #ifdef ABC_INHERIT |
---|
[2385] | 168 | abcSimplex_ = this; |
---|
[1910] | 169 | #endif |
---|
[2385] | 170 | assert(numberRows_ == rhs.numberRows_); |
---|
| 171 | assert(numberColumns_ == rhs.numberColumns_); |
---|
[1878] | 172 | numberTotal_ = rhs.numberTotal_; |
---|
| 173 | maximumNumberTotal_ = rhs.maximumNumberTotal_; |
---|
| 174 | // special options here to be safe |
---|
[2385] | 175 | specialOptions_ = rhs.specialOptions_; |
---|
[1878] | 176 | //assert (maximumInternalRows_ >= numberRows_); |
---|
| 177 | //assert (maximumInternalColumns_ >= numberColumns_); |
---|
| 178 | // Copy all scalars |
---|
[2385] | 179 | CoinAbcMemcpy(reinterpret_cast< int * >(&sumNonBasicCosts_), |
---|
| 180 | reinterpret_cast< const int * >(&rhs.sumNonBasicCosts_), |
---|
| 181 | (&swappedAlgorithm_ - reinterpret_cast< int * >(&sumNonBasicCosts_)) + 1); |
---|
[1878] | 182 | // could add 2 so can go off end |
---|
[2385] | 183 | int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_; |
---|
| 184 | internalStatus_ = ClpCopyOfArray(rhs.internalStatus_, sizeArray + maximumNumberTotal_); |
---|
[1878] | 185 | abcLower_ = ClpCopyOfArray(rhs.abcLower_, sizeArray); |
---|
| 186 | abcUpper_ = ClpCopyOfArray(rhs.abcUpper_, sizeArray); |
---|
[2385] | 187 | abcCost_ = ClpCopyOfArray(rhs.abcCost_, sizeArray + maximumNumberTotal_); |
---|
[1878] | 188 | abcDj_ = ClpCopyOfArray(rhs.abcDj_, sizeArray); |
---|
[1990] | 189 | |
---|
[2385] | 190 | abcSolution_ = ClpCopyOfArray(rhs.abcSolution_, sizeArray + maximumNumberTotal_); |
---|
| 191 | abcPerturbation_ = ClpCopyOfArray(rhs.abcPerturbation_, sizeArray); |
---|
| 192 | abcPivotVariable_ = ClpCopyOfArray(rhs.abcPivotVariable_, maximumAbcNumberRows_); |
---|
[1878] | 193 | //fromExternal_ = ClpCopyOfArray(rhs.fromExternal_,sizeArray); |
---|
| 194 | //toExternal_ = ClpCopyOfArray(rhs.toExternal_,sizeArray); |
---|
[2385] | 195 | scaleFromExternal_ = ClpCopyOfArray(rhs.scaleFromExternal_, sizeArray); |
---|
| 196 | offset_ = ClpCopyOfArray(rhs.offset_, sizeArray); |
---|
| 197 | setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_); |
---|
[1878] | 198 | if (rhs.abcMatrix_) |
---|
| 199 | abcMatrix_ = new AbcMatrix(*rhs.abcMatrix_); |
---|
| 200 | else |
---|
| 201 | abcMatrix_ = NULL; |
---|
| 202 | for (int i = 0; i < ABC_NUMBER_USEFUL; i++) { |
---|
| 203 | usefulArray_[i] = rhs.usefulArray_[i]; |
---|
| 204 | } |
---|
| 205 | if (rhs.abcFactorization_) { |
---|
| 206 | setFactorization(*rhs.abcFactorization_); |
---|
| 207 | } else { |
---|
| 208 | delete abcFactorization_; |
---|
| 209 | abcFactorization_ = NULL; |
---|
| 210 | } |
---|
| 211 | #ifdef EARLY_FACTORIZE |
---|
| 212 | delete abcEarlyFactorization_; |
---|
| 213 | if (rhs.abcEarlyFactorization_) { |
---|
| 214 | abcEarlyFactorization_ = new AbcSimplexFactorization(*rhs.abcEarlyFactorization_); |
---|
| 215 | } else { |
---|
| 216 | abcEarlyFactorization_ = NULL; |
---|
| 217 | } |
---|
| 218 | #endif |
---|
| 219 | abcDualRowPivot_ = rhs.abcDualRowPivot_->clone(true); |
---|
| 220 | abcDualRowPivot_->setModel(this); |
---|
| 221 | abcPrimalColumnPivot_ = rhs.abcPrimalColumnPivot_->clone(true); |
---|
| 222 | abcPrimalColumnPivot_->setModel(this); |
---|
| 223 | if (rhs.abcBaseModel_) { |
---|
| 224 | abcBaseModel_ = new AbcSimplex(*rhs.abcBaseModel_); |
---|
| 225 | } else { |
---|
| 226 | abcBaseModel_ = NULL; |
---|
| 227 | } |
---|
| 228 | if (rhs.clpModel_) { |
---|
| 229 | clpModel_ = new ClpSimplex(*rhs.clpModel_); |
---|
| 230 | } else { |
---|
| 231 | clpModel_ = NULL; |
---|
| 232 | } |
---|
| 233 | abcProgress_ = rhs.abcProgress_; |
---|
| 234 | solveType_ = rhs.solveType_; |
---|
| 235 | } |
---|
| 236 | // type == 0 nullify, 1 also delete |
---|
[2385] | 237 | void AbcSimplex::gutsOfDelete(int type) |
---|
[1878] | 238 | { |
---|
| 239 | if (type) { |
---|
[2385] | 240 | delete[] abcLower_; |
---|
| 241 | delete[] abcUpper_; |
---|
| 242 | delete[] abcCost_; |
---|
| 243 | delete[] abcDj_; |
---|
| 244 | delete[] abcSolution_; |
---|
[1878] | 245 | //delete [] fromExternal_ ; |
---|
| 246 | //delete [] toExternal_ ; |
---|
[2385] | 247 | delete[] scaleFromExternal_; |
---|
[1878] | 248 | //delete [] scaleToExternal_ ; |
---|
[2385] | 249 | delete[] offset_; |
---|
| 250 | delete[] internalStatus_; |
---|
| 251 | delete[] abcPerturbation_; |
---|
| 252 | delete abcMatrix_; |
---|
[1878] | 253 | delete abcFactorization_; |
---|
| 254 | #ifdef EARLY_FACTORIZE |
---|
| 255 | delete abcEarlyFactorization_; |
---|
| 256 | #endif |
---|
[2385] | 257 | delete[] abcPivotVariable_; |
---|
[1878] | 258 | delete abcDualRowPivot_; |
---|
| 259 | delete abcPrimalColumnPivot_; |
---|
| 260 | delete abcBaseModel_; |
---|
| 261 | delete clpModel_; |
---|
| 262 | delete abcNonLinearCost_; |
---|
| 263 | } |
---|
[2385] | 264 | CoinAbcMemset0(reinterpret_cast< char * >(&scaleToExternal_), |
---|
| 265 | reinterpret_cast< char * >(&usefulArray_[0]) - reinterpret_cast< char * >(&scaleToExternal_)); |
---|
[1878] | 266 | } |
---|
[2385] | 267 | template < class T > |
---|
| 268 | T *newArray(T * /*nullArray*/, int size) |
---|
[1878] | 269 | { |
---|
| 270 | if (size) { |
---|
[2385] | 271 | T *arrayNew = new T[size]; |
---|
[1878] | 272 | #ifndef NDEBUG |
---|
[2385] | 273 | memset(arrayNew, 'A', (size) * sizeof(T)); |
---|
[1878] | 274 | #endif |
---|
| 275 | return arrayNew; |
---|
| 276 | } else { |
---|
| 277 | return NULL; |
---|
| 278 | } |
---|
| 279 | } |
---|
[2385] | 280 | template < class T > |
---|
| 281 | T *resizeArray(T *array, int oldSize1, int oldSize2, int newSize1, int newSize2, int extra) |
---|
[1878] | 282 | { |
---|
[2385] | 283 | if ((array || !oldSize1) && (newSize1 != oldSize1 || newSize2 != oldSize2)) { |
---|
| 284 | int newTotal = newSize1 + newSize2; |
---|
| 285 | int oldTotal = oldSize1 + oldSize2; |
---|
| 286 | T *arrayNew; |
---|
| 287 | if (newSize1 > oldSize1 || newSize2 > oldSize2) { |
---|
| 288 | arrayNew = new T[2 * newTotal + newSize1 + extra]; |
---|
[1878] | 289 | #ifndef NDEBUG |
---|
[2385] | 290 | memset(arrayNew, 'A', (2 * newTotal + newSize1 + extra) * sizeof(T)); |
---|
[1878] | 291 | #endif |
---|
[2385] | 292 | CoinAbcMemcpy(arrayNew, array, oldSize1); |
---|
| 293 | CoinAbcMemcpy(arrayNew + newSize1, array + oldSize1, oldSize2); |
---|
[1878] | 294 | // and second half |
---|
[2385] | 295 | CoinAbcMemcpy(arrayNew + newTotal, array, oldSize1 + oldTotal); |
---|
| 296 | CoinAbcMemcpy(arrayNew + newSize1 + newTotal, array + oldSize1 + oldTotal, oldSize2); |
---|
| 297 | delete[] array; |
---|
[1878] | 298 | } else { |
---|
[2385] | 299 | arrayNew = array; |
---|
| 300 | for (int i = 0; i < newSize2; i++) |
---|
| 301 | array[newSize1 + i] = array[oldSize1 + i]; |
---|
| 302 | for (int i = 0; i < newSize1; i++) |
---|
| 303 | array[newTotal + i] = array[oldTotal + i]; |
---|
| 304 | for (int i = 0; i < newSize2; i++) |
---|
| 305 | array[newSize1 + newTotal + i] = array[oldSize1 + oldTotal + i]; |
---|
[1878] | 306 | } |
---|
| 307 | return arrayNew; |
---|
| 308 | } else { |
---|
| 309 | return array; |
---|
| 310 | } |
---|
| 311 | } |
---|
| 312 | // Initializes arrays |
---|
[2385] | 313 | void AbcSimplex::gutsOfInitialize(int numberRows, int numberColumns, bool doMore) |
---|
[1878] | 314 | { |
---|
[1910] | 315 | #ifdef ABC_INHERIT |
---|
[2385] | 316 | abcSimplex_ = this; |
---|
[1910] | 317 | #endif |
---|
[1878] | 318 | // Zero all |
---|
[2385] | 319 | CoinAbcMemset0(&sumNonBasicCosts_, (reinterpret_cast< double * >(&usefulArray_[0]) - &sumNonBasicCosts_)); |
---|
[1878] | 320 | zeroTolerance_ = 1.0e-13; |
---|
| 321 | bestObjectiveValue_ = -COIN_DBL_MAX; |
---|
| 322 | primalToleranceToGetOptimal_ = -1.0; |
---|
| 323 | primalTolerance_ = 1.0e-6; |
---|
| 324 | //dualTolerance_ = 1.0e-6; |
---|
| 325 | alphaAccuracy_ = -1.0; |
---|
| 326 | upperIn_ = -COIN_DBL_MAX; |
---|
| 327 | lowerOut_ = -1; |
---|
| 328 | valueOut_ = -1; |
---|
| 329 | upperOut_ = -1; |
---|
| 330 | dualOut_ = -1; |
---|
| 331 | acceptablePivot_ = 1.0e-8; |
---|
| 332 | //dualBound_=1.0e9; |
---|
| 333 | sequenceIn_ = -1; |
---|
| 334 | directionIn_ = -1; |
---|
| 335 | sequenceOut_ = -1; |
---|
| 336 | directionOut_ = -1; |
---|
| 337 | pivotRow_ = -1; |
---|
| 338 | lastGoodIteration_ = -100; |
---|
| 339 | numberPrimalInfeasibilities_ = 100; |
---|
[2385] | 340 | numberRefinements_ = 1; |
---|
[1878] | 341 | changeMade_ = 1; |
---|
| 342 | forceFactorization_ = -1; |
---|
[2385] | 343 | if (perturbation_ < 50 || (perturbation_ > 60 && perturbation_ < 100)) |
---|
[1878] | 344 | perturbation_ = 50; |
---|
| 345 | lastBadIteration_ = -999999; |
---|
| 346 | lastFlaggedIteration_ = -999999; // doesn't seem to be used |
---|
| 347 | firstFree_ = -1; |
---|
| 348 | incomingInfeasibility_ = 1.0; |
---|
| 349 | allowedInfeasibility_ = 10.0; |
---|
| 350 | solveType_ = 1; // say simplex based life form |
---|
| 351 | //specialOptions_|=65536; |
---|
| 352 | //ClpModel::startPermanentArrays(); |
---|
[2385] | 353 | maximumInternalRows_ = 0; |
---|
| 354 | maximumInternalColumns_ = 0; |
---|
| 355 | numberRows_ = numberRows; |
---|
| 356 | numberColumns_ = numberColumns; |
---|
| 357 | numberTotal_ = numberRows_ + numberColumns_; |
---|
| 358 | maximumAbcNumberRows_ = numberRows; |
---|
| 359 | maximumAbcNumberColumns_ = numberColumns; |
---|
| 360 | maximumNumberTotal_ = numberTotal_; |
---|
| 361 | int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_; |
---|
[1878] | 362 | if (doMore) { |
---|
| 363 | // say Steepest pricing |
---|
| 364 | abcDualRowPivot_ = new AbcDualRowSteepest(); |
---|
| 365 | abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest(); |
---|
[2429] | 366 | internalStatus_ = newArray((unsigned char *)NULL, |
---|
[2385] | 367 | sizeArray + maximumNumberTotal_); |
---|
[2429] | 368 | abcLower_ = newArray((double *)NULL, sizeArray); |
---|
| 369 | abcUpper_ = newArray((double *)NULL, sizeArray); |
---|
| 370 | abcCost_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_); |
---|
| 371 | abcDj_ = newArray((double *)NULL, sizeArray); |
---|
| 372 | abcSolution_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_); |
---|
| 373 | //fromExternal_ = newArray((int *)NULL,sizeArray); |
---|
| 374 | //toExternal_ = newArray((int *)NULL,sizeArray); |
---|
| 375 | scaleFromExternal_ = newArray((double *)NULL, sizeArray); |
---|
| 376 | offset_ = newArray((double *)NULL, sizeArray); |
---|
| 377 | abcPerturbation_ = newArray((double *)NULL, sizeArray); |
---|
| 378 | abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_); |
---|
[1878] | 379 | // Fill perturbation array |
---|
[2385] | 380 | setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_); |
---|
| 381 | fillPerturbation(0, maximumNumberTotal_); |
---|
[1878] | 382 | } |
---|
| 383 | // get an empty factorization so we can set tolerances etc |
---|
| 384 | getEmptyFactorization(); |
---|
[2385] | 385 | for (int i = 0; i < ABC_NUMBER_USEFUL; i++) |
---|
| 386 | usefulArray_[i].reserve(CoinMax(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200), 2 * numberRows_)); |
---|
[1878] | 387 | //savedStatus_ = internalStatus_+maximumNumberTotal_; |
---|
| 388 | //startPermanentArrays(); |
---|
| 389 | } |
---|
| 390 | // resizes arrays |
---|
[2385] | 391 | void AbcSimplex::gutsOfResize(int numberRows, int numberColumns) |
---|
[1878] | 392 | { |
---|
| 393 | if (!abcSolution_) { |
---|
[2385] | 394 | numberRows_ = 0; |
---|
| 395 | numberColumns_ = 0; |
---|
| 396 | numberTotal_ = 0; |
---|
| 397 | maximumAbcNumberRows_ = 0; |
---|
| 398 | maximumAbcNumberColumns_ = 0; |
---|
| 399 | maximumNumberTotal_ = 0; |
---|
[1878] | 400 | } |
---|
[2385] | 401 | if (numberRows == numberRows_ && numberColumns == numberColumns_) |
---|
[1878] | 402 | return; |
---|
| 403 | // can do on state bit |
---|
[2385] | 404 | int newSize1 = CoinMax(numberRows, maximumAbcNumberRows_); |
---|
| 405 | if ((stateOfProblem_ & ADD_A_BIT) != 0 && numberRows > maximumAbcNumberRows_) |
---|
| 406 | newSize1 = CoinMax(numberRows, maximumAbcNumberRows_ + CoinMin(100, numberRows_ / 10)); |
---|
| 407 | int newSize2 = CoinMax(numberColumns, maximumAbcNumberColumns_); |
---|
| 408 | numberRows_ = numberRows; |
---|
| 409 | numberColumns_ = numberColumns; |
---|
| 410 | numberTotal_ = numberRows_ + numberColumns_; |
---|
[1878] | 411 | //fromExternal_ = resizeArray(fromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1); |
---|
| 412 | //toExternal_ = resizeArray(toExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0); |
---|
[2385] | 413 | scaleFromExternal_ = resizeArray(scaleFromExternal_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
[1878] | 414 | //scaleToExternal_ = resizeArray(scaleToExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0); |
---|
[2385] | 415 | internalStatus_ = resizeArray(internalStatus_, maximumAbcNumberRows_, |
---|
| 416 | maximumAbcNumberColumns_, |
---|
| 417 | newSize1, newSize2, numberTotal_); |
---|
| 418 | abcLower_ = resizeArray(abcLower_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
| 419 | abcUpper_ = resizeArray(abcUpper_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
| 420 | abcCost_ = resizeArray(abcCost_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_); |
---|
| 421 | abcDj_ = resizeArray(abcDj_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
| 422 | abcSolution_ = resizeArray(abcSolution_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_); |
---|
| 423 | abcPerturbation_ = resizeArray(abcPerturbation_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
| 424 | offset_ = resizeArray(offset_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0); |
---|
| 425 | setupPointers(newSize1, newSize2); |
---|
[1878] | 426 | // Fill gaps in perturbation array |
---|
[2385] | 427 | fillPerturbation(maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_); |
---|
| 428 | fillPerturbation(newSize1 + maximumAbcNumberColumns_, newSize2 - maximumAbcNumberColumns_); |
---|
[1878] | 429 | // Clean gap |
---|
| 430 | //CoinIotaN(fromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_); |
---|
| 431 | //CoinIotaN(toExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_); |
---|
[2385] | 432 | CoinFillN(scaleFromExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0); |
---|
| 433 | CoinFillN(scaleToExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0); |
---|
| 434 | CoinFillN(internalStatusSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, static_cast< unsigned char >(1)); |
---|
| 435 | CoinFillN(lowerSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, -COIN_DBL_MAX); |
---|
| 436 | CoinFillN(upperSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, COIN_DBL_MAX); |
---|
| 437 | CoinFillN(costSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0); |
---|
| 438 | CoinFillN(djSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0); |
---|
| 439 | CoinFillN(solutionSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0); |
---|
| 440 | CoinFillN(offset_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0); |
---|
| 441 | maximumAbcNumberRows_ = newSize1; |
---|
| 442 | maximumAbcNumberColumns_ = newSize2; |
---|
| 443 | maximumNumberTotal_ = newSize1 + newSize2; |
---|
| 444 | delete[] abcPivotVariable_; |
---|
[1878] | 445 | abcPivotVariable_ = new int[maximumAbcNumberRows_]; |
---|
[2385] | 446 | for (int i = 0; i < ABC_NUMBER_USEFUL; i++) |
---|
| 447 | usefulArray_[i].reserve(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200)); |
---|
[1878] | 448 | } |
---|
[2385] | 449 | void AbcSimplex::translate(int type) |
---|
[1878] | 450 | { |
---|
| 451 | // clear bottom bits |
---|
[2385] | 452 | stateOfProblem_ &= ~(VALUES_PASS - 1); |
---|
| 453 | if ((type & DO_SCALE_AND_MATRIX) != 0) { |
---|
[1878] | 454 | //stateOfProblem_ |= DO_SCALE_AND_MATRIX; |
---|
| 455 | stateOfProblem_ |= DO_SCALE_AND_MATRIX; |
---|
| 456 | delete abcMatrix_; |
---|
[2385] | 457 | abcMatrix_ = new AbcMatrix(*matrix()); |
---|
[1878] | 458 | abcMatrix_->setModel(this); |
---|
| 459 | abcMatrix_->scale(scalingFlag_ ? 0 : -1); |
---|
| 460 | } |
---|
[2385] | 461 | if ((type & DO_STATUS) != 0 && (type & DO_BASIS_AND_ORDER) == 0) { |
---|
[1878] | 462 | // from Clp enum to Abc enum (and bound flip) |
---|
[2385] | 463 | unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 }; |
---|
| 464 | unsigned char *COIN_RESTRICT statusAbc = internalStatus_; |
---|
| 465 | const unsigned char *COIN_RESTRICT statusClp = status_; |
---|
[1878] | 466 | int i; |
---|
[2385] | 467 | for (i = numberRows_ - 1; i >= 0; i--) { |
---|
| 468 | unsigned char status = statusClp[i] & 7; |
---|
| 469 | bool basicClp = status == 1; |
---|
| 470 | bool basicAbc = (statusAbc[i] & 7) == 4; |
---|
| 471 | if (basicClp == basicAbc) |
---|
| 472 | statusAbc[i] = lookupToAbcSlack[status]; |
---|
[1878] | 473 | else |
---|
[2385] | 474 | break; |
---|
[1878] | 475 | } |
---|
| 476 | if (!i) { |
---|
| 477 | // from Clp enum to Abc enum |
---|
[2385] | 478 | unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 }; |
---|
| 479 | statusAbc += maximumAbcNumberRows_; |
---|
| 480 | statusClp += maximumAbcNumberRows_; |
---|
| 481 | for (i = numberColumns_ - 1; i >= 0; i--) { |
---|
| 482 | unsigned char status = statusClp[i] & 7; |
---|
| 483 | bool basicClp = status == 1; |
---|
| 484 | bool basicAbc = (statusAbc[i] & 7) == 4; |
---|
| 485 | if (basicClp == basicAbc) |
---|
| 486 | statusAbc[i] = lookupToAbc[status]; |
---|
| 487 | else |
---|
| 488 | break; |
---|
[1878] | 489 | } |
---|
| 490 | if (i) |
---|
[2385] | 491 | type |= DO_BASIS_AND_ORDER; |
---|
[1878] | 492 | } else { |
---|
[2385] | 493 | type |= DO_BASIS_AND_ORDER; |
---|
[1878] | 494 | } |
---|
| 495 | stateOfProblem_ |= DO_STATUS; |
---|
| 496 | } |
---|
[2385] | 497 | if ((type & DO_BASIS_AND_ORDER) != 0) { |
---|
[1878] | 498 | permuteIn(); |
---|
| 499 | permuteBasis(); |
---|
| 500 | stateOfProblem_ |= DO_BASIS_AND_ORDER; |
---|
| 501 | type &= ~DO_SOLUTION; |
---|
| 502 | } |
---|
[2385] | 503 | if ((type & DO_SOLUTION) != 0) { |
---|
[1878] | 504 | permuteBasis(); |
---|
| 505 | stateOfProblem_ |= DO_SOLUTION; |
---|
| 506 | } |
---|
[2385] | 507 | if ((type & DO_JUST_BOUNDS) != 0) { |
---|
[1878] | 508 | stateOfProblem_ |= DO_JUST_BOUNDS; |
---|
| 509 | } |
---|
| 510 | if (!type) { |
---|
| 511 | // just move stuff down |
---|
[2385] | 512 | CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_); |
---|
| 513 | CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_); |
---|
| 514 | CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_); |
---|
[1878] | 515 | } |
---|
| 516 | } |
---|
[2024] | 517 | #ifdef ABC_SPRINT |
---|
| 518 | // Overwrite to create sub problem (just internal arrays) - save full stuff |
---|
[2385] | 519 | AbcSimplex * |
---|
| 520 | AbcSimplex::createSubProblem(int numberColumns, const int *whichColumn) |
---|
[2024] | 521 | { |
---|
| 522 | int numberFullColumns = numberColumns_; |
---|
| 523 | numberColumns_ = numberColumns; |
---|
[2385] | 524 | AbcSimplex *subProblem = this; |
---|
| 525 | AbcSimplex *fullProblem = reinterpret_cast< AbcSimplex * >(new char[sizeof(AbcSimplex)]); |
---|
| 526 | memset(fullProblem, 0, sizeof(AbcSimplex)); |
---|
[2024] | 527 | fullProblem->maximumAbcNumberRows_ = maximumAbcNumberRows_; |
---|
| 528 | fullProblem->numberColumns_ = numberFullColumns; |
---|
| 529 | fullProblem->numberTotal_ = numberTotal_; |
---|
| 530 | fullProblem->maximumNumberTotal_ = maximumNumberTotal_; |
---|
| 531 | fullProblem->numberTotalWithoutFixed_ = numberTotalWithoutFixed_; |
---|
[2385] | 532 | fullProblem->abcPrimalColumnPivot_ = abcPrimalColumnPivot_; |
---|
| 533 | fullProblem->internalStatus_ = internalStatus_; |
---|
| 534 | fullProblem->abcLower_ = abcLower_; |
---|
| 535 | fullProblem->abcUpper_ = abcUpper_; |
---|
| 536 | fullProblem->abcCost_ = abcCost_; |
---|
| 537 | fullProblem->abcDj_ = abcDj_; |
---|
| 538 | fullProblem->abcSolution_ = abcSolution_; |
---|
| 539 | fullProblem->scaleFromExternal_ = scaleFromExternal_; |
---|
| 540 | fullProblem->offset_ = offset_; |
---|
| 541 | fullProblem->abcPerturbation_ = abcPerturbation_; |
---|
| 542 | fullProblem->abcPivotVariable_ = abcPivotVariable_; |
---|
| 543 | fullProblem->abcMatrix_ = abcMatrix_; |
---|
[2024] | 544 | fullProblem->abcNonLinearCost_ = abcNonLinearCost_; |
---|
[2385] | 545 | fullProblem->setupPointers(maximumAbcNumberRows_, numberFullColumns); |
---|
| 546 | subProblem->numberTotal_ = maximumAbcNumberRows_ + numberColumns; |
---|
| 547 | subProblem->maximumNumberTotal_ = maximumAbcNumberRows_ + numberColumns; |
---|
| 548 | subProblem->numberTotalWithoutFixed_ = subProblem->numberTotal_; |
---|
| 549 | int sizeArray = 2 * subProblem->maximumNumberTotal_ + maximumAbcNumberRows_; |
---|
[2429] | 550 | subProblem->internalStatus_ = newArray((unsigned char *)NULL, |
---|
[2385] | 551 | sizeArray + subProblem->maximumNumberTotal_); |
---|
[2429] | 552 | subProblem->abcLower_ = newArray((double *)NULL, sizeArray); |
---|
| 553 | subProblem->abcUpper_ = newArray((double *)NULL, sizeArray); |
---|
| 554 | subProblem->abcCost_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_); |
---|
| 555 | subProblem->abcDj_ = newArray((double *)NULL, sizeArray); |
---|
| 556 | subProblem->abcSolution_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_); |
---|
| 557 | //fromExternal_ = newArray((int *)NULL,sizeArray); |
---|
| 558 | //toExternal_ = newArray((int *)NULL,sizeArray); |
---|
| 559 | subProblem->scaleFromExternal_ = newArray((double *)NULL, sizeArray); |
---|
| 560 | subProblem->offset_ = newArray((double *)NULL, sizeArray); |
---|
| 561 | subProblem->abcPerturbation_ = newArray((double *)NULL, sizeArray); |
---|
| 562 | subProblem->abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_); |
---|
[2385] | 563 | subProblem->setupPointers(maximumAbcNumberRows_, numberColumns); |
---|
[2024] | 564 | // could use arrays - but for now be safe |
---|
[2385] | 565 | int *backward = new int[numberFullColumns + numberRows_]; |
---|
| 566 | int *whichRow = backward + numberFullColumns; |
---|
| 567 | for (int i = 0; i < numberFullColumns; i++) |
---|
| 568 | backward[i] = -1; |
---|
| 569 | for (int i = 0; i < numberColumns; i++) |
---|
| 570 | backward[whichColumn[i]] = i + numberRows_; |
---|
| 571 | for (int i = 0; i < numberRows_; i++) { |
---|
| 572 | whichRow[i] = i; |
---|
[2024] | 573 | int iPivot = fullProblem->abcPivotVariable_[i]; |
---|
[2385] | 574 | if (iPivot < numberRows_) { |
---|
| 575 | subProblem->abcPivotVariable_[i] = iPivot; |
---|
[2024] | 576 | } else { |
---|
[2385] | 577 | assert(backward[iPivot - numberRows_] >= 0); |
---|
| 578 | subProblem->abcPivotVariable_[i] = backward[iPivot - numberRows_]; |
---|
[2024] | 579 | } |
---|
| 580 | subProblem->internalStatus_[i] = fullProblem->internalStatus_[i]; |
---|
| 581 | subProblem->abcLower_[i] = fullProblem->abcLower_[i]; |
---|
| 582 | subProblem->abcUpper_[i] = fullProblem->abcUpper_[i]; |
---|
| 583 | subProblem->abcCost_[i] = fullProblem->abcCost_[i]; |
---|
| 584 | subProblem->abcDj_[i] = fullProblem->abcDj_[i]; |
---|
| 585 | subProblem->abcSolution_[i] = fullProblem->abcSolution_[i]; |
---|
| 586 | subProblem->abcPerturbation_[i] = fullProblem->abcPerturbation_[i]; |
---|
| 587 | subProblem->internalStatusSaved_[i] = fullProblem->internalStatusSaved_[i]; |
---|
| 588 | subProblem->perturbationSaved_[i] = fullProblem->perturbationSaved_[i]; |
---|
| 589 | subProblem->lowerSaved_[i] = fullProblem->lowerSaved_[i]; |
---|
| 590 | subProblem->upperSaved_[i] = fullProblem->upperSaved_[i]; |
---|
| 591 | subProblem->costSaved_[i] = fullProblem->costSaved_[i]; |
---|
| 592 | subProblem->djSaved_[i] = fullProblem->djSaved_[i]; |
---|
| 593 | subProblem->solutionSaved_[i] = fullProblem->solutionSaved_[i]; |
---|
| 594 | subProblem->offset_[i] = fullProblem->offset_[i]; |
---|
| 595 | subProblem->lowerBasic_[i] = fullProblem->lowerBasic_[i]; |
---|
| 596 | subProblem->upperBasic_[i] = fullProblem->upperBasic_[i]; |
---|
| 597 | subProblem->costBasic_[i] = fullProblem->costBasic_[i]; |
---|
| 598 | subProblem->solutionBasic_[i] = fullProblem->solutionBasic_[i]; |
---|
| 599 | subProblem->djBasic_[i] = fullProblem->djBasic_[i]; |
---|
| 600 | } |
---|
[2385] | 601 | for (int i = 0; i < numberColumns; i++) { |
---|
| 602 | int k = whichColumn[i]; |
---|
| 603 | subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k]; |
---|
| 604 | subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k]; |
---|
| 605 | subProblem->abcLower_[maximumAbcNumberRows_ + i] = fullProblem->abcLower_[maximumAbcNumberRows_ + k]; |
---|
| 606 | subProblem->abcUpper_[maximumAbcNumberRows_ + i] = fullProblem->abcUpper_[maximumAbcNumberRows_ + k]; |
---|
| 607 | subProblem->abcCost_[maximumAbcNumberRows_ + i] = fullProblem->abcCost_[maximumAbcNumberRows_ + k]; |
---|
| 608 | subProblem->abcDj_[maximumAbcNumberRows_ + i] = fullProblem->abcDj_[maximumAbcNumberRows_ + k]; |
---|
| 609 | subProblem->abcSolution_[maximumAbcNumberRows_ + i] = fullProblem->abcSolution_[maximumAbcNumberRows_ + k]; |
---|
| 610 | subProblem->abcPerturbation_[maximumAbcNumberRows_ + i] = fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k]; |
---|
| 611 | subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i] = fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k]; |
---|
| 612 | subProblem->perturbationSaved_[maximumAbcNumberRows_ + i] = fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k]; |
---|
| 613 | subProblem->lowerSaved_[maximumAbcNumberRows_ + i] = fullProblem->lowerSaved_[maximumAbcNumberRows_ + k]; |
---|
| 614 | subProblem->upperSaved_[maximumAbcNumberRows_ + i] = fullProblem->upperSaved_[maximumAbcNumberRows_ + k]; |
---|
| 615 | subProblem->costSaved_[maximumAbcNumberRows_ + i] = fullProblem->costSaved_[maximumAbcNumberRows_ + k]; |
---|
| 616 | subProblem->djSaved_[maximumAbcNumberRows_ + i] = fullProblem->djSaved_[maximumAbcNumberRows_ + k]; |
---|
| 617 | subProblem->solutionSaved_[maximumAbcNumberRows_ + i] = fullProblem->solutionSaved_[maximumAbcNumberRows_ + k]; |
---|
| 618 | subProblem->offset_[maximumAbcNumberRows_ + i] = fullProblem->offset_[maximumAbcNumberRows_ + k]; |
---|
[2024] | 619 | } |
---|
[2385] | 620 | subProblem->abcNonLinearCost_ = new AbcNonLinearCost(subProblem); |
---|
[2024] | 621 | subProblem->abcNonLinearCost_->checkInfeasibilities(0.0); |
---|
[2385] | 622 | subProblem->abcMatrix_ = new AbcMatrix(*fullProblem->abcMatrix_, numberRows_, whichRow, |
---|
| 623 | numberColumns, whichColumn); |
---|
[2024] | 624 | subProblem->abcMatrix_->setModel(subProblem); |
---|
| 625 | subProblem->abcMatrix_->rebalance(); |
---|
| 626 | subProblem->abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest(); |
---|
[2385] | 627 | subProblem->abcPrimalColumnPivot_->saveWeights(subProblem, 2); |
---|
| 628 | delete[] backward; |
---|
[2024] | 629 | // swap |
---|
| 630 | return fullProblem; |
---|
| 631 | } |
---|
| 632 | // Restore stuff from sub problem (and delete sub problem) |
---|
[2385] | 633 | void AbcSimplex::restoreFromSubProblem(AbcSimplex *fullProblem, const int *whichColumn) |
---|
[2024] | 634 | { |
---|
[2385] | 635 | AbcSimplex *subProblem = this; |
---|
| 636 | for (int i = 0; i < numberRows_; i++) { |
---|
[2024] | 637 | int iPivot = subProblem->abcPivotVariable_[i]; |
---|
[2385] | 638 | if (iPivot < numberRows_) { |
---|
| 639 | fullProblem->abcPivotVariable_[i] = iPivot; |
---|
[2024] | 640 | } else { |
---|
[2385] | 641 | fullProblem->abcPivotVariable_[i] = whichColumn[iPivot - numberRows_] + numberRows_; |
---|
[2024] | 642 | } |
---|
| 643 | fullProblem->internalStatus_[i] = subProblem->internalStatus_[i]; |
---|
| 644 | fullProblem->abcLower_[i] = subProblem->abcLower_[i]; |
---|
| 645 | fullProblem->abcUpper_[i] = subProblem->abcUpper_[i]; |
---|
| 646 | fullProblem->abcCost_[i] = subProblem->abcCost_[i]; |
---|
| 647 | fullProblem->abcDj_[i] = subProblem->abcDj_[i]; |
---|
| 648 | fullProblem->abcSolution_[i] = subProblem->abcSolution_[i]; |
---|
| 649 | fullProblem->abcPerturbation_[i] = subProblem->abcPerturbation_[i]; |
---|
| 650 | fullProblem->internalStatusSaved_[i] = subProblem->internalStatusSaved_[i]; |
---|
| 651 | fullProblem->perturbationSaved_[i] = subProblem->perturbationSaved_[i]; |
---|
| 652 | fullProblem->lowerSaved_[i] = subProblem->lowerSaved_[i]; |
---|
| 653 | fullProblem->upperSaved_[i] = subProblem->upperSaved_[i]; |
---|
| 654 | fullProblem->costSaved_[i] = subProblem->costSaved_[i]; |
---|
| 655 | fullProblem->djSaved_[i] = subProblem->djSaved_[i]; |
---|
| 656 | fullProblem->solutionSaved_[i] = subProblem->solutionSaved_[i]; |
---|
| 657 | fullProblem->offset_[i] = subProblem->offset_[i]; |
---|
| 658 | fullProblem->lowerBasic_[i] = subProblem->lowerBasic_[i]; |
---|
| 659 | fullProblem->upperBasic_[i] = subProblem->upperBasic_[i]; |
---|
| 660 | fullProblem->costBasic_[i] = subProblem->costBasic_[i]; |
---|
| 661 | fullProblem->solutionBasic_[i] = subProblem->solutionBasic_[i]; |
---|
| 662 | fullProblem->djBasic_[i] = subProblem->djBasic_[i]; |
---|
| 663 | } |
---|
| 664 | int numberColumns = subProblem->numberColumns_; |
---|
[2385] | 665 | for (int i = 0; i < numberColumns; i++) { |
---|
| 666 | int k = whichColumn[i]; |
---|
| 667 | fullProblem->internalStatus_[maximumAbcNumberRows_ + k] = subProblem->internalStatus_[maximumAbcNumberRows_ + i]; |
---|
| 668 | fullProblem->abcLower_[maximumAbcNumberRows_ + k] = subProblem->abcLower_[maximumAbcNumberRows_ + i]; |
---|
| 669 | fullProblem->abcUpper_[maximumAbcNumberRows_ + k] = subProblem->abcUpper_[maximumAbcNumberRows_ + i]; |
---|
| 670 | fullProblem->abcCost_[maximumAbcNumberRows_ + k] = subProblem->abcCost_[maximumAbcNumberRows_ + i]; |
---|
| 671 | fullProblem->abcDj_[maximumAbcNumberRows_ + k] = subProblem->abcDj_[maximumAbcNumberRows_ + i]; |
---|
| 672 | fullProblem->abcSolution_[maximumAbcNumberRows_ + k] = subProblem->abcSolution_[maximumAbcNumberRows_ + i]; |
---|
| 673 | fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k] = subProblem->abcPerturbation_[maximumAbcNumberRows_ + i]; |
---|
| 674 | fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k] = subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i]; |
---|
| 675 | fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k] = subProblem->perturbationSaved_[maximumAbcNumberRows_ + i]; |
---|
| 676 | fullProblem->lowerSaved_[maximumAbcNumberRows_ + k] = subProblem->lowerSaved_[maximumAbcNumberRows_ + i]; |
---|
| 677 | fullProblem->upperSaved_[maximumAbcNumberRows_ + k] = subProblem->upperSaved_[maximumAbcNumberRows_ + i]; |
---|
| 678 | fullProblem->costSaved_[maximumAbcNumberRows_ + k] = subProblem->costSaved_[maximumAbcNumberRows_ + i]; |
---|
| 679 | fullProblem->djSaved_[maximumAbcNumberRows_ + k] = subProblem->djSaved_[maximumAbcNumberRows_ + i]; |
---|
| 680 | fullProblem->solutionSaved_[maximumAbcNumberRows_ + k] = subProblem->solutionSaved_[maximumAbcNumberRows_ + i]; |
---|
| 681 | fullProblem->offset_[maximumAbcNumberRows_ + k] = subProblem->offset_[maximumAbcNumberRows_ + i]; |
---|
[2024] | 682 | } |
---|
[2385] | 683 | delete[] subProblem->internalStatus_; |
---|
| 684 | delete[] subProblem->abcPerturbation_; |
---|
[2024] | 685 | delete subProblem->abcMatrix_; |
---|
[2385] | 686 | delete[] subProblem->abcLower_; |
---|
| 687 | delete[] subProblem->abcUpper_; |
---|
| 688 | delete[] subProblem->abcCost_; |
---|
| 689 | delete[] subProblem->abcSolution_; |
---|
| 690 | delete[] subProblem->abcDj_; |
---|
[2024] | 691 | delete subProblem->abcPrimalColumnPivot_; |
---|
[2385] | 692 | delete[] subProblem->scaleFromExternal_; |
---|
| 693 | delete[] subProblem->offset_; |
---|
| 694 | delete[] subProblem->abcPivotVariable_; |
---|
| 695 | delete[] subProblem->reversePivotVariable_; |
---|
[2024] | 696 | delete subProblem->abcNonLinearCost_; |
---|
| 697 | numberColumns_ = fullProblem->numberColumns_; |
---|
[2385] | 698 | numberTotal_ = fullProblem->numberTotal_; |
---|
| 699 | maximumNumberTotal_ = fullProblem->maximumNumberTotal_; |
---|
| 700 | numberTotalWithoutFixed_ = fullProblem->numberTotalWithoutFixed_; |
---|
| 701 | abcPrimalColumnPivot_ = fullProblem->abcPrimalColumnPivot_; |
---|
| 702 | internalStatus_ = fullProblem->internalStatus_; |
---|
| 703 | abcLower_ = fullProblem->abcLower_; |
---|
| 704 | abcUpper_ = fullProblem->abcUpper_; |
---|
| 705 | abcCost_ = fullProblem->abcCost_; |
---|
| 706 | abcDj_ = fullProblem->abcDj_; |
---|
| 707 | abcSolution_ = fullProblem->abcSolution_; |
---|
| 708 | scaleFromExternal_ = fullProblem->scaleFromExternal_; |
---|
| 709 | offset_ = fullProblem->offset_; |
---|
| 710 | abcPerturbation_ = fullProblem->abcPerturbation_; |
---|
| 711 | abcPivotVariable_ = fullProblem->abcPivotVariable_; |
---|
| 712 | abcMatrix_ = fullProblem->abcMatrix_; |
---|
| 713 | setupPointers(maximumAbcNumberRows_, numberColumns); |
---|
[2024] | 714 | // ? redo nonlinearcost |
---|
| 715 | abcNonLinearCost_ = fullProblem->abcNonLinearCost_; |
---|
| 716 | abcNonLinearCost_->refresh(); |
---|
[2385] | 717 | delete[] reinterpret_cast< char * >(fullProblem); |
---|
[2024] | 718 | } |
---|
| 719 | #endif |
---|
[1878] | 720 | /* Sets dual values pass djs using unscaled duals |
---|
| 721 | type 1 - values pass |
---|
| 722 | type 2 - just use as infeasibility weights |
---|
| 723 | type 3 - as 2 but crash |
---|
| 724 | */ |
---|
[2385] | 725 | void AbcSimplex::setupDualValuesPass(const double *fakeDuals, |
---|
| 726 | const double *fakePrimals, |
---|
| 727 | int type) |
---|
[1878] | 728 | { |
---|
| 729 | // allslack basis |
---|
[2385] | 730 | memset(internalStatus_, 6, numberRows_); |
---|
[1878] | 731 | // temp |
---|
[2385] | 732 | if (type == 1) { |
---|
| 733 | bool allEqual = true; |
---|
| 734 | for (int i = 0; i < numberRows_; i++) { |
---|
| 735 | if (rowUpper_[i] > rowLower_[i]) { |
---|
| 736 | allEqual = false; |
---|
| 737 | break; |
---|
[1878] | 738 | } |
---|
| 739 | } |
---|
| 740 | if (allEqual) { |
---|
| 741 | // just modify costs |
---|
[2385] | 742 | transposeTimes(-1.0, fakeDuals, objective()); |
---|
[1878] | 743 | return; |
---|
| 744 | } |
---|
| 745 | } |
---|
| 746 | // compute unscaled djs |
---|
[2385] | 747 | CoinAbcMemcpy(djSaved_ + maximumAbcNumberRows_, objective(), numberColumns_); |
---|
| 748 | matrix_->transposeTimes(-1.0, fakeDuals, djSaved_ + maximumAbcNumberRows_); |
---|
[1878] | 749 | // save fake solution |
---|
[2385] | 750 | assert(solution_); |
---|
[1878] | 751 | //solution_ = new double[numberTotal_]; |
---|
[2385] | 752 | CoinAbcMemset0(solution_, numberRows_); |
---|
| 753 | CoinAbcMemcpy(solution_ + maximumAbcNumberRows_, fakePrimals, numberColumns_); |
---|
[1878] | 754 | // adjust |
---|
[2385] | 755 | for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) |
---|
| 756 | solution_[iSequence] -= offset_[iSequence]; |
---|
| 757 | matrix_->times(-1.0, solution_ + maximumAbcNumberRows_, solution_); |
---|
[1878] | 758 | double direction = optimizationDirection_; |
---|
[2385] | 759 | const double *COIN_RESTRICT rowScale = scaleFromExternal_; |
---|
| 760 | const double *COIN_RESTRICT inverseRowScale = scaleToExternal_; |
---|
| 761 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 762 | djSaved_[iRow] = direction * fakeDuals[iRow] * inverseRowScale[iRow]; |
---|
[1878] | 763 | solution_[iRow] *= rowScale[iRow]; |
---|
| 764 | } |
---|
[2385] | 765 | const double *COIN_RESTRICT columnScale = scaleToExternal_; |
---|
| 766 | for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++) |
---|
| 767 | djSaved_[iColumn] *= direction * columnScale[iColumn]; |
---|
[1878] | 768 | // Set solution values |
---|
[2385] | 769 | const double *lower = abcLower_ + maximumAbcNumberRows_; |
---|
| 770 | const double *upper = abcUpper_ + maximumAbcNumberRows_; |
---|
| 771 | double *solution = abcSolution_ + maximumAbcNumberRows_; |
---|
| 772 | const double *djs = djSaved_ + maximumAbcNumberRows_; |
---|
| 773 | const double *inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_; |
---|
| 774 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
| 775 | double thisValue = fakePrimals[iColumn]; |
---|
| 776 | double thisLower = columnLower_[iColumn]; |
---|
| 777 | double thisUpper = columnUpper_[iColumn]; |
---|
| 778 | double thisDj = djs[iColumn]; |
---|
| 779 | solution_[iColumn + maximumAbcNumberRows_] = solution_[iColumn + maximumAbcNumberRows_] * inverseColumnScale[iColumn]; |
---|
| 780 | if (thisLower > -1.0e30) { |
---|
| 781 | if (thisUpper < 1.0e30) { |
---|
| 782 | double gapUp = thisUpper - thisValue; |
---|
| 783 | double gapDown = thisValue - thisLower; |
---|
| 784 | bool goUp; |
---|
| 785 | if (gapUp > gapDown && thisDj > -dualTolerance_) { |
---|
| 786 | goUp = false; |
---|
| 787 | } else if (gapUp < gapDown && thisDj < dualTolerance_) { |
---|
| 788 | goUp = true; |
---|
| 789 | } else { |
---|
| 790 | // wild guess |
---|
| 791 | double badUp; |
---|
| 792 | double badDown; |
---|
| 793 | if (gapUp > gapDown) { |
---|
| 794 | badUp = gapUp * dualTolerance_; |
---|
| 795 | badDown = -gapDown * thisDj; |
---|
| 796 | } else { |
---|
| 797 | badUp = gapUp * thisDj; |
---|
| 798 | badDown = gapDown * dualTolerance_; |
---|
| 799 | } |
---|
| 800 | goUp = badDown > badUp; |
---|
| 801 | } |
---|
| 802 | if (goUp) { |
---|
| 803 | solution[iColumn] = upper[iColumn]; |
---|
| 804 | setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound); |
---|
| 805 | setStatus(iColumn, ClpSimplex::atUpperBound); |
---|
| 806 | } else { |
---|
| 807 | solution[iColumn] = lower[iColumn]; |
---|
| 808 | setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound); |
---|
| 809 | setStatus(iColumn, ClpSimplex::atLowerBound); |
---|
| 810 | } |
---|
[1878] | 811 | } else { |
---|
[2385] | 812 | solution[iColumn] = lower[iColumn]; |
---|
| 813 | setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound); |
---|
| 814 | setStatus(iColumn, ClpSimplex::atLowerBound); |
---|
[1878] | 815 | } |
---|
[2385] | 816 | } else if (thisUpper < 1.0e30) { |
---|
| 817 | solution[iColumn] = upper[iColumn]; |
---|
| 818 | setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound); |
---|
| 819 | setStatus(iColumn, ClpSimplex::atUpperBound); |
---|
[1878] | 820 | } else { |
---|
| 821 | // free |
---|
[2385] | 822 | solution[iColumn] = thisValue * inverseColumnScale[iColumn]; |
---|
| 823 | setInternalStatus(iColumn + maximumAbcNumberRows_, isFree); |
---|
| 824 | setStatus(iColumn, ClpSimplex::isFree); |
---|
[1878] | 825 | } |
---|
| 826 | } |
---|
| 827 | switch (type) { |
---|
| 828 | case 1: |
---|
| 829 | stateOfProblem_ |= VALUES_PASS; |
---|
| 830 | break; |
---|
| 831 | case 2: |
---|
| 832 | stateOfProblem_ |= VALUES_PASS2; |
---|
[2385] | 833 | delete[] solution_; |
---|
| 834 | solution_ = NULL; |
---|
[1878] | 835 | break; |
---|
| 836 | case 3: |
---|
| 837 | stateOfProblem_ |= VALUES_PASS2; |
---|
| 838 | break; |
---|
| 839 | } |
---|
| 840 | } |
---|
| 841 | //############################################################################# |
---|
[2385] | 842 | int AbcSimplex::computePrimals(CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector) |
---|
[1878] | 843 | { |
---|
[2385] | 844 | |
---|
[1878] | 845 | arrayVector->clear(); |
---|
| 846 | previousVector->clear(); |
---|
| 847 | // accumulate non basic stuff |
---|
[2385] | 848 | double *COIN_RESTRICT array = arrayVector->denseVector(); |
---|
| 849 | CoinAbcScatterZeroTo(abcSolution_, abcPivotVariable_, numberRows_); |
---|
[1878] | 850 | abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, array); |
---|
| 851 | #if 0 |
---|
| 852 | static int xxxxxx=0; |
---|
| 853 | if (xxxxxx==0) |
---|
| 854 | for (int i=0;i<numberRows_;i++) |
---|
| 855 | printf("%d %.18g\n",i,array[i]); |
---|
| 856 | #endif |
---|
| 857 | //arrayVector->scan(0,numberRows_,zeroTolerance_); |
---|
| 858 | // Ftran adjusted RHS and iterate to improve accuracy |
---|
| 859 | double lastError = COIN_DBL_MAX; |
---|
[2385] | 860 | CoinIndexedVector *thisVector = arrayVector; |
---|
| 861 | CoinIndexedVector *lastVector = previousVector; |
---|
[1878] | 862 | //if (arrayVector->getNumElements()) |
---|
| 863 | #if 0 |
---|
| 864 | double largest=0.0; |
---|
| 865 | int iLargest=-1; |
---|
| 866 | for (int i=0;i<numberRows_;i++) { |
---|
| 867 | if (fabs(array[i])>largest) { |
---|
| 868 | largest=array[i]; |
---|
| 869 | iLargest=i; |
---|
| 870 | } |
---|
| 871 | } |
---|
| 872 | printf("largest %g at row %d\n",array[iLargest],iLargest); |
---|
| 873 | #endif |
---|
| 874 | abcFactorization_->updateFullColumn(*thisVector); |
---|
| 875 | #if 0 |
---|
| 876 | largest=0.0; |
---|
| 877 | iLargest=-1; |
---|
| 878 | for (int i=0;i<numberRows_;i++) { |
---|
| 879 | if (fabs(array[i])>largest) { |
---|
| 880 | largest=array[i]; |
---|
| 881 | iLargest=i; |
---|
| 882 | } |
---|
| 883 | } |
---|
| 884 | printf("after largest %g at row %d\n",array[iLargest],iLargest); |
---|
| 885 | #endif |
---|
| 886 | #if 0 |
---|
| 887 | if (xxxxxx==0) |
---|
| 888 | for (int i=0;i<numberRows_;i++) |
---|
| 889 | printf("xx %d %.19g\n",i,array[i]); |
---|
| 890 | if (numberIterations_>300) |
---|
| 891 | exit(0); |
---|
| 892 | #endif |
---|
[2385] | 893 | int numberRefinements = 0; |
---|
[1878] | 894 | for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) { |
---|
| 895 | int numberIn = thisVector->getNumElements(); |
---|
[2385] | 896 | const int *COIN_RESTRICT indexIn = thisVector->getIndices(); |
---|
| 897 | const double *COIN_RESTRICT arrayIn = thisVector->denseVector(); |
---|
| 898 | CoinAbcScatterToList(arrayIn, abcSolution_, indexIn, abcPivotVariable_, numberIn); |
---|
[1878] | 899 | // check Ax == b (for all) |
---|
| 900 | abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, solutionBasic_); |
---|
| 901 | #if 0 |
---|
| 902 | if (xxxxxx==0) |
---|
| 903 | for (int i=0;i<numberTotal_;i++) |
---|
| 904 | printf("sol %d %.19g\n",i,abcSolution_[i]); |
---|
| 905 | if (xxxxxx==0) |
---|
| 906 | for (int i=0;i<numberRows_;i++) |
---|
| 907 | printf("basic %d %.19g\n",i,solutionBasic_[i]); |
---|
| 908 | #endif |
---|
| 909 | double multiplier = 131072.0; |
---|
[2385] | 910 | largestPrimalError_ = CoinAbcMaximumAbsElementAndScale(solutionBasic_, multiplier, numberRows_); |
---|
| 911 | double maxValue = 0.0; |
---|
| 912 | for (int i = 0; i < numberRows_; i++) { |
---|
| 913 | double value = fabs(solutionBasic_[i]); |
---|
| 914 | if (value > maxValue) { |
---|
[1878] | 915 | #if 0 |
---|
| 916 | if (xxxxxx==0) |
---|
| 917 | printf("larger %.19gg at row %d\n",value,i); |
---|
| 918 | maxValue=value; |
---|
| 919 | #endif |
---|
| 920 | } |
---|
| 921 | } |
---|
| 922 | if (largestPrimalError_ >= lastError) { |
---|
| 923 | // restore |
---|
[2385] | 924 | CoinIndexedVector *temp = thisVector; |
---|
[1878] | 925 | thisVector = lastVector; |
---|
| 926 | lastVector = temp; |
---|
| 927 | //goodSolution = false; |
---|
| 928 | break; |
---|
| 929 | } |
---|
| 930 | if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) { |
---|
| 931 | // try and make better |
---|
| 932 | numberRefinements++; |
---|
| 933 | // save this |
---|
[2385] | 934 | CoinIndexedVector *temp = thisVector; |
---|
[1878] | 935 | thisVector = lastVector; |
---|
| 936 | lastVector = temp; |
---|
[2385] | 937 | int *COIN_RESTRICT indexOut = thisVector->getIndices(); |
---|
[1878] | 938 | int number = 0; |
---|
| 939 | array = thisVector->denseVector(); |
---|
| 940 | thisVector->clear(); |
---|
| 941 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
[2385] | 942 | double value = solutionBasic_[iRow]; |
---|
| 943 | if (value) { |
---|
| 944 | array[iRow] = value; |
---|
| 945 | indexOut[number++] = iRow; |
---|
| 946 | solutionBasic_[iRow] = 0.0; |
---|
| 947 | } |
---|
[1878] | 948 | } |
---|
| 949 | thisVector->setNumElements(number); |
---|
| 950 | lastError = largestPrimalError_; |
---|
| 951 | abcFactorization_->updateFullColumn(*thisVector); |
---|
[2385] | 952 | double *previous = lastVector->denseVector(); |
---|
[1878] | 953 | number = 0; |
---|
[2385] | 954 | multiplier = 1.0 / multiplier; |
---|
[1878] | 955 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
[2385] | 956 | double value = previous[iRow] + multiplier * array[iRow]; |
---|
| 957 | if (value) { |
---|
| 958 | array[iRow] = value; |
---|
| 959 | indexOut[number++] = iRow; |
---|
| 960 | } else { |
---|
| 961 | array[iRow] = 0.0; |
---|
| 962 | } |
---|
[1878] | 963 | } |
---|
| 964 | thisVector->setNumElements(number); |
---|
| 965 | } else { |
---|
| 966 | break; |
---|
| 967 | } |
---|
| 968 | } |
---|
[2385] | 969 | |
---|
[1878] | 970 | // solution as accurate as we are going to get |
---|
| 971 | //if (!goodSolution) { |
---|
[2385] | 972 | CoinAbcMemcpy(solutionBasic_, thisVector->denseVector(), numberRows_); |
---|
| 973 | CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_); |
---|
[1878] | 974 | arrayVector->clear(); |
---|
| 975 | previousVector->clear(); |
---|
| 976 | return numberRefinements; |
---|
| 977 | } |
---|
| 978 | // Moves basic stuff to basic area |
---|
[2385] | 979 | void AbcSimplex::moveToBasic(int which) |
---|
[1878] | 980 | { |
---|
[2385] | 981 | if ((which & 1) != 0) |
---|
| 982 | CoinAbcGatherFrom(abcSolution_, solutionBasic_, abcPivotVariable_, numberRows_); |
---|
| 983 | if ((which & 2) != 0) |
---|
| 984 | CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_); |
---|
| 985 | if ((which & 4) != 0) |
---|
| 986 | CoinAbcGatherFrom(abcLower_, lowerBasic_, abcPivotVariable_, numberRows_); |
---|
| 987 | if ((which & 8) != 0) |
---|
| 988 | CoinAbcGatherFrom(abcUpper_, upperBasic_, abcPivotVariable_, numberRows_); |
---|
[1878] | 989 | } |
---|
| 990 | // now dual side |
---|
[2385] | 991 | int AbcSimplex::computeDuals(double *givenDjs, CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector) |
---|
[1878] | 992 | { |
---|
[2385] | 993 | double *COIN_RESTRICT array = arrayVector->denseVector(); |
---|
| 994 | int *COIN_RESTRICT index = arrayVector->getIndices(); |
---|
[1878] | 995 | int number = 0; |
---|
| 996 | if (!givenDjs) { |
---|
| 997 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 998 | double value = costBasic_[iRow]; |
---|
| 999 | if (value) { |
---|
[2385] | 1000 | array[iRow] = value; |
---|
| 1001 | index[number++] = iRow; |
---|
[1878] | 1002 | } |
---|
| 1003 | } |
---|
| 1004 | } else { |
---|
| 1005 | // dual values pass - djs may not be zero |
---|
| 1006 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 1007 | int iPivot = abcPivotVariable_[iRow]; |
---|
| 1008 | // make sure zero if done |
---|
| 1009 | if (!pivoted(iPivot)) |
---|
[2385] | 1010 | givenDjs[iPivot] = 0.0; |
---|
[1878] | 1011 | double value = abcCost_[iPivot] - givenDjs[iPivot]; |
---|
| 1012 | if (value) { |
---|
[2385] | 1013 | array[iRow] = value; |
---|
| 1014 | index[number++] = iRow; |
---|
[1878] | 1015 | } |
---|
| 1016 | } |
---|
| 1017 | } |
---|
| 1018 | arrayVector->setNumElements(number); |
---|
| 1019 | // Btran basic costs and get as accurate as possible |
---|
| 1020 | double lastError = COIN_DBL_MAX; |
---|
[2385] | 1021 | CoinIndexedVector *thisVector = arrayVector; |
---|
| 1022 | CoinIndexedVector *lastVector = previousVector; |
---|
[1878] | 1023 | abcFactorization_->updateFullColumnTranspose(*thisVector); |
---|
[2385] | 1024 | |
---|
| 1025 | int numberRefinements = 0; |
---|
| 1026 | for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) { |
---|
[1878] | 1027 | // check basic reduced costs zero |
---|
| 1028 | // put reduced cost of basic into djBasic_ |
---|
[2385] | 1029 | CoinAbcMemcpy(djBasic_, costBasic_, numberRows_); |
---|
| 1030 | abcMatrix_->transposeTimesBasic(-1.0, thisVector->denseVector(), djBasic_); |
---|
| 1031 | largestDualError_ = CoinAbcMaximumAbsElement(djBasic_, numberRows_); |
---|
[1878] | 1032 | if (largestDualError_ >= lastError) { |
---|
| 1033 | // restore |
---|
[2385] | 1034 | CoinIndexedVector *temp = thisVector; |
---|
[1878] | 1035 | thisVector = lastVector; |
---|
| 1036 | lastVector = temp; |
---|
| 1037 | break; |
---|
| 1038 | } |
---|
| 1039 | if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10 |
---|
[2385] | 1040 | && !givenDjs) { |
---|
[1878] | 1041 | numberRefinements++; |
---|
| 1042 | // try and make better |
---|
| 1043 | // save this |
---|
[2385] | 1044 | CoinIndexedVector *temp = thisVector; |
---|
[1878] | 1045 | thisVector = lastVector; |
---|
| 1046 | lastVector = temp; |
---|
| 1047 | array = thisVector->denseVector(); |
---|
| 1048 | double multiplier = 131072.0; |
---|
| 1049 | //array=djBasic_*multiplier |
---|
[2385] | 1050 | CoinAbcMultiplyAdd(djBasic_, numberRows_, multiplier, array, 0.0); |
---|
[1878] | 1051 | lastError = largestDualError_; |
---|
[2385] | 1052 | abcFactorization_->updateFullColumnTranspose(*thisVector); |
---|
[1878] | 1053 | #if ABC_DEBUG |
---|
| 1054 | thisVector->checkClean(); |
---|
| 1055 | #endif |
---|
| 1056 | multiplier = 1.0 / multiplier; |
---|
[2385] | 1057 | double *COIN_RESTRICT previous = lastVector->denseVector(); |
---|
[1878] | 1058 | // array = array*multiplier+previous |
---|
[2385] | 1059 | CoinAbcMultiplyAdd(previous, numberRows_, 1.0, array, multiplier); |
---|
[1878] | 1060 | } else { |
---|
| 1061 | break; |
---|
| 1062 | } |
---|
| 1063 | } |
---|
| 1064 | // now look at dual solution |
---|
[2385] | 1065 | CoinAbcMemcpy(dual_, thisVector->denseVector(), numberRows_); |
---|
| 1066 | CoinAbcMemset0(thisVector->denseVector(), numberRows_); |
---|
[1878] | 1067 | thisVector->setNumElements(0); |
---|
| 1068 | if (numberRefinements) { |
---|
[2385] | 1069 | CoinAbcMemset0(lastVector->denseVector(), numberRows_); |
---|
[1878] | 1070 | lastVector->setNumElements(0); |
---|
| 1071 | } |
---|
[2385] | 1072 | CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_); |
---|
| 1073 | abcMatrix_->transposeTimesNonBasic(-1.0, dual_, abcDj_); |
---|
[1878] | 1074 | // If necessary - override results |
---|
| 1075 | if (givenDjs) { |
---|
| 1076 | // restore accurate duals |
---|
| 1077 | CoinMemcpyN(abcDj_, numberTotal_, givenDjs); |
---|
| 1078 | } |
---|
| 1079 | //arrayVector->clear(); |
---|
| 1080 | //previousVector->clear(); |
---|
| 1081 | #if ABC_DEBUG |
---|
| 1082 | arrayVector->checkClean(); |
---|
| 1083 | previousVector->checkClean(); |
---|
| 1084 | #endif |
---|
| 1085 | return numberRefinements; |
---|
| 1086 | } |
---|
| 1087 | |
---|
| 1088 | /* Factorizes using current basis. |
---|
| 1089 | solveType - 1 iterating, 0 initial, -1 external |
---|
| 1090 | - 2 then iterating but can throw out of basis |
---|
| 1091 | If 10 added then in primal values pass |
---|
| 1092 | Return codes are as from AbcSimplexFactorization unless initial factorization |
---|
| 1093 | when total number of singularities is returned. |
---|
| 1094 | Special case is numberRows_+1 -> all slack basis. |
---|
| 1095 | */ |
---|
[2385] | 1096 | int AbcSimplex::internalFactorize(int solveType) |
---|
[1878] | 1097 | { |
---|
[2385] | 1098 | assert(status_); |
---|
| 1099 | |
---|
[1878] | 1100 | bool valuesPass = false; |
---|
| 1101 | if (solveType >= 10) { |
---|
| 1102 | valuesPass = true; |
---|
| 1103 | solveType -= 10; |
---|
| 1104 | } |
---|
| 1105 | #if 0 |
---|
| 1106 | // Make sure everything is clean |
---|
| 1107 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
| 1108 | if(getInternalStatus(iSequence) == isFixed) { |
---|
| 1109 | // double check fixed |
---|
| 1110 | assert (abcUpper_[iSequence] == abcLower_[iSequence]); |
---|
| 1111 | assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<100.0*primalTolerance_); |
---|
| 1112 | } else if (getInternalStatus(iSequence) == isFree) { |
---|
| 1113 | assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX); |
---|
| 1114 | } else if (getInternalStatus(iSequence) == atLowerBound) { |
---|
| 1115 | assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_); |
---|
| 1116 | } else if (getInternalStatus(iSequence) == atUpperBound) { |
---|
| 1117 | assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_); |
---|
| 1118 | } else if (getInternalStatus(iSequence) == superBasic) { |
---|
| 1119 | assert (!valuesPass); |
---|
| 1120 | } |
---|
| 1121 | } |
---|
| 1122 | #endif |
---|
| 1123 | #if 0 //ndef NDEBUG |
---|
| 1124 | // Make sure everything is clean |
---|
| 1125 | double sumOutside=0.0; |
---|
| 1126 | int numberOutside=0; |
---|
| 1127 | //double sumOutsideLarge=0.0; |
---|
| 1128 | int numberOutsideLarge=0; |
---|
| 1129 | double sumInside=0.0; |
---|
| 1130 | int numberInside=0; |
---|
| 1131 | //double sumInsideLarge=0.0; |
---|
| 1132 | int numberInsideLarge=0; |
---|
| 1133 | char rowcol[] = {'R', 'C'}; |
---|
| 1134 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
| 1135 | if(getInternalStatus(iSequence) == isFixed) { |
---|
| 1136 | // double check fixed |
---|
| 1137 | assert (abcUpper_[iSequence] == abcLower_[iSequence]); |
---|
| 1138 | assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<primalTolerance_); |
---|
| 1139 | } else if (getInternalStatus(iSequence) == isFree) { |
---|
| 1140 | assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX); |
---|
| 1141 | } else if (getInternalStatus(iSequence) == atLowerBound) { |
---|
| 1142 | assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_); |
---|
| 1143 | if (abcSolution_[iSequence]<abcLower_[iSequence]) { |
---|
| 1144 | numberOutside++; |
---|
[2385] | 1145 | #if ABC_NORMAL_DEBUG > 1 |
---|
[1878] | 1146 | if (handler_->logLevel()==63) |
---|
| 1147 | printf("%c%d below by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence), |
---|
| 1148 | abcLower_[iSequence]-abcSolution_[iSequence]); |
---|
| 1149 | #endif |
---|
| 1150 | sumOutside-=abcSolution_[iSequence]-abcLower_[iSequence]; |
---|
| 1151 | if (abcSolution_[iSequence]<abcLower_[iSequence]-primalTolerance_) |
---|
| 1152 | numberOutsideLarge++; |
---|
| 1153 | } else if (abcSolution_[iSequence]>abcLower_[iSequence]) { |
---|
| 1154 | numberInside++; |
---|
| 1155 | sumInside+=abcSolution_[iSequence]-abcLower_[iSequence]; |
---|
| 1156 | if (abcSolution_[iSequence]>abcLower_[iSequence]+primalTolerance_) |
---|
| 1157 | numberInsideLarge++; |
---|
| 1158 | } |
---|
| 1159 | } else if (getInternalStatus(iSequence) == atUpperBound) { |
---|
| 1160 | assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_); |
---|
| 1161 | if (abcSolution_[iSequence]>abcUpper_[iSequence]) { |
---|
| 1162 | numberOutside++; |
---|
[2385] | 1163 | #if ABC_NORMAL_DEBUG > 0 |
---|
[1878] | 1164 | if (handler_->logLevel()==63) |
---|
| 1165 | printf("%c%d above by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence), |
---|
| 1166 | -(abcUpper_[iSequence]-abcSolution_[iSequence])); |
---|
| 1167 | #endif |
---|
| 1168 | sumOutside+=abcSolution_[iSequence]-abcUpper_[iSequence]; |
---|
| 1169 | if (abcSolution_[iSequence]>abcUpper_[iSequence]+primalTolerance_) |
---|
| 1170 | numberOutsideLarge++; |
---|
| 1171 | } else if (abcSolution_[iSequence]<abcUpper_[iSequence]) { |
---|
| 1172 | numberInside++; |
---|
| 1173 | sumInside-=abcSolution_[iSequence]-abcUpper_[iSequence]; |
---|
| 1174 | if (abcSolution_[iSequence]<abcUpper_[iSequence]-primalTolerance_) |
---|
| 1175 | numberInsideLarge++; |
---|
| 1176 | } |
---|
| 1177 | } else if (getInternalStatus(iSequence) == superBasic) { |
---|
| 1178 | //assert (!valuesPass); |
---|
| 1179 | } |
---|
| 1180 | } |
---|
[2385] | 1181 | #if ABC_NORMAL_DEBUG > 0 |
---|
[1878] | 1182 | if (numberInside+numberOutside) |
---|
| 1183 | printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n", |
---|
| 1184 | numberOutside,sumOutside,numberOutsideLarge, |
---|
| 1185 | numberInside,sumInside,numberInsideLarge); |
---|
| 1186 | #endif |
---|
| 1187 | #endif |
---|
[2024] | 1188 | // *** replace below by cleanStatus |
---|
[1878] | 1189 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[2385] | 1190 | AbcSimplex::Status status = getInternalStatus(iSequence); |
---|
| 1191 | if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence]) |
---|
| 1192 | setInternalStatus(iSequence, isFixed); |
---|
[1878] | 1193 | } |
---|
[2385] | 1194 | if (numberIterations_ == baseIteration_ && !valuesPass) { |
---|
[1878] | 1195 | double largeValue = this->largeValue(); |
---|
[2385] | 1196 | double *COIN_RESTRICT solution = abcSolution_; |
---|
[1878] | 1197 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[2385] | 1198 | AbcSimplex::Status status = getInternalStatus(iSequence); |
---|
| 1199 | if (status == superBasic) { |
---|
| 1200 | double lower = abcLower_[iSequence]; |
---|
| 1201 | double upper = abcUpper_[iSequence]; |
---|
| 1202 | double value = solution[iSequence]; |
---|
| 1203 | AbcSimplex::Status thisStatus = isFree; |
---|
| 1204 | if (lower > -largeValue || upper < largeValue) { |
---|
| 1205 | if (lower != upper) { |
---|
| 1206 | if (fabs(value - lower) < fabs(value - upper)) { |
---|
| 1207 | thisStatus = AbcSimplex::atLowerBound; |
---|
| 1208 | solution[iSequence] = lower; |
---|
| 1209 | } else { |
---|
| 1210 | thisStatus = AbcSimplex::atUpperBound; |
---|
| 1211 | solution[iSequence] = upper; |
---|
| 1212 | } |
---|
| 1213 | } else { |
---|
| 1214 | thisStatus = AbcSimplex::isFixed; |
---|
| 1215 | solution[iSequence] = upper; |
---|
| 1216 | } |
---|
| 1217 | setInternalStatus(iSequence, thisStatus); |
---|
| 1218 | } |
---|
[1878] | 1219 | } |
---|
| 1220 | } |
---|
| 1221 | } |
---|
| 1222 | int status = abcFactorization_->factorize(this, solveType, valuesPass); |
---|
| 1223 | if (status) { |
---|
| 1224 | handler_->message(CLP_SIMPLEX_BADFACTOR, messages_) |
---|
| 1225 | << status |
---|
| 1226 | << CoinMessageEol; |
---|
| 1227 | return -1; |
---|
| 1228 | } else if (!solveType) { |
---|
| 1229 | // Initial basis - return number of singularities |
---|
| 1230 | int numberSlacks = 0; |
---|
| 1231 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 1232 | if (getInternalStatus(iRow) == basic) |
---|
[2385] | 1233 | numberSlacks++; |
---|
[1878] | 1234 | } |
---|
| 1235 | status = CoinMax(numberSlacks - numberRows_, 0); |
---|
[1990] | 1236 | if (status) |
---|
[2385] | 1237 | printf("%d singularities\n", status); |
---|
[1878] | 1238 | // special case if all slack |
---|
| 1239 | if (numberSlacks == numberRows_) { |
---|
| 1240 | status = numberRows_ + 1; |
---|
| 1241 | } |
---|
| 1242 | } |
---|
[2078] | 1243 | #ifdef ABC_USE_COIN_FACTORIZATION |
---|
| 1244 | // sparse methods |
---|
| 1245 | abcFactorization_->goSparse(); |
---|
| 1246 | #endif |
---|
[1878] | 1247 | return status; |
---|
| 1248 | } |
---|
[2024] | 1249 | // Make sure no superbasic etc |
---|
[2385] | 1250 | void AbcSimplex::cleanStatus(bool valuesPass) |
---|
[2024] | 1251 | { |
---|
| 1252 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[2385] | 1253 | AbcSimplex::Status status = getInternalStatus(iSequence); |
---|
| 1254 | if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence]) |
---|
| 1255 | setInternalStatus(iSequence, isFixed); |
---|
[2024] | 1256 | } |
---|
[2385] | 1257 | if (numberIterations_ == baseIteration_ && !valuesPass) { |
---|
[2024] | 1258 | double largeValue = this->largeValue(); |
---|
[2385] | 1259 | double *COIN_RESTRICT solution = abcSolution_; |
---|
[2024] | 1260 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[2385] | 1261 | AbcSimplex::Status status = getInternalStatus(iSequence); |
---|
| 1262 | if (status == superBasic) { |
---|
| 1263 | double lower = abcLower_[iSequence]; |
---|
| 1264 | double upper = abcUpper_[iSequence]; |
---|
| 1265 | double value = solution[iSequence]; |
---|
| 1266 | AbcSimplex::Status thisStatus = isFree; |
---|
| 1267 | if (lower > -largeValue || upper < largeValue) { |
---|
| 1268 | if (lower != upper) { |
---|
| 1269 | if (fabs(value - lower) < fabs(value - upper)) { |
---|
| 1270 | thisStatus = AbcSimplex::atLowerBound; |
---|
| 1271 | solution[iSequence] = lower; |
---|
| 1272 | } else { |
---|
| 1273 | thisStatus = AbcSimplex::atUpperBound; |
---|
| 1274 | solution[iSequence] = upper; |
---|
| 1275 | } |
---|
| 1276 | } else { |
---|
| 1277 | thisStatus = AbcSimplex::isFixed; |
---|
| 1278 | solution[iSequence] = upper; |
---|
| 1279 | } |
---|
| 1280 | setInternalStatus(iSequence, thisStatus); |
---|
| 1281 | } |
---|
[2024] | 1282 | } |
---|
| 1283 | } |
---|
| 1284 | } |
---|
| 1285 | } |
---|
[1878] | 1286 | // Sets objectiveValue_ from rawObjectiveValue_ |
---|
[2385] | 1287 | void AbcSimplex::setClpSimplexObjectiveValue() |
---|
[1878] | 1288 | { |
---|
[2385] | 1289 | objectiveValue_ = rawObjectiveValue_ / (objectiveScale_ * rhsScale_); |
---|
[1878] | 1290 | objectiveValue_ += objectiveOffset_; |
---|
| 1291 | } |
---|
| 1292 | /* |
---|
| 1293 | This does basis housekeeping and does values for in/out variables. |
---|
| 1294 | Can also decide to re-factorize |
---|
| 1295 | */ |
---|
[2385] | 1296 | int AbcSimplex::housekeeping() |
---|
[1878] | 1297 | { |
---|
| 1298 | #define DELAYED_UPDATE |
---|
| 1299 | #ifdef DELAYED_UPDATE |
---|
[2385] | 1300 | if (algorithm_ < 0) { |
---|
[1878] | 1301 | // modify dualout |
---|
| 1302 | dualOut_ /= alpha_; |
---|
| 1303 | dualOut_ *= -directionOut_; |
---|
| 1304 | //double oldValue = valueIn_; |
---|
| 1305 | if (directionIn_ == -1) { |
---|
| 1306 | // as if from upper bound |
---|
| 1307 | valueIn_ = upperIn_ + dualOut_; |
---|
| 1308 | #if 0 //def ABC_DEBUG |
---|
| 1309 | printf("In from upper dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n", |
---|
| 1310 | dualOut_,movement_,movementOld,valueIn_,upperIn_+movement_); |
---|
| 1311 | #endif |
---|
| 1312 | } else { |
---|
| 1313 | // as if from lower bound |
---|
| 1314 | valueIn_ = lowerIn_ + dualOut_; |
---|
| 1315 | #if 0 //def ABC_DEBUG |
---|
| 1316 | printf("In from lower dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n", |
---|
| 1317 | dualOut_,movement_,movementOld,valueIn_,lowerIn_+movement_); |
---|
| 1318 | #endif |
---|
| 1319 | } |
---|
| 1320 | // outgoing |
---|
| 1321 | if (directionOut_ > 0) { |
---|
| 1322 | valueOut_ = lowerOut_; |
---|
| 1323 | } else { |
---|
| 1324 | valueOut_ = upperOut_; |
---|
| 1325 | } |
---|
[2385] | 1326 | #if ABC_NORMAL_DEBUG > 3 |
---|
[1878] | 1327 | if (rawObjectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16)) |
---|
| 1328 | printf("obj backwards %g %g\n", rawObjectiveValue_, oldobj); |
---|
| 1329 | #endif |
---|
| 1330 | } |
---|
| 1331 | #endif |
---|
| 1332 | #if 0 //ndef NDEBUG |
---|
| 1333 | { |
---|
| 1334 | int numberFlagged=0; |
---|
| 1335 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 1336 | int iPivot = abcPivotVariable_[iRow]; |
---|
| 1337 | if (flagged(iPivot)) |
---|
| 1338 | numberFlagged++; |
---|
| 1339 | } |
---|
| 1340 | assert (numberFlagged==numberFlagged_); |
---|
| 1341 | } |
---|
| 1342 | #endif |
---|
| 1343 | // save value of incoming and outgoing |
---|
| 1344 | numberIterations_++; |
---|
| 1345 | changeMade_++; // something has happened |
---|
| 1346 | // incoming variable |
---|
| 1347 | if (handler_->logLevel() > 7) { |
---|
| 1348 | //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) { |
---|
| 1349 | handler_->message(CLP_SIMPLEX_HOUSE1, messages_) |
---|
| 1350 | << directionOut_ |
---|
| 1351 | << directionIn_ << theta_ |
---|
| 1352 | << dualOut_ << dualIn_ << alpha_ |
---|
| 1353 | << CoinMessageEol; |
---|
| 1354 | if (getInternalStatus(sequenceIn_) == isFree) { |
---|
| 1355 | handler_->message(CLP_SIMPLEX_FREEIN, messages_) |
---|
[2385] | 1356 | << sequenceIn_ |
---|
| 1357 | << CoinMessageEol; |
---|
[1878] | 1358 | } |
---|
| 1359 | } |
---|
| 1360 | // change of incoming |
---|
[2385] | 1361 | char rowcol[] = { 'R', 'C' }; |
---|
[1878] | 1362 | if (abcUpper_[sequenceIn_] > 1.0e20 && abcLower_[sequenceIn_] < -1.0e20) |
---|
| 1363 | progressFlag_ |= 2; // making real progress |
---|
| 1364 | #ifndef DELAYED_UPDATE |
---|
| 1365 | abcSolution_[sequenceIn_] = valueIn_; |
---|
| 1366 | #endif |
---|
| 1367 | if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] < 1.0e-12) |
---|
| 1368 | progressFlag_ |= 1; // making real progress |
---|
| 1369 | if (sequenceIn_ != sequenceOut_) { |
---|
| 1370 | if (alphaAccuracy_ > 0.0) { |
---|
| 1371 | double value = fabs(alpha_); |
---|
| 1372 | if (value > 1.0) |
---|
[2385] | 1373 | alphaAccuracy_ *= value; |
---|
[1878] | 1374 | else |
---|
[2385] | 1375 | alphaAccuracy_ /= value; |
---|
[1878] | 1376 | } |
---|
| 1377 | setInternalStatus(sequenceIn_, basic); |
---|
| 1378 | if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] > 0) { |
---|
| 1379 | // As Nonlinear costs may have moved bounds (to more feasible) |
---|
| 1380 | // Redo using value |
---|
| 1381 | if (fabs(valueOut_ - abcLower_[sequenceOut_]) < fabs(valueOut_ - abcUpper_[sequenceOut_])) { |
---|
[2385] | 1382 | // going to lower |
---|
| 1383 | setInternalStatus(sequenceOut_, atLowerBound); |
---|
[1878] | 1384 | } else { |
---|
[2385] | 1385 | // going to upper |
---|
| 1386 | setInternalStatus(sequenceOut_, atUpperBound); |
---|
[1878] | 1387 | } |
---|
| 1388 | } else { |
---|
| 1389 | // fixed |
---|
| 1390 | setInternalStatus(sequenceOut_, isFixed); |
---|
| 1391 | } |
---|
| 1392 | #ifndef DELAYED_UPDATE |
---|
| 1393 | abcSolution_[sequenceOut_] = valueOut_; |
---|
| 1394 | #endif |
---|
[2385] | 1395 | #if PARTITION_ROW_COPY == 1 |
---|
[1878] | 1396 | // move in row copy |
---|
[2385] | 1397 | abcMatrix_->inOutUseful(sequenceIn_, sequenceOut_); |
---|
[1878] | 1398 | #endif |
---|
| 1399 | } else { |
---|
| 1400 | //if (objective_->type()<2) |
---|
| 1401 | //assert (fabs(theta_)>1.0e-13); |
---|
| 1402 | // flip from bound to bound |
---|
| 1403 | // As Nonlinear costs may have moved bounds (to more feasible) |
---|
| 1404 | // Redo using value |
---|
| 1405 | if (fabs(valueIn_ - abcLower_[sequenceIn_]) < fabs(valueIn_ - abcUpper_[sequenceIn_])) { |
---|
| 1406 | // as if from upper bound |
---|
| 1407 | setInternalStatus(sequenceIn_, atLowerBound); |
---|
| 1408 | } else { |
---|
| 1409 | // as if from lower bound |
---|
| 1410 | setInternalStatus(sequenceIn_, atUpperBound); |
---|
| 1411 | } |
---|
| 1412 | } |
---|
| 1413 | setClpSimplexObjectiveValue(); |
---|
| 1414 | if (handler_->logLevel() > 7) { |
---|
| 1415 | //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) { |
---|
| 1416 | handler_->message(CLP_SIMPLEX_HOUSE2, messages_) |
---|
| 1417 | << numberIterations_ << objectiveValue() |
---|
| 1418 | << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_) |
---|
| 1419 | << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_); |
---|
| 1420 | handler_->printing(algorithm_ < 0) << dualOut_ << theta_; |
---|
| 1421 | handler_->printing(algorithm_ > 0) << dualIn_ << theta_; |
---|
| 1422 | handler_->message() << CoinMessageEol; |
---|
| 1423 | } |
---|
| 1424 | #if 0 |
---|
| 1425 | if (numberIterations_ > 10000) |
---|
| 1426 | printf(" it %d %g %c%d %c%d\n" |
---|
| 1427 | , numberIterations_, objectiveValue() |
---|
| 1428 | , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_) |
---|
| 1429 | , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_)); |
---|
| 1430 | #endif |
---|
| 1431 | if (hitMaximumIterations()) |
---|
| 1432 | return 2; |
---|
| 1433 | // check for small cycles |
---|
| 1434 | int in = sequenceIn_; |
---|
| 1435 | int out = sequenceOut_; |
---|
| 1436 | int cycle = abcProgress_.cycle(in, out, |
---|
[2385] | 1437 | directionIn_, directionOut_); |
---|
| 1438 | if (cycle > 0) { |
---|
| 1439 | #if ABC_NORMAL_DEBUG > 0 |
---|
[1878] | 1440 | if (handler_->logLevel() >= 63) |
---|
| 1441 | printf("Cycle of %d\n", cycle); |
---|
| 1442 | #endif |
---|
| 1443 | // reset |
---|
| 1444 | abcProgress_.startCheck(); |
---|
| 1445 | double random = randomNumberGenerator_.randomDouble(); |
---|
[2385] | 1446 | int extra = static_cast< int >(9.999 * random); |
---|
| 1447 | int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 }; |
---|
[1878] | 1448 | if (abcFactorization_->pivots() > cycle) { |
---|
| 1449 | forceFactorization_ = CoinMax(1, cycle - off[extra]); |
---|
| 1450 | } else { |
---|
| 1451 | // need to reject something |
---|
| 1452 | int iSequence; |
---|
[2385] | 1453 | if (algorithm_ < 0) { |
---|
| 1454 | iSequence = sequenceIn_; |
---|
[1878] | 1455 | } else { |
---|
[2385] | 1456 | /* should be better if don't reject incoming |
---|
[1878] | 1457 | as it is in basis */ |
---|
[2385] | 1458 | iSequence = sequenceOut_; |
---|
| 1459 | // so can't happen immediately again |
---|
| 1460 | sequenceOut_ = -1; |
---|
[1878] | 1461 | } |
---|
| 1462 | char x = isColumn(iSequence) ? 'C' : 'R'; |
---|
| 1463 | if (handler_->logLevel() >= 63) |
---|
[2385] | 1464 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
| 1465 | << x << sequenceWithin(iSequence) |
---|
| 1466 | << CoinMessageEol; |
---|
[1878] | 1467 | setFlagged(iSequence); |
---|
| 1468 | //printf("flagging %d\n",iSequence); |
---|
| 1469 | } |
---|
| 1470 | return 1; |
---|
| 1471 | } |
---|
[2385] | 1472 | int invertNow = 0; |
---|
[1878] | 1473 | // only time to re-factorize if one before real time |
---|
| 1474 | // this is so user won't be surprised that maximumPivots has exact meaning |
---|
| 1475 | int numberPivots = abcFactorization_->pivots(); |
---|
[2385] | 1476 | if (algorithm_ < 0) |
---|
[1878] | 1477 | numberPivots++; // allow for update not done |
---|
| 1478 | int maximumPivots = abcFactorization_->maximumPivots(); |
---|
| 1479 | int numberDense = 0; //abcFactorization_->numberDense(); |
---|
[2385] | 1480 | bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations()); |
---|
| 1481 | if (numberPivots == maximumPivots || maximumPivots < 2) { |
---|
[1878] | 1482 | // If dense then increase |
---|
[2074] | 1483 | if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots && false) { |
---|
[1878] | 1484 | abcFactorization_->maximumPivots(numberDense); |
---|
| 1485 | } |
---|
| 1486 | //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots); |
---|
[2385] | 1487 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
[2078] | 1488 | abcFactorization_->statsRefactor('M'); |
---|
| 1489 | #endif |
---|
[1878] | 1490 | return 1; |
---|
| 1491 | } else if ((abcFactorization_->timeToRefactorize() && !dontInvert) |
---|
[2385] | 1492 | || invertNow) { |
---|
[1878] | 1493 | //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots); |
---|
[2385] | 1494 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
[2078] | 1495 | abcFactorization_->statsRefactor('T'); |
---|
| 1496 | #endif |
---|
[1878] | 1497 | return 1; |
---|
| 1498 | } else if (forceFactorization_ > 0 && |
---|
| 1499 | #ifndef DELAYED_UPDATE |
---|
[2385] | 1500 | abcFactorization_->pivots() |
---|
[1878] | 1501 | #else |
---|
[2385] | 1502 | abcFactorization_->pivots() + 1 |
---|
| 1503 | #endif |
---|
| 1504 | >= forceFactorization_) { |
---|
[1878] | 1505 | // relax |
---|
| 1506 | forceFactorization_ = (3 + 5 * forceFactorization_) / 4; |
---|
| 1507 | if (forceFactorization_ > abcFactorization_->maximumPivots()) |
---|
| 1508 | forceFactorization_ = -1; //off |
---|
[2385] | 1509 | //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots); |
---|
| 1510 | #if CLP_FACTORIZATION_NEW_TIMING > 1 |
---|
[2078] | 1511 | abcFactorization_->statsRefactor('F'); |
---|
| 1512 | #endif |
---|
[1878] | 1513 | return 1; |
---|
| 1514 | } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) { |
---|
| 1515 | // A bit worried problem may be cycling - lets factorize at random intervals for a short period |
---|
| 1516 | int numberTooManyIterations = numberIterations_ - 10 * (numberRows_ + (numberColumns_ >> 2)); |
---|
| 1517 | double random = randomNumberGenerator_.randomDouble(); |
---|
[2385] | 1518 | int window = numberTooManyIterations % 5000; |
---|
| 1519 | if (window < 2 * maximumPivots) |
---|
| 1520 | random = 0.2 * random + 0.8; // randomly re-factorize but not too soon |
---|
[1878] | 1521 | else |
---|
[2385] | 1522 | random = 1.0; // switch off if not in window of opportunity |
---|
[1878] | 1523 | int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots); |
---|
| 1524 | if (abcFactorization_->pivots() >= random * maxNumber) { |
---|
| 1525 | //printf("ZZ cycling a %d\n",numberPivots); |
---|
| 1526 | return 1; |
---|
[2385] | 1527 | } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1000010 + 10 * (numberRows_ + (numberColumns_ >> 2))) { |
---|
[1878] | 1528 | //printf("ZZ cycling b %d\n",numberPivots); |
---|
| 1529 | return 1; |
---|
| 1530 | } else { |
---|
| 1531 | // carry on iterating |
---|
| 1532 | return 0; |
---|
| 1533 | } |
---|
| 1534 | } else { |
---|
| 1535 | // carry on iterating |
---|
| 1536 | return 0; |
---|
| 1537 | } |
---|
| 1538 | } |
---|
[2385] | 1539 | // Swaps primal stuff |
---|
| 1540 | void AbcSimplex::swapPrimalStuff() |
---|
[1878] | 1541 | { |
---|
[2385] | 1542 | if (sequenceOut_ < 0) |
---|
[1878] | 1543 | return; |
---|
[2385] | 1544 | assert(sequenceIn_ >= 0); |
---|
[1878] | 1545 | abcSolution_[sequenceOut_] = valueOut_; |
---|
| 1546 | abcSolution_[sequenceIn_] = valueIn_; |
---|
[2385] | 1547 | solutionBasic_[pivotRow_] = valueIn_; |
---|
| 1548 | if (algorithm_ < 0) { |
---|
[1878] | 1549 | // and set bounds correctly |
---|
[2385] | 1550 | static_cast< AbcSimplexDual * >(this)->originalBound(sequenceIn_); |
---|
| 1551 | static_cast< AbcSimplexDual * >(this)->changeBound(sequenceOut_); |
---|
[1878] | 1552 | } else { |
---|
[2385] | 1553 | #if ABC_NORMAL_DEBUG > 2 |
---|
| 1554 | if (handler_->logLevel() == 63) |
---|
[1878] | 1555 | printf("Code swapPrimalStuff for primal\n"); |
---|
| 1556 | #endif |
---|
| 1557 | } |
---|
[2385] | 1558 | if (pivotRow_ >= 0) { // may be flip in primal |
---|
| 1559 | lowerBasic_[pivotRow_] = abcLower_[sequenceIn_]; |
---|
| 1560 | upperBasic_[pivotRow_] = abcUpper_[sequenceIn_]; |
---|
| 1561 | costBasic_[pivotRow_] = abcCost_[sequenceIn_]; |
---|
[1878] | 1562 | abcPivotVariable_[pivotRow_] = sequenceIn_; |
---|
| 1563 | } |
---|
| 1564 | } |
---|
| 1565 | // Swaps dual stuff |
---|
[2385] | 1566 | void AbcSimplex::swapDualStuff(int lastSequenceOut, int lastDirectionOut) |
---|
[1878] | 1567 | { |
---|
| 1568 | // outgoing |
---|
| 1569 | // set dj to zero unless values pass |
---|
[2385] | 1570 | if (lastSequenceOut >= 0) { |
---|
| 1571 | if ((stateOfProblem_ & VALUES_PASS) == 0) { |
---|
[1878] | 1572 | if (lastDirectionOut > 0) { |
---|
[2385] | 1573 | abcDj_[lastSequenceOut] = theta_; |
---|
[1878] | 1574 | } else { |
---|
[2385] | 1575 | abcDj_[lastSequenceOut] = -theta_; |
---|
[1878] | 1576 | } |
---|
| 1577 | } else { |
---|
| 1578 | if (lastDirectionOut > 0) { |
---|
[2385] | 1579 | abcDj_[lastSequenceOut] -= theta_; |
---|
[1878] | 1580 | } else { |
---|
[2385] | 1581 | abcDj_[lastSequenceOut] += theta_; |
---|
[1878] | 1582 | } |
---|
| 1583 | } |
---|
| 1584 | } |
---|
[2385] | 1585 | if (sequenceIn_ >= 0) { |
---|
| 1586 | abcDj_[sequenceIn_] = 0.0; |
---|
[1878] | 1587 | //costBasic_[pivotRow_]=abcCost_[sequenceIn_]; |
---|
| 1588 | } |
---|
| 1589 | } |
---|
[2385] | 1590 | static void solveMany(int number, ClpSimplex **simplex) |
---|
[1878] | 1591 | { |
---|
[2385] | 1592 | for (int i = 0; i < number - 1; i++) |
---|
[1878] | 1593 | cilk_spawn simplex[i]->dual(0); |
---|
[2385] | 1594 | simplex[number - 1]->dual(0); |
---|
[1878] | 1595 | cilk_sync; |
---|
| 1596 | } |
---|
[2385] | 1597 | void AbcSimplex::crash(int type) |
---|
[1878] | 1598 | { |
---|
| 1599 | int i; |
---|
[2385] | 1600 | for (i = 0; i < numberRows_; i++) { |
---|
| 1601 | if (getInternalStatus(i) != basic) |
---|
[1878] | 1602 | break; |
---|
| 1603 | } |
---|
[2385] | 1604 | if (i < numberRows_) |
---|
[1878] | 1605 | return; |
---|
| 1606 | // all slack |
---|
[2385] | 1607 | if (type == 3) { |
---|
[1878] | 1608 | // decomposition crash |
---|
| 1609 | if (!abcMatrix_->gotRowCopy()) |
---|
| 1610 | abcMatrix_->createRowCopy(); |
---|
| 1611 | //const double * element = abcMatrix_->getPackedMatrix()->getElements(); |
---|
[2385] | 1612 | const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts(); |
---|
| 1613 | const int *row = abcMatrix_->getPackedMatrix()->getIndices(); |
---|
[1878] | 1614 | //const double * elementByRow = abcMatrix_->rowElements(); |
---|
[2385] | 1615 | const CoinBigIndex *rowStart = abcMatrix_->rowStart(); |
---|
| 1616 | const CoinBigIndex *rowEnd = abcMatrix_->rowEnd(); |
---|
| 1617 | const int *column = abcMatrix_->rowColumns(); |
---|
| 1618 | int *blockStart = usefulArray_[0].getIndices(); |
---|
| 1619 | int *columnBlock = blockStart + numberRows_; |
---|
| 1620 | int *nextColumn = usefulArray_[1].getIndices(); |
---|
| 1621 | int *blockCount = nextColumn + numberColumns_; |
---|
| 1622 | int direction[2] = { -1, 1 }; |
---|
| 1623 | int bestBreak = -1; |
---|
| 1624 | double bestValue = 0.0; |
---|
| 1625 | int iPass = 0; |
---|
| 1626 | int halfway = (numberRows_ + 1) / 2; |
---|
| 1627 | int firstMaster = -1; |
---|
| 1628 | int lastMaster = -2; |
---|
| 1629 | int numberBlocks = 0; |
---|
| 1630 | int largestRows = 0; |
---|
| 1631 | int largestColumns = 0; |
---|
| 1632 | int numberEmpty = 0; |
---|
| 1633 | int numberMaster = 0; |
---|
| 1634 | int numberEmptyColumns = 0; |
---|
| 1635 | int numberMasterColumns = 0; |
---|
| 1636 | while (iPass < 2) { |
---|
| 1637 | int increment = direction[iPass]; |
---|
| 1638 | int start = increment > 0 ? 0 : numberRows_ - 1; |
---|
| 1639 | int stop = increment > 0 ? numberRows_ : -1; |
---|
| 1640 | numberBlocks = 0; |
---|
| 1641 | int thisBestBreak = -1; |
---|
| 1642 | double thisBestValue = COIN_DBL_MAX; |
---|
| 1643 | int numberRowsDone = 0; |
---|
| 1644 | int numberMarkedColumns = 0; |
---|
| 1645 | int maximumBlockSize = 0; |
---|
| 1646 | for (int i = 0; i < numberRows_; i++) { |
---|
| 1647 | blockStart[i] = -1; |
---|
| 1648 | blockCount[i] = 0; |
---|
[1878] | 1649 | } |
---|
[2385] | 1650 | for (int i = 0; i < numberColumns_; i++) { |
---|
| 1651 | columnBlock[i] = -1; |
---|
| 1652 | nextColumn[i] = -1; |
---|
[1878] | 1653 | } |
---|
[2385] | 1654 | for (int iRow = start; iRow != stop; iRow += increment) { |
---|
| 1655 | int iBlock = -1; |
---|
| 1656 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 1657 | int iColumn = column[j] - maximumAbcNumberRows_; |
---|
| 1658 | int whichColumnBlock = columnBlock[iColumn]; |
---|
| 1659 | if (whichColumnBlock >= 0) { |
---|
| 1660 | // column marked |
---|
| 1661 | if (iBlock < 0) { |
---|
| 1662 | // put row in that block |
---|
| 1663 | iBlock = whichColumnBlock; |
---|
| 1664 | } else if (iBlock != whichColumnBlock) { |
---|
| 1665 | // merge |
---|
| 1666 | blockCount[iBlock] += blockCount[whichColumnBlock]; |
---|
| 1667 | blockCount[whichColumnBlock] = 0; |
---|
| 1668 | int jColumn = blockStart[whichColumnBlock]; |
---|
[1878] | 1669 | #ifndef NDEBUG |
---|
[2385] | 1670 | int iiColumn = iColumn; |
---|
[1878] | 1671 | #endif |
---|
[2385] | 1672 | while (jColumn >= 0) { |
---|
| 1673 | assert(columnBlock[jColumn] == whichColumnBlock); |
---|
| 1674 | columnBlock[jColumn] = iBlock; |
---|
[1878] | 1675 | #ifndef NDEBUG |
---|
[2385] | 1676 | if (jColumn == iiColumn) |
---|
| 1677 | iiColumn = -1; |
---|
[1878] | 1678 | #endif |
---|
[2385] | 1679 | iColumn = jColumn; |
---|
| 1680 | jColumn = nextColumn[jColumn]; |
---|
| 1681 | } |
---|
[1878] | 1682 | #ifndef NDEBUG |
---|
[2385] | 1683 | assert(iiColumn == -1); |
---|
[1878] | 1684 | #endif |
---|
[2385] | 1685 | nextColumn[iColumn] = blockStart[iBlock]; |
---|
| 1686 | blockStart[iBlock] = blockStart[whichColumnBlock]; |
---|
| 1687 | blockStart[whichColumnBlock] = -1; |
---|
| 1688 | } |
---|
| 1689 | } |
---|
| 1690 | } |
---|
| 1691 | int n = numberMarkedColumns; |
---|
| 1692 | if (iBlock < 0) { |
---|
| 1693 | //new block |
---|
| 1694 | if (rowEnd[iRow] > rowStart[iRow]) { |
---|
| 1695 | numberBlocks++; |
---|
| 1696 | iBlock = numberBlocks; |
---|
| 1697 | int jColumn = column[rowStart[iRow]] - maximumAbcNumberRows_; |
---|
| 1698 | columnBlock[jColumn] = iBlock; |
---|
| 1699 | blockStart[iBlock] = jColumn; |
---|
| 1700 | numberMarkedColumns++; |
---|
| 1701 | for (CoinBigIndex j = rowStart[iRow] + 1; j < rowEnd[iRow]; j++) { |
---|
| 1702 | int iColumn = column[j] - maximumAbcNumberRows_; |
---|
| 1703 | columnBlock[iColumn] = iBlock; |
---|
| 1704 | numberMarkedColumns++; |
---|
| 1705 | nextColumn[jColumn] = iColumn; |
---|
| 1706 | jColumn = iColumn; |
---|
| 1707 | } |
---|
| 1708 | blockCount[iBlock] = numberMarkedColumns - n; |
---|
| 1709 | } else { |
---|
| 1710 | // empty |
---|
| 1711 | } |
---|
| 1712 | } else { |
---|
| 1713 | // put in existing block |
---|
| 1714 | int jColumn = blockStart[iBlock]; |
---|
| 1715 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 1716 | int iColumn = column[j] - maximumAbcNumberRows_; |
---|
| 1717 | assert(columnBlock[iColumn] < 0 || columnBlock[iColumn] == iBlock); |
---|
| 1718 | if (columnBlock[iColumn] < 0) { |
---|
| 1719 | columnBlock[iColumn] = iBlock; |
---|
| 1720 | numberMarkedColumns++; |
---|
| 1721 | nextColumn[iColumn] = jColumn; |
---|
| 1722 | jColumn = iColumn; |
---|
| 1723 | } |
---|
| 1724 | } |
---|
| 1725 | blockStart[iBlock] = jColumn; |
---|
| 1726 | blockCount[iBlock] += numberMarkedColumns - n; |
---|
| 1727 | } |
---|
| 1728 | if (iBlock >= 0) |
---|
| 1729 | maximumBlockSize = CoinMax(maximumBlockSize, blockCount[iBlock]); |
---|
| 1730 | numberRowsDone++; |
---|
| 1731 | if (thisBestValue * numberRowsDone > maximumBlockSize && numberRowsDone > halfway) { |
---|
| 1732 | thisBestBreak = iRow; |
---|
| 1733 | thisBestValue = static_cast< double >(maximumBlockSize) / static_cast< double >(numberRowsDone); |
---|
| 1734 | } |
---|
[1878] | 1735 | } |
---|
[2385] | 1736 | if (thisBestBreak == stop) |
---|
| 1737 | thisBestValue = COIN_DBL_MAX; |
---|
[1878] | 1738 | iPass++; |
---|
[2385] | 1739 | if (iPass == 1) { |
---|
| 1740 | bestBreak = thisBestBreak; |
---|
| 1741 | bestValue = thisBestValue; |
---|
[1878] | 1742 | } else { |
---|
[2385] | 1743 | if (bestValue < thisBestValue) { |
---|
| 1744 | firstMaster = 0; |
---|
| 1745 | lastMaster = bestBreak; |
---|
| 1746 | } else { |
---|
| 1747 | firstMaster = thisBestBreak; // ? +1 |
---|
| 1748 | lastMaster = numberRows_; |
---|
| 1749 | } |
---|
[1878] | 1750 | } |
---|
| 1751 | } |
---|
[2385] | 1752 | bool useful = false; |
---|
| 1753 | if (firstMaster < lastMaster) { |
---|
| 1754 | for (int i = 0; i < numberRows_; i++) |
---|
| 1755 | blockStart[i] = -1; |
---|
| 1756 | for (int i = firstMaster; i < lastMaster; i++) |
---|
| 1757 | blockStart[i] = -2; |
---|
| 1758 | for (int i = 0; i < numberColumns_; i++) |
---|
| 1759 | columnBlock[i] = -1; |
---|
| 1760 | int firstRow = 0; |
---|
| 1761 | numberBlocks = -1; |
---|
[1878] | 1762 | while (true) { |
---|
[2385] | 1763 | for (; firstRow < numberRows_; firstRow++) { |
---|
| 1764 | if (blockStart[firstRow] == -1) |
---|
| 1765 | break; |
---|
| 1766 | } |
---|
| 1767 | if (firstRow == numberRows_) |
---|
| 1768 | break; |
---|
| 1769 | int nRows = 0; |
---|
| 1770 | numberBlocks++; |
---|
| 1771 | int numberStack = 1; |
---|
| 1772 | blockCount[0] = firstRow; |
---|
| 1773 | while (numberStack) { |
---|
| 1774 | int iRow = blockCount[--numberStack]; |
---|
| 1775 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 1776 | int iColumn = column[j] - maximumAbcNumberRows_; |
---|
| 1777 | int iBlock = columnBlock[iColumn]; |
---|
| 1778 | if (iBlock < 0) { |
---|
| 1779 | columnBlock[iColumn] = numberBlocks; |
---|
| 1780 | for (CoinBigIndex k = columnStart[iColumn]; |
---|
| 1781 | k < columnStart[iColumn + 1]; k++) { |
---|
| 1782 | int jRow = row[k]; |
---|
| 1783 | int rowBlock = blockStart[jRow]; |
---|
| 1784 | if (rowBlock == -1) { |
---|
| 1785 | nRows++; |
---|
| 1786 | blockStart[jRow] = numberBlocks; |
---|
| 1787 | blockCount[numberStack++] = jRow; |
---|
| 1788 | } |
---|
| 1789 | } |
---|
| 1790 | } |
---|
| 1791 | } |
---|
| 1792 | } |
---|
| 1793 | if (!nRows) { |
---|
| 1794 | // empty!! |
---|
| 1795 | numberBlocks--; |
---|
| 1796 | } |
---|
| 1797 | firstRow++; |
---|
[1878] | 1798 | } |
---|
[2385] | 1799 | // adjust |
---|
[1878] | 1800 | numberBlocks++; |
---|
[2385] | 1801 | for (int i = 0; i < numberBlocks; i++) { |
---|
| 1802 | blockCount[i] = 0; |
---|
| 1803 | nextColumn[i] = 0; |
---|
[1878] | 1804 | } |
---|
| 1805 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
[2385] | 1806 | int iBlock = blockStart[iRow]; |
---|
| 1807 | if (iBlock >= 0) { |
---|
| 1808 | blockCount[iBlock]++; |
---|
| 1809 | } else { |
---|
| 1810 | if (iBlock == -2) |
---|
| 1811 | numberMaster++; |
---|
| 1812 | else |
---|
| 1813 | numberEmpty++; |
---|
| 1814 | } |
---|
[1878] | 1815 | } |
---|
| 1816 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
[2385] | 1817 | int iBlock = columnBlock[iColumn]; |
---|
| 1818 | if (iBlock >= 0) { |
---|
| 1819 | nextColumn[iBlock]++; |
---|
| 1820 | } else { |
---|
| 1821 | if (columnStart[iColumn + 1] > columnStart[iColumn]) |
---|
| 1822 | numberMasterColumns++; |
---|
| 1823 | else |
---|
| 1824 | numberEmptyColumns++; |
---|
| 1825 | } |
---|
[1878] | 1826 | } |
---|
[2385] | 1827 | for (int i = 0; i < numberBlocks; i++) { |
---|
| 1828 | if (blockCount[i] + nextColumn[i] > largestRows + largestColumns) { |
---|
| 1829 | largestRows = blockCount[i]; |
---|
| 1830 | largestColumns = nextColumn[i]; |
---|
| 1831 | } |
---|
[1878] | 1832 | } |
---|
[2385] | 1833 | useful = true; |
---|
| 1834 | if (numberMaster > halfway || largestRows * 3 > numberRows_) |
---|
| 1835 | useful = false; |
---|
[1878] | 1836 | } |
---|
| 1837 | if (useful) { |
---|
[2385] | 1838 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 1839 | if (logLevel() >= 2) |
---|
| 1840 | printf("%d master rows %d <= < %d\n", lastMaster - firstMaster, |
---|
| 1841 | firstMaster, lastMaster); |
---|
[1878] | 1842 | printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n", |
---|
[2385] | 1843 | useful ? "**Useful" : "NoGood", |
---|
| 1844 | numberBlocks, largestRows, largestColumns, numberMaster, numberEmpty, numberRows_, |
---|
| 1845 | numberMasterColumns, numberEmptyColumns, numberColumns_); |
---|
| 1846 | if (logLevel() >= 2) { |
---|
| 1847 | for (int i = 0; i < numberBlocks; i++) |
---|
| 1848 | printf("Block %d has %d rows and %d columns\n", |
---|
| 1849 | i, blockCount[i], nextColumn[i]); |
---|
[1878] | 1850 | } |
---|
| 1851 | #endif |
---|
| 1852 | #define NUMBER_DW_BLOCKS 20 |
---|
[2385] | 1853 | int minSize1 = (numberRows_ - numberMaster + NUMBER_DW_BLOCKS - 1) / NUMBER_DW_BLOCKS; |
---|
| 1854 | int minSize2 = (numberRows_ - numberMaster + 2 * NUMBER_DW_BLOCKS - 1) / (2 * NUMBER_DW_BLOCKS); |
---|
| 1855 | int *backRow = usefulArray_[1].getIndices(); |
---|
[1878] | 1856 | // first sort |
---|
[2385] | 1857 | for (int i = 0; i < numberBlocks; i++) { |
---|
| 1858 | backRow[i] = -(3 * blockCount[i] + 0 * nextColumn[i]); |
---|
| 1859 | nextColumn[i] = i; |
---|
[1878] | 1860 | } |
---|
[2385] | 1861 | CoinSort_2(backRow, backRow + numberBlocks, nextColumn); |
---|
[1878] | 1862 | // keep if >minSize2 or sum >minSize1; |
---|
[2385] | 1863 | int n = 0; |
---|
| 1864 | for (int i = 0; i < numberBlocks; i++) { |
---|
| 1865 | int originalBlock = nextColumn[i]; |
---|
| 1866 | if (blockCount[originalBlock] < minSize2) |
---|
| 1867 | break; |
---|
| 1868 | n++; |
---|
[1878] | 1869 | } |
---|
[2385] | 1870 | int size = minSize1; |
---|
| 1871 | for (int i = n; i < numberBlocks; i++) { |
---|
| 1872 | int originalBlock = nextColumn[i]; |
---|
| 1873 | size -= blockCount[originalBlock]; |
---|
| 1874 | if (size > 0 && i < numberBlocks - 1) { |
---|
| 1875 | blockCount[originalBlock] = -1; |
---|
| 1876 | } else { |
---|
| 1877 | size = minSize1; |
---|
| 1878 | n++; |
---|
| 1879 | } |
---|
[1878] | 1880 | } |
---|
[2385] | 1881 | int n2 = numberBlocks; |
---|
| 1882 | numberBlocks = n; |
---|
| 1883 | for (int i = n2 - 1; i >= 0; i--) { |
---|
| 1884 | int originalBlock = nextColumn[i]; |
---|
| 1885 | if (blockCount[originalBlock] > 0) |
---|
| 1886 | n--; |
---|
| 1887 | blockCount[originalBlock] = n; |
---|
[1878] | 1888 | } |
---|
[2385] | 1889 | assert(!n); |
---|
[1878] | 1890 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
[2385] | 1891 | int iBlock = blockStart[iRow]; |
---|
| 1892 | if (iBlock >= 0) |
---|
| 1893 | blockStart[iRow] = blockCount[iBlock]; |
---|
[1878] | 1894 | } |
---|
| 1895 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
[2385] | 1896 | int iBlock = columnBlock[iColumn]; |
---|
| 1897 | if (iBlock >= 0) |
---|
| 1898 | columnBlock[iColumn] = blockCount[iBlock]; |
---|
[1878] | 1899 | } |
---|
| 1900 | // stick to Clp for now |
---|
[2385] | 1901 | ClpSimplex **simplex = new ClpSimplex *[numberBlocks]; |
---|
| 1902 | for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { |
---|
| 1903 | int nRow = 0; |
---|
| 1904 | int nColumn = 0; |
---|
| 1905 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 1906 | if (blockStart[iRow] == iBlock) |
---|
| 1907 | blockCount[nRow++] = iRow; |
---|
| 1908 | } |
---|
| 1909 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
| 1910 | if (columnBlock[iColumn] == iBlock) |
---|
| 1911 | nextColumn[nColumn++] = iColumn; |
---|
| 1912 | } |
---|
| 1913 | simplex[iBlock] = new ClpSimplex(this, nRow, blockCount, nColumn, nextColumn); |
---|
| 1914 | simplex[iBlock]->setSpecialOptions(simplex[iBlock]->specialOptions() & (~65536)); |
---|
| 1915 | if (logLevel() < 2) |
---|
| 1916 | simplex[iBlock]->setLogLevel(0); |
---|
[1878] | 1917 | } |
---|
[2385] | 1918 | solveMany(numberBlocks, simplex); |
---|
| 1919 | int numberBasic = numberMaster; |
---|
| 1920 | int numberStructurals = 0; |
---|
| 1921 | for (int i = 0; i < numberMaster; i++) |
---|
| 1922 | abcPivotVariable_[i] = i + firstMaster; |
---|
| 1923 | for (int iBlock = 0; iBlock < numberBlocks; iBlock++) { |
---|
| 1924 | int nRow = 0; |
---|
| 1925 | int nColumn = 0; |
---|
| 1926 | // from Clp enum to Abc enum (and bound flip) |
---|
| 1927 | unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 }; |
---|
| 1928 | unsigned char *COIN_RESTRICT getStatus = simplex[iBlock]->statusArray() + simplex[iBlock]->numberColumns(); |
---|
| 1929 | double *COIN_RESTRICT solutionSaved = solutionSaved_; |
---|
| 1930 | double *COIN_RESTRICT lowerSaved = lowerSaved_; |
---|
| 1931 | double *COIN_RESTRICT upperSaved = upperSaved_; |
---|
| 1932 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 1933 | if (blockStart[iRow] == iBlock) { |
---|
| 1934 | unsigned char status = getStatus[nRow++] & 7; |
---|
| 1935 | AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbcSlack[status]); |
---|
| 1936 | if (status != ClpSimplex::basic) { |
---|
| 1937 | double lowerValue = lowerSaved[iRow]; |
---|
| 1938 | double upperValue = upperSaved[iRow]; |
---|
| 1939 | if (lowerValue == -COIN_DBL_MAX) { |
---|
| 1940 | if (upperValue == COIN_DBL_MAX) { |
---|
| 1941 | // free |
---|
| 1942 | abcStatus = isFree; |
---|
| 1943 | } else { |
---|
| 1944 | abcStatus = atUpperBound; |
---|
| 1945 | } |
---|
| 1946 | } else if (upperValue == COIN_DBL_MAX) { |
---|
| 1947 | abcStatus = atLowerBound; |
---|
| 1948 | } else if (lowerValue == upperValue) { |
---|
| 1949 | abcStatus = isFixed; |
---|
| 1950 | } |
---|
| 1951 | switch (abcStatus) { |
---|
| 1952 | case isFixed: |
---|
| 1953 | case atLowerBound: |
---|
| 1954 | solutionSaved[iRow] = lowerValue; |
---|
| 1955 | break; |
---|
| 1956 | case atUpperBound: |
---|
| 1957 | solutionSaved[iRow] = upperValue; |
---|
| 1958 | break; |
---|
| 1959 | default: |
---|
| 1960 | break; |
---|
| 1961 | } |
---|
| 1962 | } else { |
---|
| 1963 | // basic |
---|
| 1964 | abcPivotVariable_[numberBasic++] = iRow; |
---|
| 1965 | } |
---|
| 1966 | internalStatus_[iRow] = abcStatus; |
---|
| 1967 | } |
---|
| 1968 | } |
---|
| 1969 | // from Clp enum to Abc enum |
---|
| 1970 | unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 }; |
---|
| 1971 | unsigned char *COIN_RESTRICT putStatus = internalStatus_ + maximumAbcNumberRows_; |
---|
| 1972 | getStatus = simplex[iBlock]->statusArray(); |
---|
| 1973 | solutionSaved += maximumAbcNumberRows_; |
---|
| 1974 | lowerSaved += maximumAbcNumberRows_; |
---|
| 1975 | upperSaved += maximumAbcNumberRows_; |
---|
| 1976 | int numberSaved = numberBasic; |
---|
| 1977 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
| 1978 | if (columnBlock[iColumn] == iBlock) { |
---|
| 1979 | unsigned char status = getStatus[nColumn++] & 7; |
---|
| 1980 | AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbc[status]); |
---|
| 1981 | if (status != ClpSimplex::basic) { |
---|
| 1982 | double lowerValue = lowerSaved[iColumn]; |
---|
| 1983 | double upperValue = upperSaved[iColumn]; |
---|
| 1984 | if (lowerValue == -COIN_DBL_MAX) { |
---|
| 1985 | if (upperValue == COIN_DBL_MAX) { |
---|
| 1986 | // free |
---|
| 1987 | abcStatus = isFree; |
---|
| 1988 | } else { |
---|
| 1989 | abcStatus = atUpperBound; |
---|
| 1990 | } |
---|
| 1991 | } else if (upperValue == COIN_DBL_MAX) { |
---|
| 1992 | abcStatus = atLowerBound; |
---|
| 1993 | } else if (lowerValue == upperValue) { |
---|
| 1994 | abcStatus = isFixed; |
---|
| 1995 | } else if (abcStatus == isFree) { |
---|
| 1996 | abcStatus = superBasic; |
---|
| 1997 | } |
---|
| 1998 | switch (abcStatus) { |
---|
| 1999 | case isFixed: |
---|
| 2000 | case atLowerBound: |
---|
| 2001 | solutionSaved[iColumn] = lowerValue; |
---|
| 2002 | break; |
---|
| 2003 | case atUpperBound: |
---|
| 2004 | solutionSaved[iColumn] = upperValue; |
---|
| 2005 | break; |
---|
| 2006 | default: |
---|
| 2007 | break; |
---|
| 2008 | } |
---|
| 2009 | } else { |
---|
| 2010 | // basic |
---|
| 2011 | if (numberBasic < numberRows_) |
---|
| 2012 | abcPivotVariable_[numberBasic++] = iColumn + maximumAbcNumberRows_; |
---|
| 2013 | else |
---|
| 2014 | abcStatus = superBasic; |
---|
| 2015 | } |
---|
| 2016 | putStatus[iColumn] = abcStatus; |
---|
| 2017 | } |
---|
| 2018 | } |
---|
| 2019 | numberStructurals += numberBasic - numberSaved; |
---|
| 2020 | delete simplex[iBlock]; |
---|
[1878] | 2021 | } |
---|
[2385] | 2022 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 2023 | printf("%d structurals put into basis\n", numberStructurals); |
---|
[1878] | 2024 | #endif |
---|
[2385] | 2025 | if (numberBasic < numberRows_) { |
---|
| 2026 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2027 | AbcSimplex::Status status = getInternalStatus(iRow); |
---|
| 2028 | if (status != AbcSimplex::basic) { |
---|
| 2029 | abcPivotVariable_[numberBasic++] = iRow; |
---|
| 2030 | setInternalStatus(iRow, basic); |
---|
| 2031 | if (numberBasic == numberRows_) |
---|
| 2032 | break; |
---|
| 2033 | } |
---|
| 2034 | } |
---|
[1878] | 2035 | } |
---|
[2385] | 2036 | assert(numberBasic == numberRows_); |
---|
[1878] | 2037 | #if 0 |
---|
| 2038 | int n2=0; |
---|
| 2039 | for (int i=0;i<numberTotal_;i++) { |
---|
| 2040 | if (getInternalStatus(i)==basic) |
---|
| 2041 | n2++; |
---|
| 2042 | } |
---|
| 2043 | assert (n2==numberRows_); |
---|
| 2044 | std::sort(abcPivotVariable_,abcPivotVariable_+numberRows_); |
---|
| 2045 | n2=-1; |
---|
| 2046 | for (int i=0;i<numberRows_;i++) { |
---|
| 2047 | assert (abcPivotVariable_[i]>n2); |
---|
| 2048 | n2=abcPivotVariable_[i]; |
---|
| 2049 | } |
---|
| 2050 | #endif |
---|
[2385] | 2051 | delete[] simplex; |
---|
[1878] | 2052 | return; |
---|
| 2053 | } else { |
---|
| 2054 | // try another |
---|
[2385] | 2055 | type = 2; |
---|
[1878] | 2056 | } |
---|
| 2057 | } |
---|
[2385] | 2058 | if (type == 1) { |
---|
| 2059 | const double *linearObjective = abcCost_ + maximumAbcNumberRows_; |
---|
| 2060 | double gamma = 0.0; |
---|
| 2061 | int numberTotal = numberRows_ + numberColumns_; |
---|
| 2062 | double *q = new double[numberTotal]; |
---|
| 2063 | double *v = q + numberColumns_; |
---|
| 2064 | int *which = new int[numberTotal + 3 * numberRows_]; |
---|
| 2065 | int *ii = which + numberColumns_; |
---|
| 2066 | int *r = ii + numberRows_; |
---|
| 2067 | int *pivoted = r + numberRows_; |
---|
[1878] | 2068 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
[2385] | 2069 | gamma = CoinMax(gamma, linearObjective[iColumn]); |
---|
[1878] | 2070 | } |
---|
| 2071 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2072 | double lowerBound = abcLower_[iRow]; |
---|
| 2073 | double upperBound = abcUpper_[iRow]; |
---|
[2385] | 2074 | pivoted[iRow] = -1; |
---|
| 2075 | ii[iRow] = 0; |
---|
| 2076 | r[iRow] = 0; |
---|
| 2077 | v[iRow] = COIN_DBL_MAX; |
---|
| 2078 | if (lowerBound == upperBound) |
---|
| 2079 | continue; |
---|
| 2080 | if (lowerBound <= 0.0 && upperBound >= 0.0) { |
---|
| 2081 | pivoted[iRow] = iRow; |
---|
| 2082 | ii[iRow] = 1; |
---|
| 2083 | r[iRow] = 1; |
---|
[1878] | 2084 | } |
---|
| 2085 | } |
---|
[2385] | 2086 | int nPossible = 0; |
---|
| 2087 | int lastPossible = 0; |
---|
[1878] | 2088 | double cMaxDiv; |
---|
| 2089 | if (gamma) |
---|
[2385] | 2090 | cMaxDiv = 1.0 / (1000.0 * gamma); |
---|
[1878] | 2091 | else |
---|
[2385] | 2092 | cMaxDiv = 1.0; |
---|
| 2093 | const double *columnLower = abcLower_ + maximumAbcNumberRows_; |
---|
| 2094 | const double *columnUpper = abcUpper_ + maximumAbcNumberRows_; |
---|
| 2095 | for (int iPass = 0; iPass < 3; iPass++) { |
---|
[1878] | 2096 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
[2385] | 2097 | double lowerBound = columnLower[iColumn]; |
---|
| 2098 | double upperBound = columnUpper[iColumn]; |
---|
| 2099 | if (lowerBound == upperBound) |
---|
| 2100 | continue; |
---|
| 2101 | double qValue; |
---|
| 2102 | if (lowerBound > -1.0e20) { |
---|
| 2103 | if (upperBound < 1.0e20) { |
---|
| 2104 | // both |
---|
| 2105 | qValue = lowerBound - upperBound; |
---|
| 2106 | if (iPass != 2) |
---|
| 2107 | qValue = COIN_DBL_MAX; |
---|
| 2108 | } else { |
---|
| 2109 | // just lower |
---|
| 2110 | qValue = lowerBound; |
---|
| 2111 | if (iPass != 1) |
---|
| 2112 | qValue = COIN_DBL_MAX; |
---|
| 2113 | } |
---|
| 2114 | } else { |
---|
| 2115 | if (upperBound < 1.0e20) { |
---|
| 2116 | // just upper |
---|
| 2117 | qValue = -upperBound; |
---|
| 2118 | if (iPass != 1) |
---|
| 2119 | qValue = COIN_DBL_MAX; |
---|
| 2120 | } else { |
---|
| 2121 | // free |
---|
| 2122 | qValue = 0.0; |
---|
| 2123 | if (iPass != 0) |
---|
| 2124 | qValue = COIN_DBL_MAX; |
---|
| 2125 | } |
---|
| 2126 | } |
---|
| 2127 | if (qValue != COIN_DBL_MAX) { |
---|
| 2128 | which[nPossible] = iColumn; |
---|
| 2129 | q[nPossible++] = qValue + linearObjective[iColumn] * cMaxDiv; |
---|
| 2130 | ; |
---|
| 2131 | } |
---|
[1878] | 2132 | } |
---|
[2385] | 2133 | CoinSort_2(q + lastPossible, q + nPossible, which + lastPossible); |
---|
| 2134 | lastPossible = nPossible; |
---|
[1878] | 2135 | } |
---|
[2385] | 2136 | const double *element = abcMatrix_->getPackedMatrix()->getElements(); |
---|
| 2137 | const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts(); |
---|
[1878] | 2138 | //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths(); |
---|
[2385] | 2139 | const int *row = abcMatrix_->getPackedMatrix()->getIndices(); |
---|
| 2140 | int nPut = 0; |
---|
| 2141 | for (int i = 0; i < nPossible; i++) { |
---|
| 2142 | int iColumn = which[i]; |
---|
| 2143 | double maxAlpha = 0.0; |
---|
| 2144 | int kRow = -1; |
---|
| 2145 | double alternativeAlpha = 0.0; |
---|
| 2146 | int jRow = -1; |
---|
| 2147 | bool canTake = true; |
---|
| 2148 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { |
---|
| 2149 | int iRow = row[j]; |
---|
| 2150 | double alpha = fabs(element[j]); |
---|
| 2151 | if (alpha > 0.01 * v[iRow]) { |
---|
| 2152 | canTake = false; |
---|
| 2153 | } else if (!ii[iRow] && alpha > alternativeAlpha) { |
---|
| 2154 | alternativeAlpha = alpha; |
---|
| 2155 | jRow = iRow; |
---|
| 2156 | } |
---|
| 2157 | if (!r[iRow] && alpha > maxAlpha) { |
---|
| 2158 | maxAlpha = alpha; |
---|
| 2159 | kRow = iRow; |
---|
| 2160 | } |
---|
[1878] | 2161 | } |
---|
| 2162 | // only really works if scaled |
---|
[2385] | 2163 | if (maxAlpha > 0.99) { |
---|
| 2164 | pivoted[kRow] = iColumn + maximumAbcNumberRows_; |
---|
| 2165 | v[kRow] = maxAlpha; |
---|
| 2166 | ii[kRow] = 1; |
---|
| 2167 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { |
---|
| 2168 | int iRow = row[j]; |
---|
| 2169 | r[iRow]++; |
---|
| 2170 | } |
---|
| 2171 | nPut++; |
---|
| 2172 | } else if (canTake && jRow >= 0) { |
---|
| 2173 | pivoted[jRow] = iColumn + maximumAbcNumberRows_; |
---|
| 2174 | v[jRow] = maxAlpha; |
---|
| 2175 | ii[jRow] = 1; |
---|
| 2176 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { |
---|
| 2177 | int iRow = row[j]; |
---|
| 2178 | r[iRow]++; |
---|
| 2179 | } |
---|
| 2180 | nPut++; |
---|
[1878] | 2181 | } |
---|
| 2182 | } |
---|
[2385] | 2183 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2184 | int iSequence = pivoted[iRow]; |
---|
| 2185 | if (iSequence >= 0 && iSequence < numberColumns_) { |
---|
| 2186 | abcPivotVariable_[iRow] = iSequence; |
---|
| 2187 | if (fabs(abcLower_[iRow]) < fabs(abcUpper_[iRow])) { |
---|
| 2188 | setInternalStatus(iRow, atLowerBound); |
---|
| 2189 | abcSolution_[iRow] = abcLower_[iRow]; |
---|
| 2190 | solutionSaved_[iRow] = abcLower_[iRow]; |
---|
| 2191 | } else { |
---|
| 2192 | setInternalStatus(iRow, atUpperBound); |
---|
| 2193 | abcSolution_[iRow] = abcUpper_[iRow]; |
---|
| 2194 | solutionSaved_[iRow] = abcUpper_[iRow]; |
---|
| 2195 | } |
---|
| 2196 | setInternalStatus(iSequence, basic); |
---|
[1878] | 2197 | } |
---|
| 2198 | } |
---|
[2385] | 2199 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 2200 | printf("%d put into basis\n", nPut); |
---|
[1878] | 2201 | #endif |
---|
[2385] | 2202 | delete[] q; |
---|
| 2203 | delete[] which; |
---|
| 2204 | } else if (type == 2) { |
---|
[1878] | 2205 | //return; |
---|
[2385] | 2206 | int numberBad = 0; |
---|
| 2207 | CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_); |
---|
[1878] | 2208 | // Work on savedSolution and current |
---|
[2385] | 2209 | int iVector = getAvailableArray(); |
---|
[1878] | 2210 | #define ALLOW_BAD_DJS |
---|
| 2211 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2212 | double *modDj = usefulArray_[iVector].denseVector(); |
---|
[1878] | 2213 | #endif |
---|
[2385] | 2214 | for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) { |
---|
| 2215 | double dj = abcDj_[iSequence]; |
---|
| 2216 | modDj[iSequence] = 0.0; |
---|
| 2217 | if (getInternalStatus(iSequence) == atLowerBound) { |
---|
| 2218 | if (dj < -dualTolerance_) { |
---|
| 2219 | double costThru = -dj * (abcUpper_[iSequence] - abcLower_[iSequence]); |
---|
| 2220 | if (costThru < dualBound_) { |
---|
| 2221 | // flip |
---|
| 2222 | setInternalStatus(iSequence, atUpperBound); |
---|
| 2223 | solutionSaved_[iSequence] = abcUpper_[iSequence]; |
---|
| 2224 | abcSolution_[iSequence] = abcUpper_[iSequence]; |
---|
| 2225 | } else { |
---|
| 2226 | numberBad++; |
---|
[1878] | 2227 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2228 | modDj[iSequence] = dj; |
---|
| 2229 | dj = 0.0; |
---|
[1878] | 2230 | #else |
---|
[2385] | 2231 | break; |
---|
[1878] | 2232 | #endif |
---|
[2385] | 2233 | } |
---|
| 2234 | } else { |
---|
| 2235 | dj = CoinMax(dj, 0.0); |
---|
| 2236 | } |
---|
| 2237 | } else if (getInternalStatus(iSequence) == atLowerBound) { |
---|
| 2238 | if (dj > dualTolerance_) { |
---|
| 2239 | double costThru = dj * (abcUpper_[iSequence] - abcLower_[iSequence]); |
---|
| 2240 | if (costThru < dualBound_) { |
---|
| 2241 | assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10); |
---|
| 2242 | // flip |
---|
| 2243 | setInternalStatus(iSequence, atLowerBound); |
---|
| 2244 | solutionSaved_[iSequence] = abcLower_[iSequence]; |
---|
| 2245 | abcSolution_[iSequence] = abcLower_[iSequence]; |
---|
| 2246 | } else { |
---|
| 2247 | numberBad++; |
---|
[1878] | 2248 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2249 | modDj[iSequence] = dj; |
---|
| 2250 | dj = 0.0; |
---|
[1878] | 2251 | #else |
---|
[2385] | 2252 | break; |
---|
[1878] | 2253 | #endif |
---|
[2385] | 2254 | } |
---|
| 2255 | } else { |
---|
| 2256 | dj = CoinMin(dj, 0.0); |
---|
| 2257 | } |
---|
[1878] | 2258 | } else { |
---|
[2385] | 2259 | if (fabs(dj) < dualTolerance_) { |
---|
| 2260 | dj = 0.0; |
---|
| 2261 | } else { |
---|
| 2262 | numberBad++; |
---|
[1878] | 2263 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2264 | modDj[iSequence] = dj; |
---|
| 2265 | dj = 0.0; |
---|
[1878] | 2266 | #else |
---|
[2385] | 2267 | break; |
---|
[1878] | 2268 | #endif |
---|
[2385] | 2269 | } |
---|
[1878] | 2270 | } |
---|
[2385] | 2271 | abcDj_[iSequence] = dj; |
---|
[1878] | 2272 | } |
---|
| 2273 | #ifndef ALLOW_BAD_DJS |
---|
| 2274 | if (numberBad) { |
---|
| 2275 | //CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_); |
---|
[2385] | 2276 | return; |
---|
[1878] | 2277 | } |
---|
| 2278 | #endif |
---|
| 2279 | if (!abcMatrix_->gotRowCopy()) |
---|
| 2280 | abcMatrix_->createRowCopy(); |
---|
| 2281 | //const double * element = abcMatrix_->getPackedMatrix()->getElements(); |
---|
[2385] | 2282 | const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_; |
---|
| 2283 | const int *row = abcMatrix_->getPackedMatrix()->getIndices(); |
---|
| 2284 | const double *elementByRow = abcMatrix_->rowElements(); |
---|
| 2285 | const CoinBigIndex *rowStart = abcMatrix_->rowStart(); |
---|
| 2286 | const CoinBigIndex *rowEnd = abcMatrix_->rowEnd(); |
---|
| 2287 | const int *column = abcMatrix_->rowColumns(); |
---|
| 2288 | CoinAbcMemset0(solutionBasic_, numberRows_); |
---|
| 2289 | CoinAbcMemcpy(lowerBasic_, abcLower_, numberRows_); |
---|
| 2290 | CoinAbcMemcpy(upperBasic_, abcUpper_, numberRows_); |
---|
| 2291 | abcMatrix_->timesIncludingSlacks(-1.0, solutionSaved_, solutionBasic_); |
---|
| 2292 | const double multiplier[] = { 1.0, -1.0 }; |
---|
| 2293 | int nBasic = 0; |
---|
| 2294 | int *index = usefulArray_[iVector].getIndices(); |
---|
| 2295 | int iVector2 = getAvailableArray(); |
---|
| 2296 | int *index2 = usefulArray_[iVector2].getIndices(); |
---|
| 2297 | double *sort = usefulArray_[iVector2].denseVector(); |
---|
| 2298 | int average = CoinMax(5, rowEnd[numberRows_ - 1] / (8 * numberRows_)); |
---|
| 2299 | int nPossible = 0; |
---|
| 2300 | if (numberRows_ > 10000) { |
---|
[1878] | 2301 | // allow more |
---|
[2385] | 2302 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2303 | double solutionValue = solutionBasic_[iRow]; |
---|
| 2304 | double infeasibility = 0.0; |
---|
| 2305 | double lowerValue = lowerBasic_[iRow]; |
---|
| 2306 | double upperValue = upperBasic_[iRow]; |
---|
| 2307 | if (solutionValue < lowerValue - primalTolerance_) { |
---|
| 2308 | infeasibility = -(lowerValue - solutionValue); |
---|
| 2309 | } else if (solutionValue > upperValue + primalTolerance_) { |
---|
| 2310 | infeasibility = upperValue - solutionValue; |
---|
| 2311 | } |
---|
| 2312 | int length = rowEnd[iRow] - rowStart[iRow]; |
---|
| 2313 | if (infeasibility) |
---|
| 2314 | index2[nPossible++] = length; |
---|
[1878] | 2315 | } |
---|
[2385] | 2316 | std::sort(index2, index2 + nPossible); |
---|
[1878] | 2317 | // see how much we need to get numberRows/10 or nPossible/3 |
---|
[2385] | 2318 | average = CoinMax(average, index2[CoinMin(numberRows_ / 10, nPossible / 3)]); |
---|
| 2319 | nPossible = 0; |
---|
[1878] | 2320 | } |
---|
[2385] | 2321 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2322 | double solutionValue = solutionBasic_[iRow]; |
---|
| 2323 | double infeasibility = 0.0; |
---|
| 2324 | double lowerValue = lowerBasic_[iRow]; |
---|
| 2325 | double upperValue = upperBasic_[iRow]; |
---|
| 2326 | if (solutionValue < lowerValue - primalTolerance_) { |
---|
| 2327 | infeasibility = -(lowerValue - solutionValue); |
---|
| 2328 | } else if (solutionValue > upperValue + primalTolerance_) { |
---|
| 2329 | infeasibility = upperValue - solutionValue; |
---|
[1878] | 2330 | } |
---|
[2385] | 2331 | int length = rowEnd[iRow] - rowStart[iRow]; |
---|
| 2332 | if (infeasibility && length < average) { |
---|
| 2333 | index2[nPossible] = iRow; |
---|
| 2334 | sort[nPossible++] = 1.0e5 * infeasibility - iRow; |
---|
| 2335 | //sort[nPossible++]=1.0e9*length-iRow;//infeasibility; |
---|
[1878] | 2336 | } |
---|
| 2337 | } |
---|
[2385] | 2338 | CoinSort_2(sort, sort + nPossible, index2); |
---|
| 2339 | for (int iWhich = 0; iWhich < nPossible; iWhich++) { |
---|
[1878] | 2340 | int iRow = index2[iWhich]; |
---|
[2385] | 2341 | sort[iWhich] = 0.0; |
---|
[1878] | 2342 | if (abcDj_[iRow]) |
---|
[2385] | 2343 | continue; // marked as invalid |
---|
| 2344 | double solutionValue = solutionBasic_[iRow]; |
---|
| 2345 | double infeasibility = 0.0; |
---|
| 2346 | double lowerValue = lowerBasic_[iRow]; |
---|
| 2347 | double upperValue = upperBasic_[iRow]; |
---|
| 2348 | if (solutionValue < lowerValue - primalTolerance_) { |
---|
| 2349 | infeasibility = lowerValue - solutionValue; |
---|
| 2350 | } else if (solutionValue > upperValue + primalTolerance_) { |
---|
| 2351 | infeasibility = upperValue - solutionValue; |
---|
[1878] | 2352 | } |
---|
[2385] | 2353 | assert(infeasibility); |
---|
| 2354 | double direction = infeasibility > 0 ? 1.0 : -1.0; |
---|
[1878] | 2355 | infeasibility *= direction; |
---|
[2385] | 2356 | int whichColumn = -1; |
---|
| 2357 | double upperTheta = 1.0e30; |
---|
| 2358 | int whichColumn2 = -1; |
---|
| 2359 | double upperTheta2 = 1.0e30; |
---|
| 2360 | double costThru = 0.0; |
---|
| 2361 | int nThru = 0; |
---|
| 2362 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 2363 | int iSequence = column[j]; |
---|
| 2364 | assert(iSequence >= maximumAbcNumberRows_); |
---|
| 2365 | double dj = abcDj_[iSequence]; |
---|
| 2366 | double tableauValue = -elementByRow[j] * direction; |
---|
| 2367 | unsigned char iStatus = internalStatus_[iSequence] & 7; |
---|
| 2368 | if ((iStatus & 4) == 0) { |
---|
| 2369 | double mult = multiplier[iStatus]; |
---|
| 2370 | double alpha = tableauValue * mult; |
---|
| 2371 | double oldValue = dj * mult; |
---|
| 2372 | assert(oldValue > -1.0e-2); |
---|
| 2373 | double value = oldValue - upperTheta * alpha; |
---|
| 2374 | if (value < 0.0) { |
---|
| 2375 | upperTheta2 = upperTheta; |
---|
| 2376 | whichColumn2 = whichColumn; |
---|
| 2377 | costThru = alpha * (abcUpper_[iSequence] - abcLower_[iSequence]); |
---|
| 2378 | nThru = 0; |
---|
| 2379 | upperTheta = oldValue / alpha; |
---|
| 2380 | whichColumn = iSequence; |
---|
| 2381 | } else if (oldValue - upperTheta2 * alpha < 0.0) { |
---|
| 2382 | costThru += alpha * (abcUpper_[iSequence] - abcLower_[iSequence]); |
---|
| 2383 | index[nThru++] = iSequence; |
---|
| 2384 | } |
---|
| 2385 | } else if (iStatus < 6) { |
---|
| 2386 | upperTheta = -1.0; |
---|
| 2387 | upperTheta2 = elementByRow[j]; |
---|
| 2388 | whichColumn = iSequence; |
---|
| 2389 | break; |
---|
| 2390 | } |
---|
[1878] | 2391 | } |
---|
[2385] | 2392 | if (whichColumn < 0) |
---|
| 2393 | continue; |
---|
| 2394 | if (upperTheta != -1.0) { |
---|
| 2395 | assert(upperTheta >= 0.0); |
---|
| 2396 | if (costThru < infeasibility && whichColumn2 >= 0) { |
---|
| 2397 | index[nThru++] = whichColumn; |
---|
| 2398 | for (int i = 0; i < nThru; i++) { |
---|
| 2399 | int iSequence = index[i]; |
---|
| 2400 | assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10); |
---|
| 2401 | // flip and use previous |
---|
| 2402 | if (getInternalStatus(iSequence) == atLowerBound) { |
---|
| 2403 | setInternalStatus(iSequence, atUpperBound); |
---|
| 2404 | solutionSaved_[iSequence] = abcUpper_[iSequence]; |
---|
| 2405 | abcSolution_[iSequence] = abcUpper_[iSequence]; |
---|
| 2406 | } else { |
---|
| 2407 | setInternalStatus(iSequence, atLowerBound); |
---|
| 2408 | solutionSaved_[iSequence] = abcLower_[iSequence]; |
---|
| 2409 | abcSolution_[iSequence] = abcLower_[iSequence]; |
---|
| 2410 | } |
---|
| 2411 | } |
---|
| 2412 | whichColumn = whichColumn2; |
---|
| 2413 | upperTheta = upperTheta2; |
---|
| 2414 | } |
---|
[1878] | 2415 | } else { |
---|
[2385] | 2416 | // free coming in |
---|
| 2417 | upperTheta = (abcDj_[whichColumn] * direction) / upperTheta2; |
---|
[1878] | 2418 | } |
---|
| 2419 | // update djs |
---|
| 2420 | upperTheta *= -direction; |
---|
[2385] | 2421 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 2422 | int iSequence = column[j]; |
---|
| 2423 | double dj = abcDj_[iSequence]; |
---|
| 2424 | double tableauValue = elementByRow[j]; |
---|
| 2425 | dj -= upperTheta * tableauValue; |
---|
| 2426 | unsigned char iStatus = internalStatus_[iSequence] & 7; |
---|
| 2427 | if ((iStatus & 4) == 0) { |
---|
| 2428 | if (!iStatus) { |
---|
| 2429 | assert(dj > -1.0e-2); |
---|
| 2430 | dj = CoinMax(dj, 0.0); |
---|
[1878] | 2431 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2432 | if (numberBad && modDj[iSequence]) { |
---|
| 2433 | double bad = modDj[iSequence]; |
---|
| 2434 | if (dj + bad >= 0.0) { |
---|
| 2435 | numberBad--; |
---|
| 2436 | modDj[iSequence] = 0.0; |
---|
| 2437 | dj += bad; |
---|
| 2438 | } else { |
---|
| 2439 | modDj[iSequence] += dj; |
---|
| 2440 | dj = 0.0; |
---|
| 2441 | } |
---|
| 2442 | } |
---|
[1878] | 2443 | #endif |
---|
[2385] | 2444 | } else { |
---|
| 2445 | assert(dj < 1.0e-2); |
---|
| 2446 | dj = CoinMin(dj, 0.0); |
---|
[1878] | 2447 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2448 | if (numberBad && modDj[iSequence]) { |
---|
| 2449 | double bad = modDj[iSequence]; |
---|
| 2450 | if (dj + bad <= 0.0) { |
---|
| 2451 | numberBad--; |
---|
| 2452 | modDj[iSequence] = 0.0; |
---|
| 2453 | dj += bad; |
---|
| 2454 | } else { |
---|
| 2455 | modDj[iSequence] += dj; |
---|
| 2456 | dj = 0.0; |
---|
| 2457 | } |
---|
| 2458 | } |
---|
[1878] | 2459 | #endif |
---|
[2385] | 2460 | } |
---|
| 2461 | } else if (iStatus < 6) { |
---|
| 2462 | assert(fabs(dj) < 1.0e-4); |
---|
| 2463 | dj = 0.0; |
---|
| 2464 | } |
---|
| 2465 | abcDj_[iSequence] = dj; |
---|
[1878] | 2466 | } |
---|
| 2467 | // do basis |
---|
[2385] | 2468 | if (direction > 0.0) { |
---|
| 2469 | if (upperBasic_[iRow] > lowerBasic_[iRow]) |
---|
| 2470 | setInternalStatus(iRow, atLowerBound); |
---|
| 2471 | else |
---|
| 2472 | setInternalStatus(iRow, isFixed); |
---|
| 2473 | solutionSaved_[iRow] = abcLower_[iRow]; |
---|
| 2474 | abcSolution_[iRow] = abcLower_[iRow]; |
---|
[1878] | 2475 | } else { |
---|
[2385] | 2476 | if (upperBasic_[iRow] > lowerBasic_[iRow]) |
---|
| 2477 | setInternalStatus(iRow, atUpperBound); |
---|
| 2478 | else |
---|
| 2479 | setInternalStatus(iRow, isFixed); |
---|
| 2480 | solutionSaved_[iRow] = abcUpper_[iRow]; |
---|
| 2481 | abcSolution_[iRow] = abcUpper_[iRow]; |
---|
[1878] | 2482 | } |
---|
[2385] | 2483 | setInternalStatus(whichColumn, basic); |
---|
| 2484 | abcPivotVariable_[iRow] = whichColumn; |
---|
[1878] | 2485 | nBasic++; |
---|
| 2486 | // mark rows |
---|
[2385] | 2487 | for (CoinBigIndex j = columnStart[whichColumn]; j < columnStart[whichColumn + 1]; j++) { |
---|
| 2488 | int jRow = row[j]; |
---|
| 2489 | abcDj_[jRow] = 1.0; |
---|
[1878] | 2490 | } |
---|
| 2491 | } |
---|
| 2492 | #ifdef ALLOW_BAD_DJS |
---|
[2385] | 2493 | CoinAbcMemset0(modDj + maximumAbcNumberRows_, numberColumns_); |
---|
[1878] | 2494 | #endif |
---|
| 2495 | setAvailableArray(iVector); |
---|
| 2496 | setAvailableArray(iVector2); |
---|
[2385] | 2497 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 2498 | printf("dual crash put %d in basis\n", nBasic); |
---|
[1878] | 2499 | #endif |
---|
| 2500 | } else { |
---|
[2385] | 2501 | assert((stateOfProblem_ & VALUES_PASS2) != 0); |
---|
[1878] | 2502 | // The idea is to put as many likely variables into basis as possible |
---|
[2385] | 2503 | int n = 0; |
---|
| 2504 | int iVector = getAvailableArray(); |
---|
| 2505 | int *index = usefulArray_[iVector].getIndices(); |
---|
| 2506 | double *array = usefulArray_[iVector].denseVector(); |
---|
| 2507 | int iVector2 = getAvailableArray(); |
---|
| 2508 | int *index2 = usefulArray_[iVector].getIndices(); |
---|
| 2509 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[1878] | 2510 | double dj = djSaved_[iSequence]; |
---|
| 2511 | double value = solution_[iSequence]; |
---|
| 2512 | double lower = abcLower_[iSequence]; |
---|
| 2513 | double upper = abcUpper_[iSequence]; |
---|
[2385] | 2514 | double gapUp = CoinMin(1.0e3, upper - value); |
---|
| 2515 | assert(gapUp >= -1.0e-3); |
---|
| 2516 | gapUp = CoinMax(gapUp, 0.0); |
---|
| 2517 | double gapDown = CoinMin(1.0e3, value - lower); |
---|
| 2518 | assert(gapDown >= -1.0e-3); |
---|
| 2519 | gapDown = CoinMax(gapDown, 0.0); |
---|
| 2520 | double measure = (CoinMin(gapUp, gapDown) + 1.0e-6) / (fabs(dj) + 1.0e-6); |
---|
| 2521 | if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) { |
---|
| 2522 | // set to ub |
---|
| 2523 | setInternalStatus(iSequence, atUpperBound); |
---|
| 2524 | solutionSaved_[iSequence] = abcUpper_[iSequence]; |
---|
| 2525 | abcSolution_[iSequence] = abcUpper_[iSequence]; |
---|
| 2526 | } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) { |
---|
| 2527 | // set to lb |
---|
| 2528 | setInternalStatus(iSequence, atLowerBound); |
---|
| 2529 | solutionSaved_[iSequence] = abcLower_[iSequence]; |
---|
| 2530 | abcSolution_[iSequence] = abcLower_[iSequence]; |
---|
| 2531 | } else if (upper > lower) { |
---|
| 2532 | // set to nearest |
---|
| 2533 | if (gapUp < gapDown) { |
---|
| 2534 | // set to ub |
---|
| 2535 | setInternalStatus(iSequence, atUpperBound); |
---|
| 2536 | solutionSaved_[iSequence] = abcUpper_[iSequence]; |
---|
| 2537 | abcSolution_[iSequence] = abcUpper_[iSequence]; |
---|
| 2538 | } else { |
---|
| 2539 | // set to lb |
---|
| 2540 | setInternalStatus(iSequence, atLowerBound); |
---|
| 2541 | solutionSaved_[iSequence] = abcLower_[iSequence]; |
---|
| 2542 | abcSolution_[iSequence] = abcLower_[iSequence]; |
---|
| 2543 | } |
---|
| 2544 | array[n] = -measure; |
---|
| 2545 | index[n++] = iSequence; |
---|
[1878] | 2546 | } |
---|
| 2547 | } |
---|
| 2548 | // set slacks basic |
---|
[2385] | 2549 | memset(internalStatus_, 6, numberRows_); |
---|
| 2550 | CoinSort_2(solution_, solution_ + n, index); |
---|
| 2551 | CoinAbcMemset0(array, n); |
---|
| 2552 | for (int i = 0; i < numberRows_; i++) |
---|
| 2553 | index2[i] = 0; |
---|
[1878] | 2554 | // first put all possible slacks in |
---|
[2385] | 2555 | int n2 = 0; |
---|
| 2556 | for (int i = 0; i < n; i++) { |
---|
| 2557 | int iSequence = index[i]; |
---|
| 2558 | if (iSequence < numberRows_) { |
---|
| 2559 | index2[iSequence] = numberRows_; |
---|
[1878] | 2560 | } else { |
---|
[2385] | 2561 | index[n2++] = iSequence; |
---|
[1878] | 2562 | } |
---|
| 2563 | } |
---|
[2385] | 2564 | n = n2; |
---|
| 2565 | int numberIn = 0; |
---|
| 2566 | const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_; |
---|
[1878] | 2567 | //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths(); |
---|
[2385] | 2568 | const int *row = abcMatrix_->getPackedMatrix()->getIndices(); |
---|
[1878] | 2569 | if (!abcMatrix_->gotRowCopy()) |
---|
| 2570 | abcMatrix_->createRowCopy(); |
---|
| 2571 | //const CoinBigIndex * rowStart = abcMatrix_->rowStart(); |
---|
| 2572 | //const CoinBigIndex * rowEnd = abcMatrix_->rowEnd(); |
---|
[2385] | 2573 | for (int i = 0; i < n; i++) { |
---|
| 2574 | int iSequence = index[i]; |
---|
| 2575 | int bestRow = -1; |
---|
| 2576 | int bestCount = numberRows_ + 1; |
---|
| 2577 | for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) { |
---|
| 2578 | int iRow = row[j]; |
---|
| 2579 | if (!abcMatrix_->gotRowCopy()) |
---|
| 2580 | abcMatrix_->createRowCopy(); |
---|
| 2581 | const CoinBigIndex *rowStart = abcMatrix_->rowStart(); |
---|
| 2582 | const CoinBigIndex *rowEnd = abcMatrix_->rowEnd(); |
---|
| 2583 | if (!index2[iRow]) { |
---|
| 2584 | int length = rowEnd[iRow] - rowStart[iRow]; |
---|
| 2585 | if (length < bestCount) { |
---|
| 2586 | bestCount = length; |
---|
| 2587 | bestRow = iRow; |
---|
| 2588 | } |
---|
| 2589 | } |
---|
[1878] | 2590 | } |
---|
[2385] | 2591 | if (bestRow >= 0) { |
---|
| 2592 | numberIn++; |
---|
| 2593 | for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) { |
---|
| 2594 | int iRow = row[j]; |
---|
| 2595 | index2[iRow]++; |
---|
| 2596 | } |
---|
| 2597 | setInternalStatus(iSequence, basic); |
---|
| 2598 | abcPivotVariable_[bestRow] = iSequence; |
---|
| 2599 | double dj = djSaved_[bestRow]; |
---|
| 2600 | double value = solution_[bestRow]; |
---|
| 2601 | double lower = abcLower_[bestRow]; |
---|
| 2602 | double upper = abcUpper_[bestRow]; |
---|
| 2603 | double gapUp = CoinMax(CoinMin(1.0e3, upper - value), 0.0); |
---|
| 2604 | double gapDown = CoinMax(CoinMin(1.0e3, value - lower), 0.0); |
---|
| 2605 | //double measure = (CoinMin(gapUp,gapDown)+1.0e-6)/(fabs(dj)+1.0e-6); |
---|
| 2606 | if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) { |
---|
| 2607 | // set to ub |
---|
| 2608 | setInternalStatus(bestRow, atUpperBound); |
---|
| 2609 | solutionSaved_[bestRow] = abcUpper_[bestRow]; |
---|
| 2610 | abcSolution_[bestRow] = abcUpper_[bestRow]; |
---|
| 2611 | } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) { |
---|
| 2612 | // set to lb |
---|
| 2613 | setInternalStatus(bestRow, atLowerBound); |
---|
| 2614 | solutionSaved_[bestRow] = abcLower_[bestRow]; |
---|
| 2615 | abcSolution_[bestRow] = abcLower_[bestRow]; |
---|
| 2616 | } else if (upper > lower) { |
---|
| 2617 | // set to nearest |
---|
| 2618 | if (gapUp < gapDown) { |
---|
| 2619 | // set to ub |
---|
| 2620 | setInternalStatus(bestRow, atUpperBound); |
---|
| 2621 | solutionSaved_[bestRow] = abcUpper_[bestRow]; |
---|
| 2622 | abcSolution_[bestRow] = abcUpper_[bestRow]; |
---|
| 2623 | } else { |
---|
| 2624 | // set to lb |
---|
| 2625 | setInternalStatus(bestRow, atLowerBound); |
---|
| 2626 | solutionSaved_[bestRow] = abcLower_[bestRow]; |
---|
| 2627 | abcSolution_[bestRow] = abcLower_[bestRow]; |
---|
| 2628 | } |
---|
| 2629 | } |
---|
[1878] | 2630 | } |
---|
| 2631 | } |
---|
[2385] | 2632 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 2633 | printf("Idiot crash put %d in basis\n", numberIn); |
---|
[1878] | 2634 | #endif |
---|
| 2635 | setAvailableArray(iVector); |
---|
| 2636 | setAvailableArray(iVector2); |
---|
[2385] | 2637 | delete[] solution_; |
---|
| 2638 | solution_ = NULL; |
---|
[1878] | 2639 | } |
---|
| 2640 | } |
---|
| 2641 | /* Puts more stuff in basis |
---|
| 2642 | 1 bit set - do even if basis exists |
---|
| 2643 | 2 bit set - don't bother staying triangular |
---|
| 2644 | */ |
---|
[2385] | 2645 | void AbcSimplex::putStuffInBasis(int type) |
---|
[1878] | 2646 | { |
---|
| 2647 | int i; |
---|
[2385] | 2648 | for (i = 0; i < numberRows_; i++) { |
---|
| 2649 | if (getInternalStatus(i) != basic) |
---|
[1878] | 2650 | break; |
---|
| 2651 | } |
---|
[2385] | 2652 | if (i < numberRows_ && (type & 1) == 0) |
---|
[1878] | 2653 | return; |
---|
[2385] | 2654 | int iVector = getAvailableArray(); |
---|
[1878] | 2655 | // Column part is mustColumnIn |
---|
[2385] | 2656 | int *COIN_RESTRICT mustRowOut = usefulArray_[iVector].getIndices(); |
---|
[1878] | 2657 | if (!abcMatrix_->gotRowCopy()) |
---|
| 2658 | abcMatrix_->createRowCopy(); |
---|
[2385] | 2659 | const double *elementByRow = abcMatrix_->rowElements(); |
---|
| 2660 | const CoinBigIndex *rowStart = abcMatrix_->rowStart(); |
---|
| 2661 | const CoinBigIndex *rowEnd = abcMatrix_->rowEnd(); |
---|
| 2662 | const int *column = abcMatrix_->rowColumns(); |
---|
| 2663 | for (int i = 0; i < numberTotal_; i++) |
---|
| 2664 | mustRowOut[i] = -1; |
---|
| 2665 | int nPossible = 0; |
---|
[1878] | 2666 | // find equality rows with effective nonzero rhs |
---|
[2385] | 2667 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2668 | if (abcUpper_[iRow] > abcLower_[iRow] || getInternalStatus(iRow) != basic) { |
---|
| 2669 | mustRowOut[iRow] = -2; |
---|
[1878] | 2670 | continue; |
---|
| 2671 | } |
---|
[2385] | 2672 | int chooseColumn[2] = { -1, -1 }; |
---|
| 2673 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 2674 | int iColumn = column[j]; |
---|
| 2675 | if (elementByRow[j] > 0.0) { |
---|
| 2676 | if (chooseColumn[0] == -1) |
---|
| 2677 | chooseColumn[0] = iColumn; |
---|
| 2678 | else |
---|
| 2679 | chooseColumn[0] = -2; |
---|
[1878] | 2680 | } else { |
---|
[2385] | 2681 | if (chooseColumn[1] == -1) |
---|
| 2682 | chooseColumn[1] = iColumn; |
---|
| 2683 | else |
---|
| 2684 | chooseColumn[1] = -2; |
---|
[1878] | 2685 | } |
---|
| 2686 | } |
---|
[2385] | 2687 | for (int iTry = 0; iTry < 2; iTry++) { |
---|
| 2688 | int jColumn = chooseColumn[iTry]; |
---|
| 2689 | if (jColumn >= 0 && getInternalStatus(jColumn) != basic) { |
---|
| 2690 | // see if has to be basic |
---|
| 2691 | double lowerValue = -abcUpper_[iRow]; // check swap |
---|
| 2692 | double upperValue = -abcLower_[iRow]; |
---|
| 2693 | int lowerInf = 0; |
---|
| 2694 | int upperInf = 0; |
---|
| 2695 | double alpha = 0.0; |
---|
| 2696 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 2697 | int iColumn = column[j]; |
---|
| 2698 | if (iColumn != jColumn) { |
---|
| 2699 | if (abcLower_[iColumn] > -1.0e20) |
---|
| 2700 | lowerValue -= abcLower_[iColumn] * elementByRow[j]; |
---|
| 2701 | else |
---|
| 2702 | lowerInf++; |
---|
| 2703 | if (abcUpper_[iColumn] < 1.0e20) |
---|
| 2704 | upperValue -= abcUpper_[iColumn] * elementByRow[j]; |
---|
| 2705 | else |
---|
| 2706 | upperInf++; |
---|
| 2707 | } else { |
---|
| 2708 | alpha = elementByRow[j]; |
---|
| 2709 | } |
---|
| 2710 | } |
---|
| 2711 | // find values column must lie between (signs again) |
---|
| 2712 | if (upperInf) |
---|
| 2713 | upperValue = COIN_DBL_MAX; |
---|
| 2714 | else |
---|
| 2715 | upperValue /= alpha; |
---|
| 2716 | if (lowerInf) |
---|
| 2717 | lowerValue = -COIN_DBL_MAX; |
---|
| 2718 | else |
---|
| 2719 | lowerValue /= alpha; |
---|
| 2720 | if (iTry) { |
---|
| 2721 | // swap |
---|
| 2722 | double temp = lowerValue; |
---|
| 2723 | lowerValue = upperValue; |
---|
| 2724 | upperValue = temp; |
---|
| 2725 | } |
---|
| 2726 | if (lowerValue > abcLower_[jColumn] + 10.0 * primalTolerance_ && upperValue < abcUpper_[jColumn] - 10.0 * primalTolerance_) { |
---|
| 2727 | nPossible++; |
---|
| 2728 | if (mustRowOut[jColumn] >= 0) { |
---|
| 2729 | // choose one ??? |
---|
| 2730 | //printf("Column %d already selected on row %d now on %d\n", |
---|
| 2731 | // jColumn,mustRowOut[jColumn],iRow); |
---|
| 2732 | continue; |
---|
| 2733 | } |
---|
| 2734 | mustRowOut[jColumn] = iRow; |
---|
| 2735 | mustRowOut[iRow] = jColumn; |
---|
| 2736 | } |
---|
[1878] | 2737 | } |
---|
| 2738 | } |
---|
| 2739 | } |
---|
| 2740 | if (nPossible) { |
---|
[2385] | 2741 | #if ABC_NORMAL_DEBUG > 0 |
---|
| 2742 | printf("%d possible candidates\n", nPossible); |
---|
[1878] | 2743 | #endif |
---|
[2385] | 2744 | if ((type & 2) == 0) { |
---|
[1878] | 2745 | // triangular |
---|
[2385] | 2746 | int iVector2 = getAvailableArray(); |
---|
| 2747 | int *COIN_RESTRICT counts = usefulArray_[iVector2].getIndices(); |
---|
| 2748 | CoinAbcMemset0(counts, numberRows_); |
---|
| 2749 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2750 | int n = 0; |
---|
| 2751 | for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) { |
---|
| 2752 | int iColumn = column[j]; |
---|
| 2753 | if (getInternalStatus(iColumn) == basic) |
---|
| 2754 | n++; |
---|
| 2755 | } |
---|
| 2756 | counts[iRow] = n; |
---|
[1878] | 2757 | } |
---|
[2385] | 2758 | const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() |
---|
| 2759 | - maximumAbcNumberRows_; |
---|
| 2760 | const int *row = abcMatrix_->getPackedMatrix()->getIndices(); |
---|
| 2761 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2762 | if (!counts[iRow]) { |
---|
| 2763 | int iColumn = mustRowOut[iRow]; |
---|
| 2764 | if (iColumn >= 0) { |
---|
| 2765 | setInternalStatus(iRow, isFixed); |
---|
| 2766 | solutionSaved_[iRow] = abcLower_[iRow]; |
---|
| 2767 | setInternalStatus(iColumn, basic); |
---|
| 2768 | for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) { |
---|
| 2769 | int jRow = row[j]; |
---|
| 2770 | counts[jRow]++; |
---|
| 2771 | } |
---|
| 2772 | } |
---|
| 2773 | } |
---|
[1878] | 2774 | } |
---|
| 2775 | setAvailableArray(iVector2); |
---|
| 2776 | } else { |
---|
| 2777 | // all |
---|
[2385] | 2778 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2779 | int iColumn = mustRowOut[iRow]; |
---|
| 2780 | if (iColumn >= 0) { |
---|
| 2781 | setInternalStatus(iRow, isFixed); |
---|
| 2782 | solutionSaved_[iRow] = abcLower_[iRow]; |
---|
| 2783 | setInternalStatus(iColumn, basic); |
---|
| 2784 | } |
---|
[1878] | 2785 | } |
---|
| 2786 | } |
---|
| 2787 | // redo pivot array |
---|
[2385] | 2788 | int numberBasic = 0; |
---|
| 2789 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
| 2790 | if (getInternalStatus(iSequence) == basic) |
---|
| 2791 | abcPivotVariable_[numberBasic++] = iSequence; |
---|
[1878] | 2792 | } |
---|
[2385] | 2793 | assert(numberBasic == numberRows_); |
---|
[1878] | 2794 | } |
---|
| 2795 | setAvailableArray(iVector); |
---|
| 2796 | } |
---|
| 2797 | // Computes nonbasic cost and total cost |
---|
[2385] | 2798 | void AbcSimplex::computeObjective() |
---|
[1878] | 2799 | { |
---|
| 2800 | sumNonBasicCosts_ = 0.0; |
---|
| 2801 | rawObjectiveValue_ = 0.0; |
---|
| 2802 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
[2385] | 2803 | double value = abcSolution_[iSequence] * abcCost_[iSequence]; |
---|
[1878] | 2804 | rawObjectiveValue_ += value; |
---|
[2385] | 2805 | if (getInternalStatus(iSequence) != basic) |
---|
[1878] | 2806 | sumNonBasicCosts_ += value; |
---|
| 2807 | } |
---|
| 2808 | setClpSimplexObjectiveValue(); |
---|
| 2809 | } |
---|
| 2810 | // This sets largest infeasibility and most infeasible |
---|
[2385] | 2811 | void AbcSimplex::checkPrimalSolution(bool justBasic) |
---|
[1878] | 2812 | { |
---|
[2385] | 2813 | rawObjectiveValue_ = CoinAbcInnerProduct(costBasic_, numberRows_, solutionBasic_); |
---|
| 2814 | //rawObjectiveValue_ += sumNonBasicCosts_; |
---|
[1878] | 2815 | setClpSimplexObjectiveValue(); |
---|
| 2816 | // now look at primal solution |
---|
| 2817 | sumPrimalInfeasibilities_ = 0.0; |
---|
| 2818 | numberPrimalInfeasibilities_ = 0; |
---|
| 2819 | double primalTolerance = primalTolerance_; |
---|
| 2820 | double relaxedTolerance = primalTolerance_; |
---|
| 2821 | // we can't really trust infeasibilities if there is primal error |
---|
[2385] | 2822 | double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 5.0 * primalTolerance_)); |
---|
[1878] | 2823 | // allow tolerance at least slightly bigger than standard |
---|
[2385] | 2824 | relaxedTolerance = relaxedTolerance + error; |
---|
[1878] | 2825 | sumOfRelaxedPrimalInfeasibilities_ = 0.0; |
---|
[2385] | 2826 | const double *COIN_RESTRICT lowerBasic = lowerBasic_; |
---|
| 2827 | const double *COIN_RESTRICT upperBasic = upperBasic_; |
---|
| 2828 | const double *COIN_RESTRICT solutionBasic = solutionBasic_; |
---|
[1878] | 2829 | if (justBasic) { |
---|
| 2830 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 2831 | double infeasibility = 0.0; |
---|
| 2832 | if (solutionBasic[iRow] > upperBasic[iRow]) { |
---|
[2385] | 2833 | infeasibility = solutionBasic[iRow] - upperBasic[iRow]; |
---|
[1878] | 2834 | } else if (solutionBasic[iRow] < lowerBasic[iRow]) { |
---|
[2385] | 2835 | infeasibility = lowerBasic[iRow] - solutionBasic[iRow]; |
---|
[1878] | 2836 | } |
---|
| 2837 | if (infeasibility > primalTolerance) { |
---|
[2385] | 2838 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
| 2839 | if (infeasibility > relaxedTolerance) |
---|
| 2840 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance; |
---|
| 2841 | numberPrimalInfeasibilities_++; |
---|
[1878] | 2842 | } |
---|
| 2843 | } |
---|
| 2844 | } else { |
---|
| 2845 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
| 2846 | double infeasibility = 0.0; |
---|
| 2847 | if (abcSolution_[iSequence] > abcUpper_[iSequence]) { |
---|
[2385] | 2848 | infeasibility = abcSolution_[iSequence] - abcUpper_[iSequence]; |
---|
[1878] | 2849 | } else if (abcSolution_[iSequence] < abcLower_[iSequence]) { |
---|
[2385] | 2850 | infeasibility = abcLower_[iSequence] - abcSolution_[iSequence]; |
---|
[1878] | 2851 | } |
---|
| 2852 | if (infeasibility > primalTolerance) { |
---|
[2385] | 2853 | //assert (getInternalStatus(iSequence)==basic); |
---|
| 2854 | sumPrimalInfeasibilities_ += infeasibility - primalTolerance_; |
---|
| 2855 | if (infeasibility > relaxedTolerance) |
---|
| 2856 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance; |
---|
| 2857 | numberPrimalInfeasibilities_++; |
---|
[1878] | 2858 | } |
---|
| 2859 | } |
---|
| 2860 | } |
---|
| 2861 | } |
---|
[2385] | 2862 | void AbcSimplex::checkDualSolution() |
---|
[1878] | 2863 | { |
---|
[2385] | 2864 | |
---|
[1878] | 2865 | sumDualInfeasibilities_ = 0.0; |
---|
| 2866 | numberDualInfeasibilities_ = 0; |
---|
| 2867 | int firstFreePrimal = -1; |
---|
| 2868 | int firstFreeDual = -1; |
---|
| 2869 | int numberSuperBasicWithDj = 0; |
---|
| 2870 | bestPossibleImprovement_ = 0.0; |
---|
| 2871 | // we can't really trust infeasibilities if there is dual error |
---|
[2385] | 2872 | double error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_)); |
---|
[1878] | 2873 | // allow tolerance at least slightly bigger than standard |
---|
[2385] | 2874 | double relaxedTolerance = currentDualTolerance_ + error; |
---|
[1878] | 2875 | // allow bigger tolerance for possible improvement |
---|
| 2876 | double possTolerance = 5.0 * relaxedTolerance; |
---|
| 2877 | sumOfRelaxedDualInfeasibilities_ = 0.0; |
---|
[2385] | 2878 | int numberNeedToMove = 0; |
---|
[1878] | 2879 | for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { |
---|
| 2880 | if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) { |
---|
| 2881 | // not basic |
---|
| 2882 | double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence]; |
---|
| 2883 | double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence]; |
---|
| 2884 | double value = abcDj_[iSequence]; |
---|
| 2885 | if (distanceUp > primalTolerance_) { |
---|
[2385] | 2886 | // Check if "free" |
---|
| 2887 | if (distanceDown <= primalTolerance_) { |
---|
| 2888 | // should not be negative |
---|
| 2889 | if (value < 0.0) { |
---|
| 2890 | value = -value; |
---|
| 2891 | if (value > currentDualTolerance_) { |
---|
| 2892 | sumDualInfeasibilities_ += value - currentDualTolerance_; |
---|
| 2893 | if (value > possTolerance) |
---|
| 2894 | bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value; |
---|
| 2895 | if (value > relaxedTolerance) |
---|
| 2896 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
| 2897 | numberDualInfeasibilities_++; |
---|
| 2898 | } |
---|
| 2899 | } |
---|
| 2900 | } else { |
---|
| 2901 | // free |
---|
| 2902 | value = fabs(value); |
---|
| 2903 | if (value > 1.0 * relaxedTolerance) { |
---|
| 2904 | numberSuperBasicWithDj++; |
---|
| 2905 | if (firstFreeDual < 0) |
---|
| 2906 | firstFreeDual = iSequence; |
---|
| 2907 | } |
---|
| 2908 | if (value > currentDualTolerance_ || getInternalStatus(iSequence) != AbcSimplex::isFree) { |
---|
| 2909 | numberNeedToMove++; |
---|
| 2910 | if (firstFreePrimal < 0) |
---|
| 2911 | firstFreePrimal = iSequence; |
---|
| 2912 | sumDualInfeasibilities_ += value - currentDualTolerance_; |
---|
| 2913 | if (value > possTolerance) |
---|
| 2914 | bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value; |
---|
| 2915 | if (value > relaxedTolerance) |
---|
| 2916 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
| 2917 | numberDualInfeasibilities_++; |
---|
| 2918 | } |
---|
| 2919 | } |
---|
[1878] | 2920 | } else if (distanceDown > primalTolerance_) { |
---|
[2385] | 2921 | // should not be positive |
---|
| 2922 | if (value > 0.0) { |
---|
| 2923 | if (value > currentDualTolerance_) { |
---|
| 2924 | sumDualInfeasibilities_ += value - currentDualTolerance_; |
---|
| 2925 | if (value > possTolerance) |
---|
| 2926 | bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10); |
---|
| 2927 | if (value > relaxedTolerance) |
---|
| 2928 | sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance; |
---|
| 2929 | numberDualInfeasibilities_++; |
---|
| 2930 | } |
---|
| 2931 | } |
---|
[1878] | 2932 | } |
---|
| 2933 | } |
---|
| 2934 | } |
---|
[2385] | 2935 | numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj; |
---|
[1878] | 2936 | if (algorithm_ < 0 && firstFreeDual >= 0) { |
---|
| 2937 | // dual |
---|
| 2938 | firstFree_ = firstFreeDual; |
---|
[2385] | 2939 | } else if (numberSuperBasicWithDj || numberNeedToMove) { |
---|
[1878] | 2940 | //(abcProgress_.lastIterationNumber(0) <= 0)) { |
---|
| 2941 | firstFree_ = firstFreePrimal; |
---|
| 2942 | if (!sumDualInfeasibilities_) |
---|
[2385] | 2943 | sumDualInfeasibilities_ = 1.0e-8; |
---|
[1878] | 2944 | } |
---|
| 2945 | } |
---|
| 2946 | /* This sets sum and number of infeasibilities (Dual and Primal) */ |
---|
[2385] | 2947 | void AbcSimplex::checkBothSolutions() |
---|
[1878] | 2948 | { |
---|
[2385] | 2949 | // old way |
---|
[1878] | 2950 | checkPrimalSolution(true); |
---|
| 2951 | checkDualSolution(); |
---|
| 2952 | } |
---|
| 2953 | /* |
---|
| 2954 | Unpacks one column of the matrix into indexed array |
---|
| 2955 | */ |
---|
[2385] | 2956 | void AbcSimplex::unpack(CoinIndexedVector &rowArray, int sequence) const |
---|
| 2957 | { |
---|
| 2958 | abcMatrix_->unpack(rowArray, sequence); |
---|
[1878] | 2959 | } |
---|
| 2960 | // Sets row pivot choice algorithm in dual |
---|
[2385] | 2961 | void AbcSimplex::setDualRowPivotAlgorithm(AbcDualRowPivot &choice) |
---|
[1878] | 2962 | { |
---|
| 2963 | delete abcDualRowPivot_; |
---|
| 2964 | abcDualRowPivot_ = choice.clone(true); |
---|
| 2965 | abcDualRowPivot_->setModel(this); |
---|
| 2966 | } |
---|
| 2967 | // Sets column pivot choice algorithm in primal |
---|
[2385] | 2968 | void AbcSimplex::setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot &choice) |
---|
[1878] | 2969 | { |
---|
| 2970 | delete abcPrimalColumnPivot_; |
---|
| 2971 | abcPrimalColumnPivot_ = choice.clone(true); |
---|
| 2972 | abcPrimalColumnPivot_->setModel(this); |
---|
| 2973 | } |
---|
[2385] | 2974 | void AbcSimplex::setFactorization(AbcSimplexFactorization &factorization) |
---|
[1878] | 2975 | { |
---|
| 2976 | if (abcFactorization_) |
---|
| 2977 | abcFactorization_->setFactorization(factorization); |
---|
| 2978 | else |
---|
| 2979 | abcFactorization_ = new AbcSimplexFactorization(factorization, |
---|
[2385] | 2980 | numberRows_); |
---|
[1878] | 2981 | } |
---|
| 2982 | |
---|
| 2983 | // Swaps factorization |
---|
| 2984 | AbcSimplexFactorization * |
---|
[2385] | 2985 | AbcSimplex::swapFactorization(AbcSimplexFactorization *factorization) |
---|
[1878] | 2986 | { |
---|
[2385] | 2987 | AbcSimplexFactorization *swap = abcFactorization_; |
---|
[1878] | 2988 | abcFactorization_ = factorization; |
---|
| 2989 | return swap; |
---|
| 2990 | } |
---|
| 2991 | |
---|
| 2992 | /* Tightens primal bounds to make dual faster. Unless |
---|
| 2993 | fixed, bounds are slightly looser than they could be. |
---|
| 2994 | This is to make dual go faster and is probably not needed |
---|
| 2995 | with a presolve. Returns non-zero if problem infeasible |
---|
| 2996 | */ |
---|
[2385] | 2997 | int AbcSimplex::tightenPrimalBounds() |
---|
[1878] | 2998 | { |
---|
[2385] | 2999 | int tightenType = 1; |
---|
| 3000 | if (maximumIterations() >= 1000000 && maximumIterations() < 1000010) |
---|
| 3001 | tightenType = maximumIterations() - 1000000; |
---|
[1878] | 3002 | if (!tightenType) |
---|
| 3003 | return 0; |
---|
| 3004 | if (integerType_) { |
---|
[2385] | 3005 | #if ABC_NORMAL_DEBUG > 0 |
---|
[1878] | 3006 | printf("Redo tighten to relax by 1 for integers (and don't be shocked by infeasibility)\n"); |
---|
| 3007 | #endif |
---|
| 3008 | return 0; |
---|
| 3009 | } |
---|
[2385] | 3010 | // This needs to work on scaled matrix - then replace |
---|
[1878] | 3011 | // Get a row copy in standard format |
---|
| 3012 | CoinPackedMatrix copy; |
---|
| 3013 | copy.setExtraGap(0.0); |
---|
| 3014 | copy.setExtraMajor(0.0); |
---|
| 3015 | copy.reverseOrderedCopyOf(*abcMatrix_->matrix()); |
---|
| 3016 | // get matrix data pointers |
---|
[2385] | 3017 | const int *COIN_RESTRICT column = copy.getIndices(); |
---|
| 3018 | const CoinBigIndex *COIN_RESTRICT rowStart = copy.getVectorStarts(); |
---|
| 3019 | const int *COIN_RESTRICT rowLength = copy.getVectorLengths(); |
---|
| 3020 | const double *COIN_RESTRICT element = copy.getElements(); |
---|
[1878] | 3021 | int numberChanged = 1, iPass = 0; |
---|
| 3022 | double large = largeValue(); // treat bounds > this as infinite |
---|
| 3023 | #ifndef NDEBUG |
---|
| 3024 | double large2 = 1.0e10 * large; |
---|
| 3025 | #endif |
---|
| 3026 | int numberInfeasible = 0; |
---|
| 3027 | int totalTightened = 0; |
---|
[2385] | 3028 | |
---|
[1878] | 3029 | double tolerance = primalTolerance(); |
---|
[2385] | 3030 | |
---|
[1878] | 3031 | // A copy of bounds is up at top |
---|
| 3032 | translate(0); // move down |
---|
[2385] | 3033 | |
---|
[1878] | 3034 | #define MAXPASS 10 |
---|
[2385] | 3035 | |
---|
[1878] | 3036 | // Loop round seeing if we can tighten bounds |
---|
| 3037 | // Would be faster to have a stack of possible rows |
---|
| 3038 | // and we put altered rows back on stack |
---|
| 3039 | int numberCheck = -1; |
---|
| 3040 | // temp swap signs so can use existing code |
---|
[2385] | 3041 | double *COIN_RESTRICT rowLower = abcLower_; |
---|
| 3042 | double *COIN_RESTRICT rowUpper = abcUpper_; |
---|
| 3043 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
| 3044 | double value = -rowLower[iRow]; |
---|
| 3045 | rowLower[iRow] = -rowUpper[iRow]; |
---|
| 3046 | rowUpper[iRow] = value; |
---|
[1878] | 3047 | } |
---|
| 3048 | // Keep lower and upper for rows |
---|
| 3049 | int whichArray[5]; |
---|
[2385] | 3050 | for (int i = 0; i < 5; i++) |
---|
| 3051 | whichArray[i] = getAvailableArray(); |
---|
| 3052 | double *COIN_RESTRICT minRhs = usefulArray_[whichArray[0]].denseVector(); |
---|
| 3053 | double *COIN_RESTRICT maxRhs = usefulArray_[whichArray[1]].denseVector(); |
---|
| 3054 | int *COIN_RESTRICT changedList = usefulArray_[whichArray[0]].getIndices(); |
---|
| 3055 | int *COIN_RESTRICT changedListNext = usefulArray_[whichArray[1]].getIndices(); |
---|
| 3056 | char *COIN_RESTRICT changed = reinterpret_cast< char * >(usefulArray_[whichArray[2]].getIndices()); |
---|
[1878] | 3057 | // -1 no infinite, -2 more than one infinite >=0 column |
---|
[2385] | 3058 | int *COIN_RESTRICT whichInfiniteDown = usefulArray_[whichArray[3]].getIndices(); |
---|
| 3059 | int *COIN_RESTRICT whichInfiniteUp = usefulArray_[whichArray[3]].getIndices(); |
---|
| 3060 | int numberToDoNext = 0; |
---|
| 3061 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
[1878] | 3062 | if (rowLower[iRow] > -large || rowUpper[iRow] < large) { |
---|
[2385] | 3063 | changedListNext[numberToDoNext++] = iRow; |
---|
| 3064 | whichInfiniteDown[iRow] = -1; |
---|
| 3065 | whichInfiniteUp[iRow] = -1; |
---|
[1878] | 3066 | } else { |
---|
[2385] | 3067 | minRhs[iRow] = -COIN_DBL_MAX; |
---|
| 3068 | maxRhs[iRow] = COIN_DBL_MAX; |
---|
| 3069 | whichInfiniteDown[iRow] = -2; |
---|
| 3070 | whichInfiniteUp[iRow] = -2; |
---|
[1878] | 3071 | } |
---|
| 3072 | } |
---|
[2385] | 3073 | const int *COIN_RESTRICT row = abcMatrix_->matrix()->getIndices(); |
---|
| 3074 | const CoinBigIndex *COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts(); |
---|
| 3075 | const double *COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements(); |
---|
| 3076 | |
---|
| 3077 | double *COIN_RESTRICT columnLower = abcLower_ + maximumAbcNumberRows_; |
---|
| 3078 | double *COIN_RESTRICT columnUpper = abcUpper_ + maximumAbcNumberRows_; |
---|
| 3079 | while (numberChanged > numberCheck) { |
---|
| 3080 | int numberToDo = numberToDoNext; |
---|
| 3081 | numberToDoNext = 0; |
---|
| 3082 | int *COIN_RESTRICT temp = changedList; |
---|
| 3083 | changedList = changedListNext; |
---|
| 3084 | changedListNext = temp; |
---|
| 3085 | CoinAbcMemset0(changed, numberRows_); |
---|
| 3086 | |
---|
[1878] | 3087 | numberChanged = 0; // Bounds tightened this pass |
---|
[2385] | 3088 | |
---|
| 3089 | if (iPass == MAXPASS) |
---|
| 3090 | break; |
---|
[1878] | 3091 | iPass++; |
---|
[2385] | 3092 | for (int k = 0; k < numberToDo; k++) { |
---|
| 3093 | int iRow = changedList[k]; |
---|
[1878] | 3094 | // possible row |
---|
| 3095 | int infiniteUpper = 0; |
---|
| 3096 | int infiniteLower = 0; |
---|
[2385] | 3097 | int firstInfiniteUpper = -1; |
---|
| 3098 | int firstInfiniteLower = -1; |
---|
[1878] | 3099 | double maximumUp = 0.0; |
---|
| 3100 | double maximumDown = 0.0; |
---|
| 3101 | double newBound; |
---|
| 3102 | CoinBigIndex rStart = rowStart[iRow]; |
---|
| 3103 | CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow]; |
---|
| 3104 | CoinBigIndex j; |
---|
| 3105 | // Compute possible lower and upper ranges |
---|
[2385] | 3106 | |
---|
[1878] | 3107 | for (j = rStart; j < rEnd; ++j) { |
---|
[2385] | 3108 | double value = element[j]; |
---|
| 3109 | int iColumn = column[j]; |
---|
| 3110 | if (value > 0.0) { |
---|
| 3111 | if (columnUpper[iColumn] >= large) { |
---|
| 3112 | firstInfiniteUpper = iColumn; |
---|
| 3113 | ++infiniteUpper; |
---|
| 3114 | } else { |
---|
| 3115 | maximumUp += columnUpper[iColumn] * value; |
---|
| 3116 | } |
---|
| 3117 | if (columnLower[iColumn] <= -large) { |
---|
| 3118 | firstInfiniteLower = iColumn; |
---|
| 3119 | ++infiniteLower; |
---|
| 3120 | } else { |
---|
| 3121 | maximumDown += columnLower[iColumn] * value; |
---|
| 3122 | } |
---|
| 3123 | } else if (value < 0.0) { |
---|
| 3124 | if (columnUpper[iColumn] >= large) { |
---|
| 3125 | firstInfiniteLower = iColumn; |
---|
| 3126 | ++infiniteLower; |
---|
| 3127 | } else { |
---|
| 3128 | maximumDown += columnUpper[iColumn] * value; |
---|
| 3129 | } |
---|
| 3130 | if (columnLower[iColumn] <= -large) { |
---|
| 3131 | firstInfiniteUpper = iColumn; |
---|
| 3132 | ++infiniteUpper; |
---|
| 3133 | } else { |
---|
| 3134 | maximumUp += columnLower[iColumn] * value; |
---|
| 3135 | } |
---|
| 3136 | } |
---|
[1878] | 3137 | } |
---|
| 3138 | // Build in a margin of error |
---|
| 3139 | maximumUp += 1.0e-8 * fabs(maximumUp); |
---|
| 3140 | maximumDown -= 1.0e-8 * fabs(maximumDown); |
---|
| 3141 | double maxUp = maximumUp + infiniteUpper * 1.0e31; |
---|
| 3142 | double maxDown = maximumDown - infiniteLower * 1.0e31; |
---|
[2385] | 3143 | minRhs[iRow] = infiniteLower ? -COIN_DBL_MAX : maximumDown; |
---|
| 3144 | maxRhs[iRow] = infiniteUpper ? COIN_DBL_MAX : maximumUp; |
---|
| 3145 | if (infiniteLower == 0) |
---|
| 3146 | whichInfiniteDown[iRow] = -1; |
---|
| 3147 | else if (infiniteLower == 1) |
---|
| 3148 | whichInfiniteDown[iRow] = firstInfiniteLower; |
---|
[1878] | 3149 | else |
---|
[2385] | 3150 | whichInfiniteDown[iRow] = -2; |
---|
| 3151 | if (infiniteUpper == 0) |
---|
| 3152 | whichInfiniteUp[iRow] = -1; |
---|
| 3153 | else if (infiniteUpper == 1) |
---|
| 3154 | whichInfiniteUp[iRow] = firstInfiniteUpper; |
---|
[1878] | 3155 | else |
---|
[2385] | 3156 | whichInfiniteUp[iRow] = -2; |
---|
| 3157 | if (maxUp <= rowUpper[iRow] + tolerance && maxDown >= rowLower[iRow] - tolerance) { |
---|
| 3158 | |
---|
| 3159 | // Row is redundant - make totally free |
---|
| 3160 | // NO - can't do this for postsolve |
---|
| 3161 | // rowLower[iRow]=-COIN_DBL_MAX; |
---|
| 3162 | // rowUpper[iRow]=COIN_DBL_MAX; |
---|
| 3163 | //printf("Redundant row in presolveX %d\n",iRow); |
---|
| 3164 | |
---|
[1878] | 3165 | } else { |
---|
[2385] | 3166 | if (maxUp < rowLower[iRow] - 100.0 * tolerance || maxDown > rowUpper[iRow] + 100.0 * tolerance) { |
---|
| 3167 | // problem is infeasible - exit at once |
---|
| 3168 | numberInfeasible++; |
---|
| 3169 | break; |
---|
| 3170 | } |
---|
| 3171 | double lower = rowLower[iRow]; |
---|
| 3172 | double upper = rowUpper[iRow]; |
---|
| 3173 | for (j = rStart; j < rEnd; ++j) { |
---|
| 3174 | double value = element[j]; |
---|
| 3175 | int iColumn = column[j]; |
---|
| 3176 | double nowLower = columnLower[iColumn]; |
---|
| 3177 | double nowUpper = columnUpper[iColumn]; |
---|
| 3178 | if (value > 0.0) { |
---|
| 3179 | // positive value |
---|
| 3180 | if (lower > -large) { |
---|
| 3181 | if (!infiniteUpper) { |
---|
| 3182 | assert(nowUpper < large2); |
---|
| 3183 | newBound = nowUpper + (lower - maximumUp) / value; |
---|
| 3184 | // relax if original was large |
---|
| 3185 | if (fabs(maximumUp) > 1.0e8) |
---|
| 3186 | newBound -= 1.0e-12 * fabs(maximumUp); |
---|
| 3187 | } else if (infiniteUpper == 1 && nowUpper > large) { |
---|
| 3188 | newBound = (lower - maximumUp) / value; |
---|
| 3189 | // relax if original was large |
---|
| 3190 | if (fabs(maximumUp) > 1.0e8) |
---|
| 3191 | newBound -= 1.0e-12 * fabs(maximumUp); |
---|
| 3192 | } else { |
---|
| 3193 | newBound = -COIN_DBL_MAX; |
---|
| 3194 | } |
---|
| 3195 | if (newBound > nowLower + 1.0e-12 && newBound > -large) { |
---|
| 3196 | // Tighten the lower bound |
---|
| 3197 | numberChanged++; |
---|
| 3198 | // check infeasible (relaxed) |
---|
| 3199 | if (nowUpper < newBound) { |
---|
| 3200 | if (nowUpper - newBound < -100.0 * tolerance) |
---|
| 3201 | numberInfeasible++; |
---|
| 3202 | else |
---|
| 3203 | newBound = nowUpper; |
---|
| 3204 | } |
---|
| 3205 | columnLower[iColumn] = newBound; |
---|
| 3206 | for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) { |
---|
| 3207 | int jRow = row[j1]; |
---|
| 3208 | if (!changed[jRow]) { |
---|
| 3209 | changed[jRow] = 1; |
---|
| 3210 | changedListNext[numberToDoNext++] = jRow; |
---|
| 3211 | } |
---|
| 3212 | } |
---|
| 3213 | // adjust |
---|
| 3214 | double now; |
---|
| 3215 | if (nowLower < -large) { |
---|
| 3216 | now = 0.0; |
---|
| 3217 | infiniteLower--; |
---|
| 3218 | } else { |
---|
| 3219 | now = nowLower; |
---|
| 3220 | } |
---|
| 3221 | maximumDown += (newBound - now) * value; |
---|
| 3222 | nowLower = newBound; |
---|
| 3223 | } |
---|
| 3224 | } |
---|
| 3225 | if (upper < large) { |
---|
| 3226 | if (!infiniteLower) { |
---|
| 3227 | assert(nowLower > -large2); |
---|
| 3228 | newBound = nowLower + (upper - maximumDown) / value; |
---|
| 3229 | // relax if original was large |
---|
| 3230 | if (fabs(maximumDown) > 1.0e8) |
---|
| 3231 | newBound += 1.0e-12 * fabs(maximumDown); |
---|
| 3232 | } else if (infiniteLower == 1 && nowLower < -large) { |
---|
| 3233 | newBound = (upper - maximumDown) / value; |
---|
| 3234 | // relax if original was large |
---|
| 3235 | if (fabs(maximumDown) > 1.0e8) |
---|
| 3236 | newBound += 1.0e-12 * fabs(maximumDown); |
---|
| 3237 | } else { |
---|
| 3238 | newBound = COIN_DBL_MAX; |
---|
| 3239 | } |
---|
| 3240 | if (newBound < nowUpper - 1.0e-12 && newBound < large) { |
---|
| 3241 | // Tighten the upper bound |
---|
| 3242 | numberChanged++; |
---|
| 3243 | // check infeasible (relaxed) |
---|
| 3244 | if (nowLower > newBound) { |
---|
| 3245 | if (newBound - nowLower < -100.0 * tolerance) |
---|
| 3246 | numberInfeasible++; |
---|
| 3247 | else |
---|
| 3248 | newBound = nowLower; |
---|
| 3249 | } |
---|
| 3250 | columnUpper[iColumn] = newBound; |
---|
| 3251 | for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) { |
---|
| 3252 | int jRow = row[j1]; |
---|
| 3253 | if (!changed[jRow]) { |
---|
| 3254 | changed[jRow] = 1; |
---|
| 3255 | changedListNext[numberToDoNext++] = jRow; |
---|
| 3256 | } |
---|
| 3257 | } |
---|
| 3258 | // adjust |
---|
| 3259 | double now; |
---|
| 3260 | if (nowUpper > large) { |
---|
| 3261 | now = 0.0; |
---|
| 3262 | infiniteUpper--; |
---|
| 3263 | } else { |
---|
| 3264 | now = nowUpper; |
---|
| 3265 | } |
---|
| 3266 | maximumUp += (newBound - now) * value; |
---|
| 3267 | nowUpper = newBound; |
---|
| 3268 | } |
---|
| 3269 | } |
---|
| 3270 | } else { |
---|
| 3271 | // negative value |
---|
| 3272 | if (lower > -large) { |
---|
| 3273 | if (!infiniteUpper) { |
---|
| 3274 | assert(nowLower < large2); |
---|
| 3275 | newBound = nowLower + (lower - maximumUp) / value; |
---|
| 3276 | // relax if original was large |
---|
| 3277 | if (fabs(maximumUp) > 1.0e8) |
---|
| 3278 | newBound += 1.0e-12 * fabs(maximumUp); |
---|
| 3279 | } else if (infiniteUpper == 1 && nowLower < -large) { |
---|
| 3280 | newBound = (lower - maximumUp) / value; |
---|
| 3281 | // relax if original was large |
---|
| 3282 | if (fabs(maximumUp) > 1.0e8) |
---|
| 3283 | newBound += 1.0e-12 * fabs(maximumUp); |
---|
| 3284 | } else { |
---|
| 3285 | newBound = COIN_DBL_MAX; |
---|
| 3286 | } |
---|
| 3287 | if (newBound < nowUpper - 1.0e-12 && newBound < large) { |
---|
| 3288 | // Tighten the upper bound |
---|
| 3289 | numberChanged++; |
---|
| 3290 | // check infeasible (relaxed) |
---|
| 3291 | if (nowLower > newBound) { |
---|
| 3292 | if (newBound - nowLower < -100.0 * tolerance) |
---|
| 3293 | numberInfeasible++; |
---|
| 3294 | else |
---|
| 3295 | newBound = nowLower; |
---|
| 3296 | } |
---|
| 3297 | columnUpper[iColumn] = newBound; |
---|
| 3298 | for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) { |
---|
| 3299 | int jRow = row[j1]; |
---|
| 3300 | if (!changed[jRow]) { |
---|
| 3301 | changed[jRow] = 1; |
---|
| 3302 | changedListNext[numberToDoNext++] = jRow; |
---|
| 3303 | } |
---|
| 3304 | } |
---|
| 3305 | // adjust |
---|
| 3306 | double now; |
---|
| 3307 | if (nowUpper > large) { |
---|
| 3308 | now = 0.0; |
---|
| 3309 | infiniteLower--; |
---|
| 3310 | } else { |
---|
| 3311 | now = nowUpper; |
---|
| 3312 | } |
---|
| 3313 | maximumDown += (newBound |
---|