[328] | 1 | // Copyright (C) 2004, International Business Machines |
---|
| 2 | // Corporation and others. All Rights Reserved. |
---|
| 3 | |
---|
| 4 | |
---|
| 5 | #include <cstdio> |
---|
| 6 | |
---|
| 7 | #include "CoinPragma.hpp" |
---|
| 8 | #include "CoinIndexedVector.hpp" |
---|
| 9 | #include "CoinHelperFunctions.hpp" |
---|
| 10 | |
---|
| 11 | #include "ClpSimplex.hpp" |
---|
| 12 | #include "ClpFactorization.hpp" |
---|
| 13 | #include "ClpQuadraticObjective.hpp" |
---|
| 14 | #include "ClpNonLinearCost.hpp" |
---|
| 15 | // at end to get min/max! |
---|
| 16 | #include "ClpDynamicMatrix.hpp" |
---|
| 17 | #include "ClpMessage.hpp" |
---|
[335] | 18 | //#define CLP_DEBUG |
---|
[328] | 19 | //#define CLP_DEBUG_PRINT |
---|
| 20 | //############################################################################# |
---|
| 21 | // Constructors / Destructor / Assignment |
---|
| 22 | //############################################################################# |
---|
| 23 | |
---|
| 24 | //------------------------------------------------------------------- |
---|
| 25 | // Default Constructor |
---|
| 26 | //------------------------------------------------------------------- |
---|
| 27 | ClpDynamicMatrix::ClpDynamicMatrix () |
---|
| 28 | : ClpPackedMatrix(), |
---|
| 29 | sumDualInfeasibilities_(0.0), |
---|
| 30 | sumPrimalInfeasibilities_(0.0), |
---|
| 31 | sumOfRelaxedDualInfeasibilities_(0.0), |
---|
| 32 | sumOfRelaxedPrimalInfeasibilities_(0.0), |
---|
| 33 | savedBestGubDual_(0.0), |
---|
| 34 | savedBestSet_(0), |
---|
| 35 | backToPivotRow_(NULL), |
---|
| 36 | keyVariable_(NULL), |
---|
| 37 | toIndex_(NULL), |
---|
| 38 | fromIndex_(NULL), |
---|
| 39 | numberSets_(0), |
---|
| 40 | numberActiveSets_(0), |
---|
| 41 | objectiveOffset_(0.0), |
---|
| 42 | lowerSet_(NULL), |
---|
| 43 | upperSet_(NULL), |
---|
| 44 | status_(NULL), |
---|
| 45 | model_(NULL), |
---|
| 46 | firstAvailable_(0), |
---|
| 47 | firstAvailableBefore_(0), |
---|
| 48 | firstDynamic_(0), |
---|
| 49 | lastDynamic_(0), |
---|
| 50 | numberStaticRows_(0), |
---|
| 51 | numberElements_(0), |
---|
| 52 | numberDualInfeasibilities_(0), |
---|
| 53 | numberPrimalInfeasibilities_(0), |
---|
| 54 | noCheck_(-1), |
---|
| 55 | infeasibilityWeight_(0.0), |
---|
| 56 | numberGubColumns_(0), |
---|
| 57 | maximumGubColumns_(0), |
---|
| 58 | maximumElements_(0), |
---|
| 59 | startSet_(NULL), |
---|
| 60 | next_(NULL), |
---|
| 61 | startColumn_(NULL), |
---|
| 62 | row_(NULL), |
---|
| 63 | element_(NULL), |
---|
| 64 | cost_(NULL), |
---|
| 65 | id_(NULL), |
---|
| 66 | dynamicStatus_(NULL), |
---|
| 67 | columnLower_(NULL), |
---|
| 68 | columnUpper_(NULL) |
---|
| 69 | { |
---|
| 70 | setType(15); |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | //------------------------------------------------------------------- |
---|
| 74 | // Copy constructor |
---|
| 75 | //------------------------------------------------------------------- |
---|
| 76 | ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs) |
---|
| 77 | : ClpPackedMatrix(rhs) |
---|
| 78 | { |
---|
| 79 | objectiveOffset_ = rhs.objectiveOffset_; |
---|
| 80 | numberSets_ = rhs.numberSets_; |
---|
| 81 | numberActiveSets_ = rhs.numberActiveSets_; |
---|
| 82 | firstAvailable_ = rhs.firstAvailable_; |
---|
| 83 | firstAvailableBefore_ = rhs.firstAvailableBefore_; |
---|
| 84 | firstDynamic_ = rhs.firstDynamic_; |
---|
| 85 | lastDynamic_ = rhs.lastDynamic_; |
---|
| 86 | numberStaticRows_ = rhs.numberStaticRows_; |
---|
| 87 | numberElements_ = rhs.numberElements_; |
---|
| 88 | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_); |
---|
| 89 | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_); |
---|
| 90 | toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_); |
---|
| 91 | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_); |
---|
| 92 | lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_); |
---|
| 93 | upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_); |
---|
| 94 | status_ = ClpCopyOfArray(rhs.status_,numberSets_); |
---|
| 95 | model_ = rhs.model_; |
---|
| 96 | sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; |
---|
| 97 | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
---|
| 98 | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
---|
| 99 | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
---|
| 100 | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
---|
| 101 | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
---|
| 102 | savedBestGubDual_=rhs.savedBestGubDual_; |
---|
| 103 | savedBestSet_=rhs.savedBestSet_; |
---|
| 104 | noCheck_ = rhs.noCheck_; |
---|
| 105 | infeasibilityWeight_=rhs.infeasibilityWeight_; |
---|
| 106 | // Now secondary data |
---|
| 107 | numberGubColumns_ = rhs.numberGubColumns_; |
---|
| 108 | maximumGubColumns_=rhs.maximumGubColumns_; |
---|
| 109 | maximumElements_ = rhs.maximumElements_; |
---|
[330] | 110 | startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_); |
---|
[328] | 111 | next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_); |
---|
[330] | 112 | startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1); |
---|
| 113 | row_ = ClpCopyOfArray(rhs.row_,maximumElements_); |
---|
| 114 | element_ = ClpCopyOfArray(rhs.element_,maximumElements_); |
---|
| 115 | cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_); |
---|
[328] | 116 | id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_); |
---|
[330] | 117 | columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_); |
---|
| 118 | columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_); |
---|
| 119 | dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_); |
---|
[328] | 120 | } |
---|
| 121 | |
---|
| 122 | /* This is the real constructor*/ |
---|
| 123 | ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * model, int numberSets, |
---|
| 124 | int numberGubColumns, const int * starts, |
---|
| 125 | const double * lower, const double * upper, |
---|
| 126 | const CoinBigIndex * startColumn, const int * row, |
---|
| 127 | const double * element, const double * cost, |
---|
| 128 | const double * columnLower, const double * columnUpper, |
---|
| 129 | const unsigned char * status, |
---|
| 130 | const unsigned char * dynamicStatus) |
---|
| 131 | : ClpPackedMatrix() |
---|
| 132 | { |
---|
| 133 | setType(15); |
---|
| 134 | objectiveOffset_ = model->objectiveOffset(); |
---|
| 135 | model_ = model; |
---|
| 136 | numberSets_ = numberSets; |
---|
| 137 | numberGubColumns_ = numberGubColumns; |
---|
| 138 | maximumGubColumns_=numberGubColumns_; |
---|
[330] | 139 | if (numberGubColumns_) |
---|
| 140 | maximumElements_ = startColumn[numberGubColumns_]; |
---|
| 141 | else |
---|
| 142 | maximumElements_ = 0; |
---|
[328] | 143 | startSet_ = new int [numberSets_]; |
---|
| 144 | next_ = new int [maximumGubColumns_]; |
---|
| 145 | // fill in startSet and next |
---|
| 146 | int iSet; |
---|
[330] | 147 | if (numberGubColumns_) { |
---|
| 148 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 149 | int first = starts[iSet]; |
---|
| 150 | int last = starts[iSet+1]-1; |
---|
| 151 | startSet_[iSet]=first; |
---|
| 152 | for (int i=first;i<last;i++) |
---|
| 153 | next_[i]=i+1; |
---|
| 154 | next_[last]=-iSet-1; |
---|
| 155 | } |
---|
[328] | 156 | } |
---|
| 157 | int numberColumns = model->numberColumns(); |
---|
| 158 | int numberRows = model->numberRows(); |
---|
| 159 | numberStaticRows_ = numberRows; |
---|
| 160 | savedBestGubDual_=0.0; |
---|
| 161 | savedBestSet_=0; |
---|
| 162 | // Number of columns needed |
---|
| 163 | int frequency = model->factorizationFrequency(); |
---|
[399] | 164 | int numberGubInSmall = numberRows + frequency + CoinMin(frequency,numberSets_)+4; |
---|
[328] | 165 | // for small problems this could be too big |
---|
[399] | 166 | //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_); |
---|
[328] | 167 | int numberNeeded = numberGubInSmall + numberColumns; |
---|
| 168 | firstAvailable_ = numberColumns; |
---|
| 169 | firstAvailableBefore_ = firstAvailable_; |
---|
| 170 | firstDynamic_ = numberColumns; |
---|
| 171 | lastDynamic_ = numberNeeded; |
---|
| 172 | startColumn_ = ClpCopyOfArray(startColumn,numberGubColumns_+1); |
---|
[330] | 173 | if (!numberGubColumns_) { |
---|
| 174 | startColumn_ = new CoinBigIndex [1]; |
---|
| 175 | startColumn_[0]=0; |
---|
| 176 | } |
---|
[328] | 177 | CoinBigIndex numberElements = startColumn_[numberGubColumns_]; |
---|
| 178 | row_ = ClpCopyOfArray(row,numberElements); |
---|
| 179 | element_ = new float[numberElements]; |
---|
| 180 | CoinBigIndex i; |
---|
| 181 | for (i=0;i<numberElements;i++) |
---|
| 182 | element_[i]=element[i]; |
---|
| 183 | cost_ = new float[numberGubColumns_]; |
---|
| 184 | for (i=0;i<numberGubColumns_;i++) { |
---|
| 185 | cost_[i]=cost[i]; |
---|
| 186 | // I don't think I need sorted but ... |
---|
| 187 | CoinSort_2(row_+startColumn_[i],row_+startColumn_[i+1],element_+startColumn_[i]); |
---|
| 188 | } |
---|
| 189 | if (columnLower) { |
---|
| 190 | columnLower_ = new float[numberGubColumns_]; |
---|
| 191 | for (i=0;i<numberGubColumns_;i++) |
---|
| 192 | columnLower_[i]=columnLower[i]; |
---|
| 193 | } else { |
---|
| 194 | columnLower_=NULL; |
---|
| 195 | } |
---|
| 196 | if (columnUpper) { |
---|
| 197 | columnUpper_ = new float[numberGubColumns_]; |
---|
| 198 | for (i=0;i<numberGubColumns_;i++) |
---|
| 199 | columnUpper_[i]=columnUpper[i]; |
---|
| 200 | } else { |
---|
| 201 | columnUpper_=NULL; |
---|
| 202 | } |
---|
| 203 | lowerSet_ = new float[numberSets_]; |
---|
| 204 | for (i=0;i<numberSets_;i++) { |
---|
| 205 | if (lower[i]>-1.0e20) |
---|
| 206 | lowerSet_[i]=lower[i]; |
---|
| 207 | else |
---|
| 208 | lowerSet_[i] = -1.0e30; |
---|
| 209 | } |
---|
| 210 | upperSet_ = new float[numberSets_]; |
---|
| 211 | for (i=0;i<numberSets_;i++) { |
---|
| 212 | if (upper[i]<1.0e20) |
---|
| 213 | upperSet_[i]=upper[i]; |
---|
| 214 | else |
---|
| 215 | upperSet_[i] = 1.0e30; |
---|
| 216 | } |
---|
| 217 | id_ = new int[numberGubInSmall]; |
---|
| 218 | for (i=0;i<numberGubInSmall;i++) |
---|
| 219 | id_[i]=-1; |
---|
| 220 | ClpPackedMatrix* originalMatrixA = |
---|
| 221 | dynamic_cast< ClpPackedMatrix*>(model->clpMatrix()); |
---|
| 222 | assert (originalMatrixA); |
---|
| 223 | CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix(); |
---|
[335] | 224 | originalMatrixA->setMatrixNull(); // so can be deleted safely |
---|
[328] | 225 | // guess how much space needed |
---|
| 226 | double guess = numberElements; |
---|
| 227 | guess /= (double) numberColumns; |
---|
[330] | 228 | guess *= 2*numberGubInSmall; |
---|
[328] | 229 | numberElements_ = (int) guess; |
---|
[399] | 230 | numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements(); |
---|
[328] | 231 | matrix_ = originalMatrix; |
---|
| 232 | zeroElements_ = false; |
---|
| 233 | // resize model (matrix stays same) |
---|
[399] | 234 | int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(frequency,numberRows))+1; |
---|
[328] | 235 | model->resize(newRowSize,numberNeeded); |
---|
| 236 | for (i=numberRows;i<newRowSize;i++) |
---|
| 237 | model->setRowStatus(i,ClpSimplex::basic); |
---|
| 238 | if (columnUpper_) { |
---|
| 239 | // set all upper bounds so we have enough space |
---|
| 240 | double * columnUpper = model->columnUpper(); |
---|
| 241 | for(i=firstDynamic_;i<lastDynamic_;i++) |
---|
| 242 | columnUpper[i]=1.0e10; |
---|
| 243 | } |
---|
| 244 | // resize matrix |
---|
| 245 | // extra 1 is so can keep number of elements handy |
---|
| 246 | originalMatrix->reserve(numberNeeded,numberElements_,true); |
---|
| 247 | originalMatrix->reserve(numberNeeded+1,numberElements_,false); |
---|
| 248 | originalMatrix->getMutableVectorStarts()[numberColumns]=originalMatrix->getNumElements(); |
---|
| 249 | originalMatrix->setDimensions(newRowSize,-1); |
---|
[344] | 250 | numberActiveColumns_=firstDynamic_; |
---|
[328] | 251 | // redo number of columns |
---|
| 252 | numberColumns = matrix_->getNumCols(); |
---|
| 253 | backToPivotRow_ = new int[numberNeeded]; |
---|
| 254 | keyVariable_ = new int[numberSets_]; |
---|
| 255 | if (status) { |
---|
| 256 | status_ = ClpCopyOfArray(status,numberSets_); |
---|
| 257 | assert (dynamicStatus); |
---|
| 258 | dynamicStatus_ = ClpCopyOfArray(dynamicStatus,numberGubColumns_); |
---|
| 259 | } else { |
---|
| 260 | status_= new unsigned char [numberSets_]; |
---|
| 261 | memset(status_,0,numberSets_); |
---|
| 262 | int i; |
---|
| 263 | for (i=0;i<numberSets_;i++) { |
---|
| 264 | // make slack key |
---|
| 265 | setStatus(i,ClpSimplex::basic); |
---|
| 266 | } |
---|
| 267 | dynamicStatus_ = new unsigned char [numberGubColumns_]; |
---|
| 268 | memset(dynamicStatus_,0,numberGubColumns_); // for clarity |
---|
| 269 | for (i=0;i<numberGubColumns_;i++) |
---|
| 270 | setDynamicStatus(i,atLowerBound); |
---|
| 271 | } |
---|
| 272 | toIndex_ = new int[numberSets_]; |
---|
| 273 | for (iSet=0;iSet<numberSets_;iSet++) |
---|
| 274 | toIndex_[iSet]=-1; |
---|
| 275 | fromIndex_ = new int [newRowSize-numberStaticRows_+1]; |
---|
| 276 | numberActiveSets_=0; |
---|
| 277 | rhsOffset_=NULL; |
---|
[330] | 278 | if (numberGubColumns_) { |
---|
| 279 | if (!status) { |
---|
| 280 | gubCrash(); |
---|
| 281 | } else { |
---|
| 282 | initialProblem(); |
---|
| 283 | } |
---|
[328] | 284 | } |
---|
| 285 | noCheck_ = -1; |
---|
| 286 | infeasibilityWeight_=0.0; |
---|
| 287 | } |
---|
| 288 | |
---|
| 289 | //------------------------------------------------------------------- |
---|
| 290 | // Destructor |
---|
| 291 | //------------------------------------------------------------------- |
---|
| 292 | ClpDynamicMatrix::~ClpDynamicMatrix () |
---|
| 293 | { |
---|
| 294 | delete [] backToPivotRow_; |
---|
| 295 | delete [] keyVariable_; |
---|
| 296 | delete [] toIndex_; |
---|
| 297 | delete [] fromIndex_; |
---|
| 298 | delete [] lowerSet_; |
---|
| 299 | delete [] upperSet_; |
---|
| 300 | delete [] status_; |
---|
| 301 | delete [] startSet_; |
---|
| 302 | delete [] next_; |
---|
| 303 | delete [] startColumn_; |
---|
| 304 | delete [] row_; |
---|
| 305 | delete [] element_; |
---|
| 306 | delete [] cost_; |
---|
| 307 | delete [] id_; |
---|
| 308 | delete [] dynamicStatus_; |
---|
| 309 | delete [] columnLower_; |
---|
| 310 | delete [] columnUpper_; |
---|
| 311 | } |
---|
| 312 | |
---|
| 313 | //---------------------------------------------------------------- |
---|
| 314 | // Assignment operator |
---|
| 315 | //------------------------------------------------------------------- |
---|
| 316 | ClpDynamicMatrix & |
---|
| 317 | ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs) |
---|
| 318 | { |
---|
| 319 | if (this != &rhs) { |
---|
| 320 | ClpPackedMatrix::operator=(rhs); |
---|
| 321 | delete [] backToPivotRow_; |
---|
| 322 | delete [] keyVariable_; |
---|
| 323 | delete [] toIndex_; |
---|
| 324 | delete [] fromIndex_; |
---|
| 325 | delete [] lowerSet_; |
---|
| 326 | delete [] upperSet_; |
---|
| 327 | delete [] status_; |
---|
| 328 | delete [] startSet_; |
---|
| 329 | delete [] next_; |
---|
| 330 | delete [] startColumn_; |
---|
| 331 | delete [] row_; |
---|
| 332 | delete [] element_; |
---|
| 333 | delete [] cost_; |
---|
| 334 | delete [] id_; |
---|
| 335 | delete [] dynamicStatus_; |
---|
| 336 | delete [] columnLower_; |
---|
| 337 | delete [] columnUpper_; |
---|
| 338 | objectiveOffset_ = rhs.objectiveOffset_; |
---|
| 339 | numberSets_ = rhs.numberSets_; |
---|
| 340 | numberActiveSets_ = rhs.numberActiveSets_; |
---|
| 341 | firstAvailable_ = rhs.firstAvailable_; |
---|
| 342 | firstAvailableBefore_ = rhs.firstAvailableBefore_; |
---|
| 343 | firstDynamic_ = rhs.firstDynamic_; |
---|
| 344 | lastDynamic_ = rhs.lastDynamic_; |
---|
| 345 | numberStaticRows_ = rhs.numberStaticRows_; |
---|
| 346 | numberElements_ = rhs.numberElements_; |
---|
| 347 | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_); |
---|
| 348 | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_); |
---|
| 349 | toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_); |
---|
| 350 | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_); |
---|
| 351 | lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_); |
---|
| 352 | upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_); |
---|
| 353 | status_ = ClpCopyOfArray(rhs.status_,numberSets_); |
---|
| 354 | model_ = rhs.model_; |
---|
| 355 | sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; |
---|
| 356 | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
---|
| 357 | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
---|
| 358 | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
---|
| 359 | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
---|
| 360 | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
---|
| 361 | savedBestGubDual_=rhs.savedBestGubDual_; |
---|
| 362 | savedBestSet_=rhs.savedBestSet_; |
---|
| 363 | noCheck_ = rhs.noCheck_; |
---|
| 364 | infeasibilityWeight_=rhs.infeasibilityWeight_; |
---|
| 365 | // Now secondary data |
---|
| 366 | numberGubColumns_ = rhs.numberGubColumns_; |
---|
| 367 | maximumGubColumns_=rhs.maximumGubColumns_; |
---|
| 368 | maximumElements_ = rhs.maximumElements_; |
---|
[330] | 369 | startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_); |
---|
[328] | 370 | next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_); |
---|
[330] | 371 | startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1); |
---|
| 372 | row_ = ClpCopyOfArray(rhs.row_,maximumElements_); |
---|
| 373 | element_ = ClpCopyOfArray(rhs.element_,maximumElements_); |
---|
| 374 | cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_); |
---|
[328] | 375 | id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_); |
---|
[330] | 376 | columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_); |
---|
| 377 | columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_); |
---|
| 378 | dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_); |
---|
[328] | 379 | } |
---|
| 380 | return *this; |
---|
| 381 | } |
---|
| 382 | //------------------------------------------------------------------- |
---|
| 383 | // Clone |
---|
| 384 | //------------------------------------------------------------------- |
---|
| 385 | ClpMatrixBase * ClpDynamicMatrix::clone() const |
---|
| 386 | { |
---|
| 387 | return new ClpDynamicMatrix(*this); |
---|
| 388 | } |
---|
| 389 | // Partial pricing |
---|
| 390 | void |
---|
[336] | 391 | ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, |
---|
[328] | 392 | int & bestSequence, int & numberWanted) |
---|
| 393 | { |
---|
[336] | 394 | numberWanted=currentWanted_; |
---|
[328] | 395 | assert(!model->rowScale()); |
---|
| 396 | if (numberSets_) { |
---|
| 397 | // Do packed part before gub |
---|
| 398 | // always??? |
---|
| 399 | //printf("normal packed price - start %d end %d (passed end %d, first %d)\n", |
---|
[336] | 400 | ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted); |
---|
[328] | 401 | } else { |
---|
| 402 | // no gub |
---|
[336] | 403 | ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted); |
---|
[328] | 404 | return; |
---|
| 405 | } |
---|
| 406 | if (numberWanted>0) { |
---|
| 407 | // and do some proportion of full set |
---|
[336] | 408 | int startG2 = (int) (startFraction*numberSets_); |
---|
| 409 | int endG2 = (int) (endFraction*numberSets_+0.1); |
---|
[399] | 410 | endG2 = CoinMin(endG2,numberSets_); |
---|
[328] | 411 | //printf("gub price - set start %d end %d\n", |
---|
| 412 | // startG2,endG2); |
---|
| 413 | double tolerance=model->currentDualTolerance(); |
---|
| 414 | double * reducedCost = model->djRegion(); |
---|
| 415 | const double * duals = model->dualRowSolution(); |
---|
| 416 | double bestDj; |
---|
| 417 | int numberRows = model->numberRows(); |
---|
| 418 | int slackOffset = lastDynamic_+numberRows; |
---|
| 419 | int structuralOffset = slackOffset+numberSets_; |
---|
| 420 | // If nothing found yet can go all the way to end |
---|
| 421 | int endAll = endG2; |
---|
| 422 | if (bestSequence<0&&!startG2) |
---|
| 423 | endAll = numberSets_; |
---|
| 424 | if (bestSequence>=0) { |
---|
| 425 | if (bestSequence!=savedBestSequence_) |
---|
| 426 | bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent |
---|
| 427 | else |
---|
| 428 | bestDj = savedBestDj_; |
---|
| 429 | } else { |
---|
| 430 | bestDj=tolerance; |
---|
| 431 | } |
---|
| 432 | int saveSequence = bestSequence; |
---|
| 433 | double djMod=0.0; |
---|
| 434 | double bestDjMod=0.0; |
---|
| 435 | //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(), |
---|
| 436 | // startG2,endG2,numberWanted); |
---|
| 437 | int bestSet=-1; |
---|
| 438 | #if 0 |
---|
| 439 | // make sure first available is clean (in case last iteration rejected) |
---|
| 440 | cost[firstAvailable_]=0.0; |
---|
| 441 | length[firstAvailable_]=0; |
---|
| 442 | model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0); |
---|
| 443 | model->setStatus(firstAvailable_,ClpSimplex::atLowerBound); |
---|
| 444 | { |
---|
| 445 | for (int i=firstAvailable_;i<lastDynamic_;i++) |
---|
| 446 | assert(!cost[i]); |
---|
| 447 | } |
---|
| 448 | #endif |
---|
[336] | 449 | int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_; |
---|
| 450 | int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_; |
---|
[328] | 451 | for (int iSet = startG2;iSet<endAll;iSet++) { |
---|
[336] | 452 | if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) { |
---|
[328] | 453 | // give up |
---|
| 454 | numberWanted=0; |
---|
| 455 | break; |
---|
[336] | 456 | } else if (iSet==endG2&&bestSequence>=0) { |
---|
[328] | 457 | break; |
---|
| 458 | } |
---|
| 459 | int gubRow = toIndex_[iSet]; |
---|
| 460 | if (gubRow>=0) { |
---|
| 461 | djMod = duals[gubRow+numberStaticRows_]; // have I got sign right? |
---|
| 462 | } else { |
---|
| 463 | int iBasic = keyVariable_[iSet]; |
---|
[330] | 464 | if (iBasic>=maximumGubColumns_) { |
---|
[328] | 465 | djMod = 0.0; // set not in |
---|
| 466 | } else { |
---|
| 467 | // get dj without |
---|
| 468 | djMod=0.0; |
---|
[329] | 469 | for (CoinBigIndex j=startColumn_[iBasic]; |
---|
[328] | 470 | j<startColumn_[iBasic+1];j++) { |
---|
| 471 | int jRow=row_[j]; |
---|
| 472 | djMod -= duals[jRow]*element_[j]; |
---|
| 473 | } |
---|
| 474 | djMod += cost_[iBasic]; |
---|
| 475 | // See if gub slack possible - dj is djMod |
---|
| 476 | if (getStatus(iSet)==ClpSimplex::atLowerBound) { |
---|
| 477 | double value = -djMod; |
---|
| 478 | if (value>tolerance) { |
---|
| 479 | numberWanted--; |
---|
| 480 | if (value>bestDj) { |
---|
| 481 | // check flagged variable and correct dj |
---|
| 482 | if (!flagged(iSet)) { |
---|
| 483 | bestDj=value; |
---|
| 484 | bestSequence = slackOffset+iSet; |
---|
| 485 | bestDjMod = djMod; |
---|
| 486 | bestSet = iSet; |
---|
| 487 | } else { |
---|
| 488 | // just to make sure we don't exit before got something |
---|
| 489 | numberWanted++; |
---|
| 490 | abort(); |
---|
| 491 | } |
---|
| 492 | } |
---|
| 493 | } |
---|
| 494 | } else if (getStatus(iSet)==ClpSimplex::atUpperBound) { |
---|
| 495 | double value = djMod; |
---|
| 496 | if (value>tolerance) { |
---|
| 497 | numberWanted--; |
---|
| 498 | if (value>bestDj) { |
---|
| 499 | // check flagged variable and correct dj |
---|
| 500 | if (!flagged(iSet)) { |
---|
| 501 | bestDj=value; |
---|
| 502 | bestSequence = slackOffset+iSet; |
---|
| 503 | bestDjMod = djMod; |
---|
| 504 | bestSet = iSet; |
---|
| 505 | } else { |
---|
| 506 | // just to make sure we don't exit before got something |
---|
| 507 | numberWanted++; |
---|
| 508 | abort(); |
---|
| 509 | } |
---|
| 510 | } |
---|
| 511 | } |
---|
| 512 | } |
---|
| 513 | } |
---|
| 514 | } |
---|
[329] | 515 | int iSequence= startSet_[iSet]; |
---|
| 516 | while (iSequence>=0) { |
---|
[328] | 517 | DynamicStatus status = getDynamicStatus(iSequence); |
---|
| 518 | if (status==atLowerBound||status==atUpperBound) { |
---|
| 519 | double value=cost_[iSequence]-djMod; |
---|
[329] | 520 | for (CoinBigIndex j=startColumn_[iSequence]; |
---|
[328] | 521 | j<startColumn_[iSequence+1];j++) { |
---|
| 522 | int jRow=row_[j]; |
---|
| 523 | value -= duals[jRow]*element_[j]; |
---|
| 524 | } |
---|
| 525 | // change sign if at lower bound |
---|
| 526 | if (status==atLowerBound) |
---|
| 527 | value = -value; |
---|
| 528 | if (value>tolerance) { |
---|
| 529 | numberWanted--; |
---|
| 530 | if (value>bestDj) { |
---|
| 531 | // check flagged variable and correct dj |
---|
| 532 | if (!flagged(iSequence)) { |
---|
| 533 | bestDj=value; |
---|
| 534 | bestSequence = structuralOffset+iSequence; |
---|
| 535 | bestDjMod = djMod; |
---|
| 536 | bestSet = iSet; |
---|
| 537 | } else { |
---|
| 538 | // just to make sure we don't exit before got something |
---|
| 539 | numberWanted++; |
---|
| 540 | } |
---|
| 541 | } |
---|
| 542 | } |
---|
| 543 | } |
---|
[329] | 544 | iSequence = next_[iSequence]; //onto next in set |
---|
[328] | 545 | } |
---|
| 546 | if (numberWanted<=0) { |
---|
| 547 | numberWanted=0; |
---|
| 548 | break; |
---|
| 549 | } |
---|
| 550 | } |
---|
| 551 | if (bestSequence!=saveSequence) { |
---|
| 552 | savedBestGubDual_ = bestDjMod; |
---|
| 553 | savedBestDj_ = bestDj; |
---|
| 554 | savedBestSequence_ = bestSequence; |
---|
| 555 | savedBestSet_ = bestSet; |
---|
| 556 | } |
---|
| 557 | // See if may be finished |
---|
| 558 | if (!startG2&&bestSequence<0) |
---|
| 559 | infeasibilityWeight_=model_->infeasibilityCost(); |
---|
[336] | 560 | else if (bestSequence>=0) |
---|
[328] | 561 | infeasibilityWeight_=-1.0; |
---|
| 562 | } |
---|
[336] | 563 | currentWanted_=numberWanted; |
---|
[328] | 564 | } |
---|
| 565 | /* Returns effective RHS if it is being used. This is used for long problems |
---|
| 566 | or big gub or anywhere where going through full columns is |
---|
| 567 | expensive. This may re-compute */ |
---|
| 568 | double * |
---|
| 569 | ClpDynamicMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh, |
---|
| 570 | bool check) |
---|
| 571 | { |
---|
[335] | 572 | //forceRefresh=true; |
---|
[328] | 573 | if (!model_->numberIterations()) |
---|
| 574 | forceRefresh=true; |
---|
| 575 | //check=false; |
---|
| 576 | #ifdef CLP_DEBUG |
---|
| 577 | double * saveE=NULL; |
---|
| 578 | if (rhsOffset_&&check) { |
---|
| 579 | int numberRows = model->numberRows(); |
---|
| 580 | saveE = new double[numberRows]; |
---|
| 581 | } |
---|
| 582 | #endif |
---|
| 583 | if (rhsOffset_) { |
---|
| 584 | #ifdef CLP_DEBUG |
---|
| 585 | if (check) { |
---|
| 586 | // no need - but check anyway |
---|
| 587 | int numberRows = model->numberRows(); |
---|
| 588 | double * rhs = new double[numberRows]; |
---|
| 589 | int iRow; |
---|
| 590 | int iSet; |
---|
| 591 | CoinZeroN(rhs,numberRows); |
---|
| 592 | // do ones at bounds before gub |
---|
| 593 | const double * smallSolution = model->solutionRegion(); |
---|
| 594 | const double * element =matrix_->getElements(); |
---|
| 595 | const int * row = matrix_->getIndices(); |
---|
| 596 | const CoinBigIndex * startColumn = matrix_->getVectorStarts(); |
---|
| 597 | const int * length = matrix_->getVectorLengths(); |
---|
| 598 | int iColumn; |
---|
| 599 | double objectiveOffset=0.0; |
---|
| 600 | for (iColumn=0;iColumn<firstDynamic_;iColumn++) { |
---|
| 601 | if (model->getStatus(iColumn)!=ClpSimplex::basic) { |
---|
| 602 | double value = smallSolution[iColumn]; |
---|
| 603 | for (CoinBigIndex j=startColumn[iColumn]; |
---|
| 604 | j<startColumn[iColumn]+length[iColumn];j++) { |
---|
| 605 | int jRow=row[j]; |
---|
| 606 | rhs[jRow] -= value*element[j]; |
---|
| 607 | } |
---|
| 608 | } |
---|
| 609 | } |
---|
| 610 | if (columnLower_||columnUpper_) { |
---|
| 611 | double * solution = new double [numberGubColumns_]; |
---|
| 612 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
[329] | 613 | int j= startSet_[iSet]; |
---|
| 614 | while (j>=0) { |
---|
[328] | 615 | double value=0.0; |
---|
| 616 | if (getDynamicStatus(j)!=inSmall) { |
---|
| 617 | if (getDynamicStatus(j)==atLowerBound) { |
---|
| 618 | if (columnLower_) |
---|
| 619 | value= columnLower_[j]; |
---|
| 620 | } else if (getDynamicStatus(j)==atUpperBound) { |
---|
| 621 | value= columnUpper_[j]; |
---|
| 622 | } else if (getDynamicStatus(j)==soloKey) { |
---|
| 623 | value =keyValue(iSet); |
---|
| 624 | } |
---|
| 625 | objectiveOffset += value*cost_[j]; |
---|
| 626 | } |
---|
| 627 | solution[j]=value; |
---|
[329] | 628 | j = next_[j]; //onto next in set |
---|
[328] | 629 | } |
---|
| 630 | } |
---|
| 631 | // ones in gub and in small problem |
---|
| 632 | for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) { |
---|
| 633 | if (model_->getStatus(iColumn)!=ClpSimplex::basic) { |
---|
| 634 | int jFull = id_[iColumn-firstDynamic_]; |
---|
| 635 | solution[jFull]=smallSolution[iColumn]; |
---|
| 636 | } |
---|
| 637 | } |
---|
| 638 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 639 | int kRow = toIndex_[iSet]; |
---|
| 640 | if (kRow>=0) |
---|
| 641 | kRow += numberStaticRows_; |
---|
[329] | 642 | int j= startSet_[iSet]; |
---|
| 643 | while (j>=0) { |
---|
[328] | 644 | double value = solution[j]; |
---|
| 645 | if (value) { |
---|
| 646 | for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) { |
---|
| 647 | int iRow = row_[k]; |
---|
| 648 | rhs[iRow] -= element_[k]*value; |
---|
| 649 | } |
---|
| 650 | if (kRow>=0) |
---|
| 651 | rhs[kRow] -= value; |
---|
| 652 | } |
---|
[329] | 653 | j = next_[j]; //onto next in set |
---|
[328] | 654 | } |
---|
| 655 | } |
---|
| 656 | delete [] solution; |
---|
| 657 | } else { |
---|
| 658 | // bounds |
---|
| 659 | ClpSimplex::Status iStatus; |
---|
| 660 | for (int iSet=0;iSet<numberSets_;iSet++) { |
---|
| 661 | int kRow = toIndex_[iSet]; |
---|
| 662 | if (kRow<0) { |
---|
| 663 | int iColumn = keyVariable_[iSet]; |
---|
[330] | 664 | if (iColumn<maximumGubColumns_) { |
---|
[328] | 665 | // key is not treated as basic |
---|
| 666 | double b=0.0; |
---|
| 667 | // key is structural - where is slack |
---|
| 668 | iStatus = getStatus(iSet); |
---|
| 669 | assert (iStatus!=ClpSimplex::basic); |
---|
| 670 | if (iStatus==ClpSimplex::atLowerBound) |
---|
| 671 | b=lowerSet_[iSet]; |
---|
| 672 | else |
---|
| 673 | b=upperSet_[iSet]; |
---|
| 674 | if (b) { |
---|
| 675 | objectiveOffset += b*cost_[iColumn]; |
---|
| 676 | for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) { |
---|
| 677 | int iRow = row_[j]; |
---|
| 678 | rhs[iRow] -= element_[j]*b; |
---|
| 679 | } |
---|
| 680 | } |
---|
| 681 | } |
---|
| 682 | } |
---|
| 683 | } |
---|
| 684 | } |
---|
| 685 | if (fabs(model->objectiveOffset()-objectiveOffset_-objectiveOffset)>1.0e-1) |
---|
| 686 | printf("old offset %g, true %g\n",model->objectiveOffset(), |
---|
| 687 | objectiveOffset_-objectiveOffset); |
---|
| 688 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
| 689 | if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3) |
---|
| 690 | printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]); |
---|
| 691 | } |
---|
| 692 | memcpy(saveE,rhs,numberRows*sizeof(double)); |
---|
| 693 | delete [] rhs; |
---|
| 694 | } |
---|
| 695 | #endif |
---|
| 696 | if (forceRefresh||(refreshFrequency_&&model->numberIterations()>= |
---|
| 697 | lastRefresh_+refreshFrequency_)) { |
---|
| 698 | int numberRows = model->numberRows(); |
---|
| 699 | int iSet; |
---|
| 700 | CoinZeroN(rhsOffset_,numberRows); |
---|
| 701 | // do ones at bounds before gub |
---|
| 702 | const double * smallSolution = model->solutionRegion(); |
---|
| 703 | const double * element =matrix_->getElements(); |
---|
| 704 | const int * row = matrix_->getIndices(); |
---|
| 705 | const CoinBigIndex * startColumn = matrix_->getVectorStarts(); |
---|
| 706 | const int * length = matrix_->getVectorLengths(); |
---|
| 707 | int iColumn; |
---|
| 708 | double objectiveOffset=0.0; |
---|
| 709 | for (iColumn=0;iColumn<firstDynamic_;iColumn++) { |
---|
| 710 | if (model->getStatus(iColumn)!=ClpSimplex::basic) { |
---|
| 711 | double value = smallSolution[iColumn]; |
---|
| 712 | for (CoinBigIndex j=startColumn[iColumn]; |
---|
| 713 | j<startColumn[iColumn]+length[iColumn];j++) { |
---|
| 714 | int jRow=row[j]; |
---|
| 715 | rhsOffset_[jRow] -= value*element[j]; |
---|
| 716 | } |
---|
| 717 | } |
---|
| 718 | } |
---|
| 719 | if (columnLower_||columnUpper_) { |
---|
| 720 | double * solution = new double [numberGubColumns_]; |
---|
| 721 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
[329] | 722 | int j= startSet_[iSet]; |
---|
| 723 | while (j>=0) { |
---|
[328] | 724 | double value=0.0; |
---|
| 725 | if (getDynamicStatus(j)!=inSmall) { |
---|
| 726 | if (getDynamicStatus(j)==atLowerBound) { |
---|
| 727 | if (columnLower_) |
---|
| 728 | value= columnLower_[j]; |
---|
| 729 | } else if (getDynamicStatus(j)==atUpperBound) { |
---|
| 730 | value= columnUpper_[j]; |
---|
| 731 | } else if (getDynamicStatus(j)==soloKey) { |
---|
| 732 | value =keyValue(iSet); |
---|
| 733 | } |
---|
| 734 | objectiveOffset += value*cost_[j]; |
---|
| 735 | } |
---|
| 736 | solution[j]=value; |
---|
[329] | 737 | j = next_[j]; //onto next in set |
---|
[328] | 738 | } |
---|
| 739 | } |
---|
| 740 | // ones in gub and in small problem |
---|
| 741 | for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) { |
---|
| 742 | if (model_->getStatus(iColumn)!=ClpSimplex::basic) { |
---|
| 743 | int jFull = id_[iColumn-firstDynamic_]; |
---|
| 744 | solution[jFull]=smallSolution[iColumn]; |
---|
| 745 | } |
---|
| 746 | } |
---|
| 747 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 748 | int kRow = toIndex_[iSet]; |
---|
| 749 | if (kRow>=0) |
---|
| 750 | kRow += numberStaticRows_; |
---|
[329] | 751 | int j= startSet_[iSet]; |
---|
| 752 | while (j>=0) { |
---|
[328] | 753 | double value = solution[j]; |
---|
| 754 | if (value) { |
---|
| 755 | for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) { |
---|
| 756 | int iRow = row_[k]; |
---|
| 757 | rhsOffset_[iRow] -= element_[k]*value; |
---|
| 758 | } |
---|
| 759 | if (kRow>=0) |
---|
| 760 | rhsOffset_[kRow] -= value; |
---|
| 761 | } |
---|
[329] | 762 | j = next_[j]; //onto next in set |
---|
[328] | 763 | } |
---|
| 764 | } |
---|
| 765 | delete [] solution; |
---|
| 766 | } else { |
---|
| 767 | // bounds |
---|
| 768 | ClpSimplex::Status iStatus; |
---|
| 769 | for (int iSet=0;iSet<numberSets_;iSet++) { |
---|
| 770 | int kRow = toIndex_[iSet]; |
---|
| 771 | if (kRow<0) { |
---|
| 772 | int iColumn = keyVariable_[iSet]; |
---|
[330] | 773 | if (iColumn<maximumGubColumns_) { |
---|
[328] | 774 | // key is not treated as basic |
---|
| 775 | double b=0.0; |
---|
| 776 | // key is structural - where is slack |
---|
| 777 | iStatus = getStatus(iSet); |
---|
| 778 | assert (iStatus!=ClpSimplex::basic); |
---|
| 779 | if (iStatus==ClpSimplex::atLowerBound) |
---|
| 780 | b=lowerSet_[iSet]; |
---|
| 781 | else |
---|
| 782 | b=upperSet_[iSet]; |
---|
| 783 | if (b) { |
---|
| 784 | objectiveOffset += b*cost_[iColumn]; |
---|
| 785 | for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) { |
---|
| 786 | int iRow = row_[j]; |
---|
| 787 | rhsOffset_[iRow] -= element_[j]*b; |
---|
| 788 | } |
---|
| 789 | } |
---|
| 790 | } |
---|
| 791 | } |
---|
| 792 | } |
---|
| 793 | } |
---|
| 794 | model->setObjectiveOffset(objectiveOffset_-objectiveOffset); |
---|
| 795 | #ifdef CLP_DEBUG |
---|
| 796 | if (saveE) { |
---|
| 797 | int iRow; |
---|
| 798 | for (iRow=0;iRow<numberRows;iRow++) { |
---|
| 799 | if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3) |
---|
| 800 | printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]); |
---|
| 801 | } |
---|
| 802 | delete [] saveE; |
---|
| 803 | } |
---|
| 804 | #endif |
---|
| 805 | lastRefresh_ = model->numberIterations(); |
---|
| 806 | } |
---|
| 807 | } |
---|
| 808 | return rhsOffset_; |
---|
| 809 | } |
---|
| 810 | /* |
---|
| 811 | update information for a pivot (and effective rhs) |
---|
| 812 | */ |
---|
| 813 | int |
---|
| 814 | ClpDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue) |
---|
| 815 | { |
---|
| 816 | |
---|
| 817 | // now update working model |
---|
| 818 | int sequenceIn = model->sequenceIn(); |
---|
| 819 | int sequenceOut = model->sequenceOut(); |
---|
| 820 | int numberColumns = model->numberColumns(); |
---|
| 821 | if (sequenceIn!=sequenceOut&&sequenceIn<numberColumns) |
---|
| 822 | backToPivotRow_[sequenceIn]=model->pivotRow(); |
---|
| 823 | if (sequenceIn>=firstDynamic_&&sequenceIn<numberColumns) { |
---|
| 824 | int bigSequence = id_[sequenceIn-firstDynamic_]; |
---|
| 825 | if (getDynamicStatus(bigSequence)!=inSmall) { |
---|
| 826 | firstAvailable_++; |
---|
| 827 | setDynamicStatus(bigSequence,inSmall); |
---|
| 828 | } |
---|
| 829 | } |
---|
| 830 | // make sure slack is synchronized |
---|
| 831 | if (sequenceIn>=numberColumns+numberStaticRows_) { |
---|
| 832 | int iDynamic = sequenceIn-numberColumns-numberStaticRows_; |
---|
| 833 | int iSet = fromIndex_[iDynamic]; |
---|
| 834 | setStatus(iSet,model->getStatus(sequenceIn)); |
---|
| 835 | } |
---|
| 836 | if (sequenceOut>=numberColumns+numberStaticRows_) { |
---|
| 837 | int iDynamic = sequenceOut-numberColumns-numberStaticRows_; |
---|
| 838 | int iSet = fromIndex_[iDynamic]; |
---|
[341] | 839 | // out may have gone through barrier - so check |
---|
| 840 | double valueOut = model->lowerRegion()[sequenceOut]; |
---|
| 841 | if (fabs(valueOut -(double) lowerSet_[iSet])< |
---|
| 842 | fabs(valueOut -(double) upperSet_[iSet])) |
---|
| 843 | setStatus(iSet,ClpSimplex::atLowerBound); |
---|
| 844 | else |
---|
| 845 | setStatus(iSet,ClpSimplex::atUpperBound); |
---|
[344] | 846 | if (lowerSet_[iSet]==upperSet_[iSet]) |
---|
| 847 | setStatus(iSet,ClpSimplex::isFixed); |
---|
[341] | 848 | if (getStatus(iSet)!=model->getStatus(sequenceOut)) |
---|
| 849 | printf("** set %d status %d, var status %d\n",iSet, |
---|
| 850 | getStatus(iSet),model->getStatus(sequenceOut)); |
---|
[328] | 851 | } |
---|
| 852 | ClpMatrixBase::updatePivot(model,oldInValue,oldOutValue); |
---|
| 853 | #ifdef CLP_DEBUG |
---|
| 854 | char * inSmall = new char [numberGubColumns_]; |
---|
| 855 | memset(inSmall,0,numberGubColumns_); |
---|
[341] | 856 | const double * solution = model->solutionRegion(); |
---|
[328] | 857 | for (int i=0;i<numberGubColumns_;i++) |
---|
| 858 | if (getDynamicStatus(i)==ClpDynamicMatrix::inSmall) |
---|
| 859 | inSmall[i]=1; |
---|
| 860 | for (int i=firstDynamic_;i<firstAvailable_;i++) { |
---|
| 861 | int k=id_[i-firstDynamic_]; |
---|
| 862 | inSmall[k]=0; |
---|
[341] | 863 | //if (k>=23289&&k<23357&&solution[i]) |
---|
| 864 | //printf("var %d (in small %d) has value %g\n",k,i,solution[i]); |
---|
[328] | 865 | } |
---|
| 866 | for (int i=0;i<numberGubColumns_;i++) |
---|
| 867 | assert (!inSmall[i]); |
---|
| 868 | delete [] inSmall; |
---|
[332] | 869 | #ifndef NDEBUG |
---|
[328] | 870 | for (int i=0;i<numberActiveSets_;i++) { |
---|
| 871 | int iSet = fromIndex_[i]; |
---|
| 872 | assert (toIndex_[iSet]==i); |
---|
[341] | 873 | //if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_)) |
---|
| 874 | //printf("*** set %d status %d, var status %d\n",iSet, |
---|
| 875 | // getStatus(iSet),model->getRowStatus(i+numberStaticRows_)); |
---|
| 876 | //assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet)); |
---|
| 877 | //if (iSet==1035) { |
---|
| 878 | //printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_], |
---|
| 879 | // model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]); |
---|
[344] | 880 | //} |
---|
[328] | 881 | } |
---|
| 882 | #endif |
---|
[332] | 883 | #endif |
---|
[328] | 884 | return 0; |
---|
| 885 | } |
---|
| 886 | /* |
---|
| 887 | utility dual function for dealing with dynamic constraints |
---|
| 888 | mode=n see ClpGubMatrix.hpp for definition |
---|
| 889 | Remember to update here when settled down |
---|
| 890 | */ |
---|
| 891 | void |
---|
| 892 | ClpDynamicMatrix::dualExpanded(ClpSimplex * model, |
---|
| 893 | CoinIndexedVector * array, |
---|
| 894 | double * other,int mode) |
---|
| 895 | { |
---|
| 896 | switch (mode) { |
---|
| 897 | // modify costs before transposeUpdate |
---|
| 898 | case 0: |
---|
| 899 | break; |
---|
| 900 | // create duals for key variables (without check on dual infeasible) |
---|
| 901 | case 1: |
---|
| 902 | break; |
---|
| 903 | // as 1 but check slacks and compute djs |
---|
| 904 | case 2: |
---|
| 905 | { |
---|
| 906 | // do pivots |
---|
| 907 | int * pivotVariable = model->pivotVariable(); |
---|
| 908 | int numberRows=numberStaticRows_+numberActiveSets_; |
---|
| 909 | int numberColumns = model->numberColumns(); |
---|
| 910 | for (int iRow=0;iRow<numberRows;iRow++) { |
---|
| 911 | int iPivot=pivotVariable[iRow]; |
---|
| 912 | if (iPivot<numberColumns) |
---|
| 913 | backToPivotRow_[iPivot]=iRow; |
---|
| 914 | } |
---|
| 915 | if (noCheck_>=0) { |
---|
| 916 | if (infeasibilityWeight_!=model_->infeasibilityCost()) { |
---|
| 917 | // don't bother checking |
---|
| 918 | sumDualInfeasibilities_=100.0; |
---|
| 919 | numberDualInfeasibilities_=1; |
---|
| 920 | sumOfRelaxedDualInfeasibilities_ =100.0; |
---|
| 921 | return; |
---|
| 922 | } |
---|
| 923 | } |
---|
| 924 | // In theory we should subtract out ones we have done but .... |
---|
| 925 | // If key slack then dual 0.0 |
---|
| 926 | // If not then slack could be dual infeasible |
---|
| 927 | // dj for key is zero so that defines dual on set |
---|
| 928 | int i; |
---|
| 929 | double * dual = model->dualRowSolution(); |
---|
| 930 | double dualTolerance = model->dualTolerance(); |
---|
| 931 | double relaxedTolerance=dualTolerance; |
---|
| 932 | // we can't really trust infeasibilities if there is dual error |
---|
[800] | 933 | double error = CoinMin(1.0e-2,model->largestDualError()); |
---|
[328] | 934 | // allow tolerance at least slightly bigger than standard |
---|
| 935 | relaxedTolerance = relaxedTolerance + error; |
---|
| 936 | // but we will be using difference |
---|
| 937 | relaxedTolerance -= dualTolerance; |
---|
| 938 | sumDualInfeasibilities_=0.0; |
---|
| 939 | numberDualInfeasibilities_=0; |
---|
| 940 | sumOfRelaxedDualInfeasibilities_ =0.0; |
---|
| 941 | for (i=0;i<numberSets_;i++) { |
---|
| 942 | double value=0.0; |
---|
| 943 | int gubRow = toIndex_[i]; |
---|
| 944 | if (gubRow<0) { |
---|
| 945 | int kColumn = keyVariable_[i]; |
---|
[330] | 946 | if (kColumn<maximumGubColumns_) { |
---|
[328] | 947 | // dj without set |
---|
| 948 | value = cost_[kColumn]; |
---|
| 949 | for (CoinBigIndex j=startColumn_[kColumn]; |
---|
| 950 | j<startColumn_[kColumn+1];j++) { |
---|
| 951 | int iRow = row_[j]; |
---|
| 952 | value -= dual[iRow]*element_[j]; |
---|
| 953 | } |
---|
| 954 | double infeasibility = 0.0; |
---|
| 955 | if (getStatus(i)==ClpSimplex::atLowerBound) { |
---|
| 956 | if (-value>dualTolerance) |
---|
| 957 | infeasibility=-value-dualTolerance; |
---|
| 958 | } else if (getStatus(i)==ClpSimplex::atUpperBound) { |
---|
| 959 | if (value>dualTolerance) |
---|
| 960 | infeasibility=value-dualTolerance; |
---|
| 961 | } |
---|
| 962 | if (infeasibility>0.0) { |
---|
| 963 | sumDualInfeasibilities_ += infeasibility; |
---|
| 964 | if (infeasibility>relaxedTolerance) |
---|
| 965 | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
---|
| 966 | numberDualInfeasibilities_ ++; |
---|
| 967 | } |
---|
| 968 | } |
---|
| 969 | } else { |
---|
| 970 | value = dual[gubRow+numberStaticRows_]; |
---|
| 971 | } |
---|
| 972 | // Now subtract out from all |
---|
[329] | 973 | int k= startSet_[i]; |
---|
| 974 | while (k>=0) { |
---|
[328] | 975 | if (getDynamicStatus(k)!=inSmall) { |
---|
| 976 | double djValue = cost_[k]-value; |
---|
| 977 | for (CoinBigIndex j=startColumn_[k]; |
---|
| 978 | j<startColumn_[k+1];j++) { |
---|
| 979 | int iRow = row_[j]; |
---|
| 980 | djValue -= dual[iRow]*element_[j]; |
---|
| 981 | } |
---|
| 982 | double infeasibility=0.0; |
---|
| 983 | if (getDynamicStatus(k)==atLowerBound) { |
---|
| 984 | if (djValue<-dualTolerance) |
---|
| 985 | infeasibility=-djValue-dualTolerance; |
---|
| 986 | } else if (getDynamicStatus(k)==atUpperBound) { |
---|
| 987 | // at upper bound |
---|
| 988 | if (djValue>dualTolerance) |
---|
| 989 | infeasibility=djValue-dualTolerance; |
---|
| 990 | } |
---|
| 991 | if (infeasibility>0.0) { |
---|
| 992 | sumDualInfeasibilities_ += infeasibility; |
---|
| 993 | if (infeasibility>relaxedTolerance) |
---|
| 994 | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
---|
| 995 | numberDualInfeasibilities_ ++; |
---|
| 996 | } |
---|
| 997 | } |
---|
[329] | 998 | k = next_[k]; //onto next in set |
---|
[328] | 999 | } |
---|
| 1000 | } |
---|
| 1001 | } |
---|
| 1002 | infeasibilityWeight_=-1.0; |
---|
| 1003 | break; |
---|
| 1004 | // Report on infeasibilities of key variables |
---|
| 1005 | case 3: |
---|
| 1006 | { |
---|
| 1007 | model->setSumDualInfeasibilities(model->sumDualInfeasibilities()+ |
---|
| 1008 | sumDualInfeasibilities_); |
---|
| 1009 | model->setNumberDualInfeasibilities(model->numberDualInfeasibilities()+ |
---|
| 1010 | numberDualInfeasibilities_); |
---|
| 1011 | model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities()+ |
---|
| 1012 | sumOfRelaxedDualInfeasibilities_); |
---|
| 1013 | } |
---|
| 1014 | break; |
---|
| 1015 | // modify costs before transposeUpdate for partial pricing |
---|
| 1016 | case 4: |
---|
| 1017 | break; |
---|
| 1018 | } |
---|
| 1019 | } |
---|
| 1020 | /* |
---|
| 1021 | general utility function for dealing with dynamic constraints |
---|
| 1022 | mode=n see ClpGubMatrix.hpp for definition |
---|
| 1023 | Remember to update here when settled down |
---|
| 1024 | */ |
---|
| 1025 | int |
---|
| 1026 | ClpDynamicMatrix::generalExpanded(ClpSimplex * model,int mode,int &number) |
---|
| 1027 | { |
---|
| 1028 | int returnCode=0; |
---|
| 1029 | switch (mode) { |
---|
| 1030 | // Fill in pivotVariable |
---|
| 1031 | case 0: |
---|
| 1032 | { |
---|
| 1033 | // If no effective rhs - form it |
---|
| 1034 | if (!rhsOffset_) { |
---|
| 1035 | rhsOffset_ = new double[model->numberRows()]; |
---|
| 1036 | rhsOffset(model,true); |
---|
| 1037 | } |
---|
| 1038 | int i; |
---|
| 1039 | int numberBasic=number; |
---|
| 1040 | int numberColumns = model->numberColumns(); |
---|
| 1041 | // Use different array so can build from true pivotVariable_ |
---|
| 1042 | //int * pivotVariable = model->pivotVariable(); |
---|
| 1043 | int * pivotVariable = model->rowArray(0)->getIndices(); |
---|
| 1044 | for (i=0;i<numberColumns;i++) { |
---|
| 1045 | if (model->getColumnStatus(i) == ClpSimplex::basic) |
---|
| 1046 | pivotVariable[numberBasic++]=i; |
---|
| 1047 | } |
---|
| 1048 | number = numberBasic; |
---|
| 1049 | } |
---|
| 1050 | break; |
---|
| 1051 | // Do initial extra rows + maximum basic |
---|
| 1052 | case 2: |
---|
| 1053 | { |
---|
| 1054 | number = model->numberRows(); |
---|
| 1055 | } |
---|
| 1056 | break; |
---|
[335] | 1057 | // Before normal replaceColumn |
---|
| 1058 | case 3: |
---|
| 1059 | { |
---|
| 1060 | if (numberActiveSets_+numberStaticRows_==model_->numberRows()) { |
---|
| 1061 | // no space - re-factorize |
---|
| 1062 | returnCode=4; |
---|
| 1063 | number=-1; // say no need for normal replaceColumn |
---|
| 1064 | } |
---|
| 1065 | } |
---|
| 1066 | break; |
---|
[328] | 1067 | // To see if can dual or primal |
---|
| 1068 | case 4: |
---|
| 1069 | { |
---|
| 1070 | returnCode= 1; |
---|
| 1071 | } |
---|
| 1072 | break; |
---|
| 1073 | // save status |
---|
| 1074 | case 5: |
---|
| 1075 | { |
---|
| 1076 | } |
---|
| 1077 | break; |
---|
| 1078 | // restore status |
---|
| 1079 | case 6: |
---|
| 1080 | { |
---|
| 1081 | // maybe mismatch - redo all in small model |
---|
| 1082 | //for (int i=firstDynamic_;i<firstAvailable_;i++) { |
---|
| 1083 | //int jColumn = id_[i-firstDynamic_]; |
---|
| 1084 | //setDynamicStatus(jColumn,inSmall); |
---|
| 1085 | //} |
---|
| 1086 | } |
---|
| 1087 | break; |
---|
| 1088 | // unflag all variables |
---|
| 1089 | case 8: |
---|
| 1090 | { |
---|
| 1091 | for (int i=0;i<numberGubColumns_;i++) { |
---|
| 1092 | if (flagged(i)) { |
---|
| 1093 | unsetFlagged(i); |
---|
| 1094 | returnCode++; |
---|
| 1095 | } |
---|
| 1096 | } |
---|
| 1097 | } |
---|
| 1098 | break; |
---|
| 1099 | // redo costs in primal |
---|
| 1100 | case 9: |
---|
| 1101 | { |
---|
| 1102 | double * cost = model->costRegion(); |
---|
| 1103 | double * solution = model->solutionRegion(); |
---|
| 1104 | double * columnLower = model->lowerRegion(); |
---|
| 1105 | double * columnUpper = model->upperRegion(); |
---|
| 1106 | int i; |
---|
| 1107 | bool doCosts = (number&4)!=0; |
---|
| 1108 | bool doBounds = (number&1)!=0; |
---|
| 1109 | for ( i=firstDynamic_;i<firstAvailable_;i++) { |
---|
| 1110 | int jColumn = id_[i-firstDynamic_]; |
---|
| 1111 | if (doBounds) { |
---|
| 1112 | if (!columnLower_&&!columnUpper_) { |
---|
| 1113 | columnLower[i] = 0.0; |
---|
| 1114 | columnUpper[i] = COIN_DBL_MAX; |
---|
| 1115 | } else { |
---|
| 1116 | if (columnLower_) |
---|
| 1117 | columnLower[i] = columnLower_[jColumn]; |
---|
| 1118 | else |
---|
| 1119 | columnLower[i] = 0.0; |
---|
| 1120 | if (columnUpper_) |
---|
| 1121 | columnUpper[i] = columnUpper_[jColumn]; |
---|
| 1122 | else |
---|
| 1123 | columnUpper[i] = COIN_DBL_MAX; |
---|
| 1124 | } |
---|
| 1125 | } |
---|
| 1126 | if (doCosts) { |
---|
| 1127 | cost[i]=cost_[jColumn]; |
---|
[332] | 1128 | // Original bounds |
---|
[328] | 1129 | if (model->nonLinearCost()) |
---|
| 1130 | model->nonLinearCost()->setOne(i,solution[i], |
---|
[332] | 1131 | this->columnLower(jColumn), |
---|
| 1132 | this->columnUpper(jColumn),cost_[jColumn]); |
---|
[328] | 1133 | } |
---|
| 1134 | } |
---|
| 1135 | // and active sets |
---|
| 1136 | for (i=0;i<numberActiveSets_;i++ ) { |
---|
| 1137 | int iSet = fromIndex_[i]; |
---|
| 1138 | int iSequence=lastDynamic_+numberStaticRows_+i; |
---|
| 1139 | if (doBounds) { |
---|
| 1140 | if (lowerSet_[iSet]>-1.0e20) |
---|
| 1141 | columnLower[iSequence] = lowerSet_[iSet]; |
---|
| 1142 | else |
---|
| 1143 | columnLower[iSequence]=-COIN_DBL_MAX; |
---|
| 1144 | if (upperSet_[iSet]<1.0e20) |
---|
| 1145 | columnUpper[iSequence] = upperSet_[iSet]; |
---|
| 1146 | else |
---|
| 1147 | columnUpper[iSequence]=COIN_DBL_MAX; |
---|
| 1148 | } |
---|
| 1149 | if (doCosts) { |
---|
[341] | 1150 | if (model->nonLinearCost()) { |
---|
| 1151 | double trueLower; |
---|
| 1152 | if (lowerSet_[iSet]>-1.0e20) |
---|
| 1153 | trueLower = lowerSet_[iSet]; |
---|
| 1154 | else |
---|
| 1155 | trueLower=-COIN_DBL_MAX; |
---|
| 1156 | double trueUpper; |
---|
| 1157 | if (upperSet_[iSet]<1.0e20) |
---|
| 1158 | trueUpper = upperSet_[iSet]; |
---|
| 1159 | else |
---|
| 1160 | trueUpper=COIN_DBL_MAX; |
---|
[328] | 1161 | model->nonLinearCost()->setOne(iSequence,solution[iSequence], |
---|
[341] | 1162 | trueLower,trueUpper,0.0); |
---|
| 1163 | } |
---|
[328] | 1164 | } |
---|
| 1165 | } |
---|
| 1166 | } |
---|
| 1167 | break; |
---|
| 1168 | // return 1 if there may be changing bounds on variable (column generation) |
---|
| 1169 | case 10: |
---|
| 1170 | { |
---|
| 1171 | // return 1 as bounds on rhs will change |
---|
| 1172 | returnCode = 1; |
---|
| 1173 | } |
---|
| 1174 | break; |
---|
| 1175 | // make sure set is clean |
---|
| 1176 | case 7: |
---|
| 1177 | { |
---|
| 1178 | // first flag |
---|
| 1179 | if (number>=firstDynamic_&&number<lastDynamic_) { |
---|
| 1180 | assert (number==model->sequenceIn()); |
---|
| 1181 | // id will be sitting at firstAvailable |
---|
| 1182 | int sequence = id_[firstAvailable_-firstDynamic_]; |
---|
| 1183 | setFlagged(sequence); |
---|
| 1184 | model->clearFlagged(firstAvailable_); |
---|
| 1185 | } |
---|
| 1186 | } |
---|
| 1187 | case 11: |
---|
| 1188 | { |
---|
| 1189 | int sequenceIn = model->sequenceIn(); |
---|
| 1190 | if (sequenceIn>=firstDynamic_&&sequenceIn<lastDynamic_) { |
---|
[335] | 1191 | assert (number==model->sequenceIn()); |
---|
| 1192 | // take out variable (but leave key) |
---|
[328] | 1193 | double * cost = model->costRegion(); |
---|
| 1194 | double * columnLower = model->lowerRegion(); |
---|
| 1195 | double * columnUpper = model->upperRegion(); |
---|
| 1196 | double * solution = model->solutionRegion(); |
---|
| 1197 | int * length = matrix_->getMutableVectorLengths(); |
---|
| 1198 | #ifndef NDEBUG |
---|
| 1199 | if (length[sequenceIn]) { |
---|
| 1200 | int * row = matrix_->getMutableIndices(); |
---|
| 1201 | CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1202 | int iRow = row[startColumn[sequenceIn]+length[sequenceIn]-1]; |
---|
| 1203 | int which = iRow-numberStaticRows_; |
---|
| 1204 | assert (which>=0); |
---|
| 1205 | int iSet = fromIndex_[which]; |
---|
| 1206 | assert (toIndex_[iSet]==which); |
---|
| 1207 | } |
---|
| 1208 | #endif |
---|
| 1209 | // no need firstAvailable_--; |
---|
| 1210 | solution[firstAvailable_]=0.0; |
---|
| 1211 | cost[firstAvailable_]=0.0; |
---|
| 1212 | length[firstAvailable_]=0; |
---|
| 1213 | model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0); |
---|
| 1214 | model->setStatus(firstAvailable_,ClpSimplex::atLowerBound); |
---|
| 1215 | columnLower[firstAvailable_]=0.0; |
---|
| 1216 | columnUpper[firstAvailable_]=COIN_DBL_MAX; |
---|
| 1217 | |
---|
| 1218 | // not really in small problem |
---|
| 1219 | int iBig=id_[sequenceIn-firstDynamic_]; |
---|
| 1220 | if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound) { |
---|
| 1221 | setDynamicStatus(iBig,atLowerBound); |
---|
| 1222 | if (columnLower_) |
---|
| 1223 | modifyOffset(sequenceIn,columnLower_[iBig]); |
---|
| 1224 | } else { |
---|
| 1225 | setDynamicStatus(iBig,atUpperBound); |
---|
| 1226 | modifyOffset(sequenceIn,columnUpper_[iBig]); |
---|
| 1227 | } |
---|
| 1228 | } |
---|
| 1229 | } |
---|
| 1230 | break; |
---|
| 1231 | default: |
---|
| 1232 | break; |
---|
| 1233 | } |
---|
| 1234 | return returnCode; |
---|
| 1235 | } |
---|
| 1236 | /* Purely for column generation and similar ideas. Allows |
---|
| 1237 | matrix and any bounds or costs to be updated (sensibly). |
---|
| 1238 | Returns non-zero if any changes. |
---|
| 1239 | */ |
---|
| 1240 | int |
---|
| 1241 | ClpDynamicMatrix::refresh(ClpSimplex * model) |
---|
| 1242 | { |
---|
| 1243 | // If at beginning then create initial problem |
---|
| 1244 | if (firstDynamic_==firstAvailable_) { |
---|
| 1245 | initialProblem(); |
---|
| 1246 | return 1; |
---|
| 1247 | } else if (!model->nonLinearCost()) { |
---|
| 1248 | // will be same as last time |
---|
| 1249 | return 1; |
---|
| 1250 | } |
---|
| 1251 | // lookup array |
---|
| 1252 | int * active = new int [numberActiveSets_]; |
---|
| 1253 | CoinZeroN(active,numberActiveSets_); |
---|
| 1254 | int iColumn; |
---|
| 1255 | int numberColumns = model->numberColumns(); |
---|
| 1256 | double * element = matrix_->getMutableElements(); |
---|
| 1257 | int * row = matrix_->getMutableIndices(); |
---|
| 1258 | CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1259 | int * length = matrix_->getMutableVectorLengths(); |
---|
| 1260 | double * cost = model->costRegion(); |
---|
| 1261 | double * columnLower = model->lowerRegion(); |
---|
| 1262 | double * columnUpper = model->upperRegion(); |
---|
| 1263 | CoinBigIndex numberElements=startColumn[firstDynamic_]; |
---|
| 1264 | // first just do lookup and basic stuff |
---|
| 1265 | int currentNumber = firstAvailable_; |
---|
| 1266 | firstAvailable_ = firstDynamic_; |
---|
| 1267 | double objectiveChange=0.0; |
---|
| 1268 | double * solution = model->solutionRegion(); |
---|
| 1269 | int currentNumberActiveSets=numberActiveSets_; |
---|
| 1270 | for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) { |
---|
| 1271 | int iRow = row[startColumn[iColumn]+length[iColumn]-1]; |
---|
| 1272 | int which = iRow-numberStaticRows_; |
---|
| 1273 | assert (which>=0); |
---|
| 1274 | int iSet = fromIndex_[which]; |
---|
| 1275 | assert (toIndex_[iSet]==which); |
---|
| 1276 | if (model->getStatus(iColumn)==ClpSimplex::basic) { |
---|
| 1277 | active[which]++; |
---|
| 1278 | // may as well make key |
---|
| 1279 | keyVariable_[iSet]=id_[iColumn-firstDynamic_]; |
---|
| 1280 | } |
---|
| 1281 | } |
---|
| 1282 | int i; |
---|
| 1283 | numberActiveSets_=0; |
---|
[330] | 1284 | int numberDeleted=0; |
---|
[328] | 1285 | for (i=0;i<currentNumberActiveSets;i++) { |
---|
| 1286 | int iRow = i+numberStaticRows_; |
---|
| 1287 | int numberActive = active[i]; |
---|
| 1288 | int iSet = fromIndex_[i]; |
---|
| 1289 | if (model->getRowStatus(iRow)==ClpSimplex::basic) { |
---|
| 1290 | numberActive++; |
---|
| 1291 | // may as well make key |
---|
[330] | 1292 | keyVariable_[iSet]=maximumGubColumns_+iSet; |
---|
[328] | 1293 | } |
---|
| 1294 | if (numberActive>1) { |
---|
| 1295 | // keep |
---|
| 1296 | active[i]=numberActiveSets_; |
---|
| 1297 | numberActiveSets_++; |
---|
| 1298 | } else { |
---|
| 1299 | active[i]=-1; |
---|
| 1300 | } |
---|
| 1301 | } |
---|
| 1302 | |
---|
| 1303 | for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) { |
---|
| 1304 | int iRow = row[startColumn[iColumn]+length[iColumn]-1]; |
---|
| 1305 | int which = iRow-numberStaticRows_; |
---|
| 1306 | int jColumn = id_[iColumn-firstDynamic_]; |
---|
| 1307 | if (model->getStatus(iColumn)==ClpSimplex::atLowerBound|| |
---|
| 1308 | model->getStatus(iColumn)==ClpSimplex::atUpperBound) { |
---|
| 1309 | double value = solution[iColumn]; |
---|
[330] | 1310 | bool toLowerBound=true; |
---|
[328] | 1311 | if (columnUpper_) { |
---|
| 1312 | if (!columnLower_) { |
---|
| 1313 | if (fabs(value-columnUpper_[jColumn])<fabs(value)) |
---|
| 1314 | toLowerBound=false; |
---|
| 1315 | } else if (fabs(value-columnUpper_[jColumn])<fabs(value-columnLower_[iColumn])) { |
---|
| 1316 | toLowerBound=false; |
---|
| 1317 | } |
---|
| 1318 | } |
---|
| 1319 | if (toLowerBound) { |
---|
| 1320 | // throw out to lower bound |
---|
| 1321 | if (columnLower_) { |
---|
| 1322 | setDynamicStatus(jColumn,atLowerBound); |
---|
| 1323 | // treat solution as if exactly at a bound |
---|
| 1324 | double value = columnLower[iColumn]; |
---|
| 1325 | objectiveChange += cost[iColumn]*value; |
---|
| 1326 | } else { |
---|
| 1327 | setDynamicStatus(jColumn,atLowerBound); |
---|
| 1328 | } |
---|
| 1329 | } else { |
---|
| 1330 | // throw out to upper bound |
---|
| 1331 | setDynamicStatus(jColumn,atUpperBound); |
---|
| 1332 | double value = columnUpper[iColumn]; |
---|
| 1333 | objectiveChange += cost[iColumn]*value; |
---|
| 1334 | } |
---|
[330] | 1335 | numberDeleted++; |
---|
[328] | 1336 | } else { |
---|
| 1337 | assert(model->getStatus(iColumn)!=ClpSimplex::superBasic); // deal with later |
---|
| 1338 | int iPut = active[which]; |
---|
| 1339 | if (iPut>=0) { |
---|
| 1340 | // move |
---|
| 1341 | id_[firstAvailable_-firstDynamic_] = jColumn; |
---|
| 1342 | int numberThis = startColumn_[jColumn+1]-startColumn_[jColumn]+1; |
---|
| 1343 | length[firstAvailable_]=numberThis; |
---|
| 1344 | cost[firstAvailable_]=cost[iColumn]; |
---|
| 1345 | columnLower[firstAvailable_]=columnLower[iColumn]; |
---|
| 1346 | columnUpper[firstAvailable_]=columnUpper[iColumn]; |
---|
| 1347 | model->nonLinearCost()->setOne(firstAvailable_,solution[iColumn],0.0,COIN_DBL_MAX, |
---|
| 1348 | cost_[jColumn]); |
---|
| 1349 | CoinBigIndex base = startColumn_[jColumn]; |
---|
| 1350 | for (int j=0;j<numberThis-1;j++) { |
---|
| 1351 | row[numberElements]=row_[base+j]; |
---|
| 1352 | element[numberElements++]=element_[base+j]; |
---|
| 1353 | } |
---|
| 1354 | row[numberElements]=iPut+numberStaticRows_; |
---|
| 1355 | element[numberElements++]=1.0; |
---|
| 1356 | model->setStatus(firstAvailable_,model->getStatus(iColumn)); |
---|
| 1357 | solution[firstAvailable_] = solution[iColumn]; |
---|
| 1358 | firstAvailable_++; |
---|
| 1359 | startColumn[firstAvailable_]=numberElements; |
---|
| 1360 | } |
---|
| 1361 | } |
---|
| 1362 | } |
---|
| 1363 | // zero solution |
---|
| 1364 | CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_); |
---|
| 1365 | // zero cost |
---|
| 1366 | CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_); |
---|
| 1367 | // zero lengths |
---|
| 1368 | CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_); |
---|
| 1369 | for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) { |
---|
| 1370 | model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0); |
---|
| 1371 | model->setStatus(iColumn,ClpSimplex::atLowerBound); |
---|
| 1372 | columnLower[iColumn]=0.0; |
---|
| 1373 | columnUpper[iColumn]=COIN_DBL_MAX; |
---|
| 1374 | } |
---|
| 1375 | // move rhs and set rest to infinite |
---|
| 1376 | numberActiveSets_=0; |
---|
| 1377 | for (i=0;i<currentNumberActiveSets;i++) { |
---|
| 1378 | int iSet = fromIndex_[i]; |
---|
| 1379 | assert (toIndex_[iSet]==i); |
---|
| 1380 | int iRow = i+numberStaticRows_; |
---|
| 1381 | int iPut=active[i]; |
---|
| 1382 | if (iPut>=0) { |
---|
| 1383 | int in = iRow+numberColumns; |
---|
| 1384 | int out = iPut+numberColumns+numberStaticRows_; |
---|
| 1385 | solution[out]=solution[in]; |
---|
| 1386 | columnLower[out]=columnLower[in]; |
---|
| 1387 | columnUpper[out]=columnUpper[in]; |
---|
| 1388 | cost[out]=cost[in]; |
---|
[341] | 1389 | double trueLower; |
---|
| 1390 | if (lowerSet_[iSet]>-1.0e20) |
---|
| 1391 | trueLower = lowerSet_[iSet]; |
---|
| 1392 | else |
---|
| 1393 | trueLower=-COIN_DBL_MAX; |
---|
| 1394 | double trueUpper; |
---|
| 1395 | if (upperSet_[iSet]<1.0e20) |
---|
| 1396 | trueUpper = upperSet_[iSet]; |
---|
| 1397 | else |
---|
| 1398 | trueUpper=COIN_DBL_MAX; |
---|
[328] | 1399 | model->nonLinearCost()->setOne(out,solution[out], |
---|
[341] | 1400 | trueLower,trueUpper,0.0); |
---|
[328] | 1401 | model->setStatus(out,model->getStatus(in)); |
---|
| 1402 | toIndex_[iSet]=numberActiveSets_; |
---|
| 1403 | rhsOffset_[numberActiveSets_+numberStaticRows_]= rhsOffset_[i+numberStaticRows_]; |
---|
| 1404 | fromIndex_[numberActiveSets_++]=fromIndex_[i]; |
---|
| 1405 | } else { |
---|
| 1406 | // adjust offset |
---|
| 1407 | // put out as key |
---|
| 1408 | int jColumn = keyVariable_[iSet]; |
---|
| 1409 | toIndex_[iSet]=-1; |
---|
[330] | 1410 | if (jColumn<maximumGubColumns_) { |
---|
[328] | 1411 | setDynamicStatus(jColumn,soloKey); |
---|
| 1412 | double value = keyValue(iSet); |
---|
| 1413 | objectiveChange += cost_[jColumn]*value; |
---|
| 1414 | modifyOffset(jColumn,-value); |
---|
| 1415 | } |
---|
| 1416 | } |
---|
| 1417 | } |
---|
| 1418 | model->setObjectiveOffset(model->objectiveOffset()-objectiveChange); |
---|
| 1419 | delete [] active; |
---|
| 1420 | for (i=numberActiveSets_;i<currentNumberActiveSets;i++) { |
---|
| 1421 | int iSequence = i+numberStaticRows_+numberColumns; |
---|
| 1422 | solution[iSequence]=0.0; |
---|
| 1423 | columnLower[iSequence]=-COIN_DBL_MAX; |
---|
| 1424 | columnUpper[iSequence]=COIN_DBL_MAX; |
---|
| 1425 | cost[iSequence]=0.0; |
---|
| 1426 | model->nonLinearCost()->setOne(iSequence,solution[iSequence], |
---|
| 1427 | columnLower[iSequence], |
---|
| 1428 | columnUpper[iSequence],0.0); |
---|
| 1429 | model->setStatus(iSequence,ClpSimplex::basic); |
---|
| 1430 | rhsOffset_[i+numberStaticRows_]=0.0; |
---|
| 1431 | } |
---|
[330] | 1432 | if (currentNumberActiveSets!=numberActiveSets_||numberDeleted) { |
---|
[328] | 1433 | // deal with pivotVariable |
---|
| 1434 | int * pivotVariable = model->pivotVariable(); |
---|
| 1435 | int sequence=firstDynamic_; |
---|
| 1436 | int iRow=0; |
---|
| 1437 | int base1 = firstDynamic_; |
---|
| 1438 | int base2 = lastDynamic_; |
---|
| 1439 | int base3 = numberColumns+numberStaticRows_; |
---|
| 1440 | int numberRows = numberStaticRows_+currentNumberActiveSets; |
---|
| 1441 | while (sequence<firstAvailable_) { |
---|
| 1442 | int iPivot=pivotVariable[iRow]; |
---|
| 1443 | while (iPivot<base1||(iPivot>=base2&&iPivot<base3)) { |
---|
| 1444 | iPivot= pivotVariable[++iRow]; |
---|
| 1445 | } |
---|
| 1446 | pivotVariable[iRow++]=sequence; |
---|
| 1447 | sequence++; |
---|
| 1448 | } |
---|
| 1449 | // move normal basic ones down |
---|
| 1450 | int iPut=iRow; |
---|
| 1451 | for (;iRow<numberRows;iRow++) { |
---|
| 1452 | int iPivot=pivotVariable[iRow]; |
---|
| 1453 | if (iPivot<base1||(iPivot>=base2&&iPivot<base3)) { |
---|
| 1454 | pivotVariable[iPut++]=iPivot; |
---|
| 1455 | } |
---|
| 1456 | } |
---|
| 1457 | // and basic others |
---|
| 1458 | for (i=0;i<numberActiveSets_;i++) { |
---|
| 1459 | if (model->getRowStatus(i+numberStaticRows_)==ClpSimplex::basic) { |
---|
| 1460 | pivotVariable[iPut++]=i+base3; |
---|
| 1461 | } |
---|
| 1462 | } |
---|
| 1463 | for (i=numberActiveSets_;i<currentNumberActiveSets;i++) { |
---|
| 1464 | pivotVariable[iPut++]=i+base3; |
---|
| 1465 | } |
---|
| 1466 | assert (iPut==numberRows); |
---|
| 1467 | } |
---|
| 1468 | #ifdef CLP_DEBUG |
---|
[330] | 1469 | #if 0 |
---|
| 1470 | printf("row for set 244 is %d, row status %d value %g ",toIndex_[244],status_[244], |
---|
| 1471 | keyValue(244)); |
---|
| 1472 | int jj=startSet_[244]; |
---|
| 1473 | while (jj>=0) { |
---|
| 1474 | printf("var %d status %d ",jj,dynamicStatus_[jj]); |
---|
| 1475 | jj=next_[jj]; |
---|
| 1476 | } |
---|
| 1477 | printf("\n"); |
---|
| 1478 | #endif |
---|
[332] | 1479 | #ifdef CLP_DEBUG |
---|
| 1480 | { |
---|
| 1481 | // debug coding to analyze set |
---|
| 1482 | int i; |
---|
| 1483 | int kSet=-6; |
---|
| 1484 | if (kSet>=0) { |
---|
| 1485 | int * back = new int [numberGubColumns_]; |
---|
| 1486 | for (i=0;i<numberGubColumns_;i++) |
---|
| 1487 | back[i]=-1; |
---|
| 1488 | for (i=firstDynamic_;i<firstAvailable_;i++) |
---|
| 1489 | back[id_[i-firstDynamic_]]= i; |
---|
| 1490 | int sequence = startSet_[kSet]; |
---|
| 1491 | if (toIndex_[kSet]<0){ |
---|
| 1492 | printf("Not in - Set %d - slack status %d, key %d\n",kSet,status_[kSet],keyVariable_[kSet]); |
---|
| 1493 | while (sequence>=0) { |
---|
| 1494 | printf("( %d status %d ) ",sequence,getDynamicStatus(sequence)); |
---|
| 1495 | sequence = next_[sequence]; |
---|
| 1496 | } |
---|
| 1497 | } else { |
---|
| 1498 | int iRow=numberStaticRows_+toIndex_[kSet]; |
---|
| 1499 | printf("In - Set %d - slack status %d, key %d offset %g slack %g\n",kSet,status_[kSet],keyVariable_[kSet], |
---|
| 1500 | rhsOffset_[iRow],model->solutionRegion(0)[iRow]); |
---|
| 1501 | while (sequence>=0) { |
---|
| 1502 | int iBack = back[sequence]; |
---|
| 1503 | printf("( %d status %d value %g) ",sequence,getDynamicStatus(sequence),model->solutionRegion()[iBack]); |
---|
| 1504 | sequence = next_[sequence]; |
---|
| 1505 | } |
---|
| 1506 | } |
---|
| 1507 | printf("\n"); |
---|
| 1508 | delete [] back; |
---|
| 1509 | } |
---|
| 1510 | } |
---|
| 1511 | #endif |
---|
[328] | 1512 | int n=numberActiveSets_; |
---|
| 1513 | for (i=0;i<numberSets_;i++) { |
---|
| 1514 | if (toIndex_[i]<0) { |
---|
[344] | 1515 | //assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]); |
---|
[328] | 1516 | n++; |
---|
| 1517 | } |
---|
| 1518 | } |
---|
| 1519 | assert (n==numberSets_); |
---|
| 1520 | #endif |
---|
| 1521 | return 1; |
---|
| 1522 | } |
---|
| 1523 | void |
---|
| 1524 | ClpDynamicMatrix::times(double scalar, |
---|
| 1525 | const double * x, double * y) const |
---|
| 1526 | { |
---|
| 1527 | if (model_->specialOptions()!=16) { |
---|
| 1528 | ClpPackedMatrix::times(scalar,x,y); |
---|
| 1529 | } else { |
---|
| 1530 | int iRow; |
---|
| 1531 | const double * element = matrix_->getElements(); |
---|
| 1532 | const int * row = matrix_->getIndices(); |
---|
| 1533 | const CoinBigIndex * startColumn = matrix_->getVectorStarts(); |
---|
| 1534 | const int * length = matrix_->getVectorLengths(); |
---|
| 1535 | int * pivotVariable = model_->pivotVariable(); |
---|
| 1536 | for (iRow=0;iRow<numberStaticRows_+numberActiveSets_;iRow++) { |
---|
| 1537 | y[iRow] -= scalar*rhsOffset_[iRow]; |
---|
| 1538 | int iColumn = pivotVariable[iRow]; |
---|
| 1539 | if (iColumn<lastDynamic_) { |
---|
| 1540 | CoinBigIndex j; |
---|
| 1541 | double value = scalar*x[iColumn]; |
---|
| 1542 | if (value) { |
---|
| 1543 | for (j=startColumn[iColumn]; |
---|
| 1544 | j<startColumn[iColumn]+length[iColumn];j++) { |
---|
| 1545 | int jRow=row[j]; |
---|
| 1546 | y[jRow] += value*element[j]; |
---|
| 1547 | } |
---|
| 1548 | } |
---|
| 1549 | } |
---|
| 1550 | } |
---|
| 1551 | } |
---|
| 1552 | } |
---|
| 1553 | // Modifies rhs offset |
---|
| 1554 | void |
---|
| 1555 | ClpDynamicMatrix::modifyOffset(int sequence, double amount) |
---|
| 1556 | { |
---|
| 1557 | if (amount) { |
---|
| 1558 | assert (rhsOffset_); |
---|
| 1559 | CoinBigIndex j; |
---|
| 1560 | for (j=startColumn_[sequence];j<startColumn_[sequence+1];j++) { |
---|
| 1561 | int iRow = row_[j]; |
---|
| 1562 | rhsOffset_[iRow] += amount*element_[j]; |
---|
| 1563 | } |
---|
| 1564 | } |
---|
| 1565 | } |
---|
| 1566 | // Gets key value when none in small |
---|
| 1567 | double |
---|
| 1568 | ClpDynamicMatrix::keyValue(int iSet) const |
---|
| 1569 | { |
---|
| 1570 | double value=0.0; |
---|
| 1571 | if (toIndex_[iSet]<0) { |
---|
| 1572 | int key = keyVariable_[iSet]; |
---|
[330] | 1573 | if (key<maximumGubColumns_) { |
---|
[328] | 1574 | if (getStatus(iSet)==ClpSimplex::atLowerBound) |
---|
| 1575 | value=lowerSet_[iSet]; |
---|
| 1576 | else |
---|
| 1577 | value=upperSet_[iSet]; |
---|
| 1578 | int numberKey=0; |
---|
[329] | 1579 | int j= startSet_[iSet]; |
---|
| 1580 | while (j>=0) { |
---|
[328] | 1581 | DynamicStatus status = getDynamicStatus(j); |
---|
| 1582 | assert (status!=inSmall); |
---|
| 1583 | if (status==soloKey) { |
---|
| 1584 | numberKey++; |
---|
| 1585 | } else if (status==atUpperBound) { |
---|
| 1586 | value -= columnUpper_[j]; |
---|
| 1587 | } else if (columnLower_) { |
---|
| 1588 | value -= columnLower_[j]; |
---|
| 1589 | } |
---|
[329] | 1590 | j = next_[j]; //onto next in set |
---|
[328] | 1591 | } |
---|
| 1592 | assert (numberKey==1); |
---|
| 1593 | } else { |
---|
[329] | 1594 | int j= startSet_[iSet]; |
---|
| 1595 | while (j>=0) { |
---|
[328] | 1596 | DynamicStatus status = getDynamicStatus(j); |
---|
| 1597 | assert (status!=inSmall); |
---|
| 1598 | assert (status!=soloKey); |
---|
| 1599 | if (status==atUpperBound) { |
---|
| 1600 | value += columnUpper_[j]; |
---|
| 1601 | } else if (columnLower_) { |
---|
| 1602 | value += columnLower_[j]; |
---|
| 1603 | } |
---|
[329] | 1604 | j = next_[j]; //onto next in set |
---|
[328] | 1605 | } |
---|
| 1606 | } |
---|
| 1607 | } |
---|
| 1608 | return value; |
---|
| 1609 | } |
---|
| 1610 | // Switches off dj checking each factorization (for BIG models) |
---|
| 1611 | void |
---|
| 1612 | ClpDynamicMatrix::switchOffCheck() |
---|
| 1613 | { |
---|
| 1614 | noCheck_=0; |
---|
| 1615 | infeasibilityWeight_=0.0; |
---|
| 1616 | } |
---|
| 1617 | /* Creates a variable. This is called after partial pricing and may modify matrix. |
---|
| 1618 | May update bestSequence. |
---|
| 1619 | */ |
---|
| 1620 | void |
---|
| 1621 | ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence) |
---|
| 1622 | { |
---|
| 1623 | int numberRows = model->numberRows(); |
---|
| 1624 | int slackOffset = lastDynamic_+numberRows; |
---|
| 1625 | int structuralOffset = slackOffset+numberSets_; |
---|
| 1626 | int bestSequence2=savedBestSequence_-structuralOffset; |
---|
| 1627 | if (bestSequence>=slackOffset) { |
---|
| 1628 | double * columnLower = model->lowerRegion(); |
---|
| 1629 | double * columnUpper = model->upperRegion(); |
---|
| 1630 | double * solution = model->solutionRegion(); |
---|
| 1631 | double * reducedCost = model->djRegion(); |
---|
| 1632 | const double * duals = model->dualRowSolution(); |
---|
| 1633 | if (toIndex_[savedBestSet_]<0) { |
---|
| 1634 | // need to put key into basis |
---|
| 1635 | int newRow = numberActiveSets_+numberStaticRows_; |
---|
| 1636 | model->dualRowSolution()[newRow]=savedBestGubDual_; |
---|
| 1637 | double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set |
---|
| 1638 | toIndex_[savedBestSet_]=numberActiveSets_; |
---|
| 1639 | fromIndex_[numberActiveSets_++]=savedBestSet_; |
---|
| 1640 | int iSequence=lastDynamic_+newRow; |
---|
| 1641 | // we need to get lower and upper correct |
---|
| 1642 | double shift=0.0; |
---|
[329] | 1643 | int j= startSet_[savedBestSet_]; |
---|
| 1644 | while (j>=0) { |
---|
[328] | 1645 | if (getDynamicStatus(j)==atUpperBound) |
---|
| 1646 | shift += columnUpper_[j]; |
---|
| 1647 | else if (getDynamicStatus(j)==atLowerBound&&columnLower_) |
---|
| 1648 | shift += columnLower_[j]; |
---|
[329] | 1649 | j = next_[j]; //onto next in set |
---|
[328] | 1650 | } |
---|
| 1651 | if (lowerSet_[savedBestSet_]>-1.0e20) |
---|
| 1652 | columnLower[iSequence] = lowerSet_[savedBestSet_]; |
---|
| 1653 | else |
---|
| 1654 | columnLower[iSequence]=-COIN_DBL_MAX; |
---|
| 1655 | if (upperSet_[savedBestSet_]<1.0e20) |
---|
| 1656 | columnUpper[iSequence] = upperSet_[savedBestSet_]; |
---|
| 1657 | else |
---|
| 1658 | columnUpper[iSequence]=COIN_DBL_MAX; |
---|
| 1659 | #ifdef CLP_DEBUG |
---|
| 1660 | if (model_->logLevel()==63) { |
---|
| 1661 | printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n", |
---|
| 1662 | bestSequence2,savedBestSet_,keyVariable_[savedBestSet_], |
---|
| 1663 | columnLower[iSequence],columnUpper[iSequence],valueOfKey); |
---|
[329] | 1664 | int j= startSet_[savedBestSet_]; |
---|
| 1665 | while (j>=0) { |
---|
[328] | 1666 | if (getDynamicStatus(j)==atUpperBound) |
---|
| 1667 | printf("%d atup ",j); |
---|
| 1668 | else if (getDynamicStatus(j)==atLowerBound) |
---|
| 1669 | printf("%d atlo ",j); |
---|
| 1670 | else if (getDynamicStatus(j)==soloKey) |
---|
| 1671 | printf("%d solo ",j); |
---|
| 1672 | else |
---|
| 1673 | abort(); |
---|
[329] | 1674 | j = next_[j]; //onto next in set |
---|
[328] | 1675 | } |
---|
| 1676 | printf("\n"); |
---|
| 1677 | } |
---|
| 1678 | #endif |
---|
[330] | 1679 | if (keyVariable_[savedBestSet_]<maximumGubColumns_) { |
---|
[328] | 1680 | // slack not key |
---|
| 1681 | model_->pivotVariable()[newRow]=firstAvailable_; |
---|
| 1682 | backToPivotRow_[firstAvailable_]=newRow; |
---|
| 1683 | model->setStatus(iSequence,getStatus(savedBestSet_)); |
---|
| 1684 | model->djRegion()[iSequence] = savedBestGubDual_; |
---|
| 1685 | if (getStatus(savedBestSet_)==ClpSimplex::atUpperBound) |
---|
| 1686 | solution[iSequence] = columnUpper[iSequence]; |
---|
| 1687 | else |
---|
| 1688 | solution[iSequence] = columnLower[iSequence]; |
---|
| 1689 | // create variable and pivot in |
---|
| 1690 | int key=keyVariable_[savedBestSet_]; |
---|
| 1691 | setDynamicStatus(key,inSmall); |
---|
| 1692 | double * element = matrix_->getMutableElements(); |
---|
| 1693 | int * row = matrix_->getMutableIndices(); |
---|
| 1694 | CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1695 | int * length = matrix_->getMutableVectorLengths(); |
---|
| 1696 | CoinBigIndex numberElements = startColumn[firstAvailable_]; |
---|
| 1697 | int numberThis = startColumn_[key+1]-startColumn_[key]+1; |
---|
| 1698 | if (numberElements+numberThis>numberElements_) { |
---|
| 1699 | // need to redo |
---|
[399] | 1700 | numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis); |
---|
[328] | 1701 | matrix_->reserve(lastDynamic_,numberElements_); |
---|
| 1702 | element = matrix_->getMutableElements(); |
---|
| 1703 | row = matrix_->getMutableIndices(); |
---|
| 1704 | // these probably okay but be safe |
---|
| 1705 | startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1706 | length = matrix_->getMutableVectorLengths(); |
---|
| 1707 | } |
---|
| 1708 | // already set startColumn[firstAvailable_]=numberElements; |
---|
| 1709 | length[firstAvailable_]=numberThis; |
---|
| 1710 | model->costRegion()[firstAvailable_]=cost_[key]; |
---|
| 1711 | CoinBigIndex base = startColumn_[key]; |
---|
| 1712 | for (int j=0;j<numberThis-1;j++) { |
---|
| 1713 | row[numberElements]=row_[base+j]; |
---|
| 1714 | element[numberElements++]=element_[base+j]; |
---|
| 1715 | } |
---|
| 1716 | row[numberElements]=newRow; |
---|
| 1717 | element[numberElements++]=1.0; |
---|
| 1718 | id_[firstAvailable_-firstDynamic_]=key; |
---|
| 1719 | model->setObjectiveOffset(model->objectiveOffset()+cost_[key]* |
---|
| 1720 | valueOfKey); |
---|
| 1721 | model->solutionRegion()[firstAvailable_]= valueOfKey; |
---|
| 1722 | model->setStatus(firstAvailable_,ClpSimplex::basic); |
---|
| 1723 | // ***** need to adjust effective rhs |
---|
| 1724 | if (!columnLower_&&!columnUpper_) { |
---|
| 1725 | columnLower[firstAvailable_] = 0.0; |
---|
| 1726 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 1727 | } else { |
---|
| 1728 | if (columnLower_) |
---|
| 1729 | columnLower[firstAvailable_] = columnLower_[key]; |
---|
| 1730 | else |
---|
| 1731 | columnLower[firstAvailable_] = 0.0; |
---|
| 1732 | if (columnUpper_) |
---|
| 1733 | columnUpper[firstAvailable_] = columnUpper_[key]; |
---|
| 1734 | else |
---|
| 1735 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 1736 | } |
---|
| 1737 | model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_], |
---|
| 1738 | columnLower[firstAvailable_], |
---|
| 1739 | columnUpper[firstAvailable_],cost_[key]); |
---|
| 1740 | startColumn[firstAvailable_+1]=numberElements; |
---|
| 1741 | reducedCost[firstAvailable_] = 0.0; |
---|
| 1742 | modifyOffset(key,valueOfKey); |
---|
| 1743 | rhsOffset_[newRow] = -shift; // sign? |
---|
[336] | 1744 | #ifdef CLP_DEBUG |
---|
| 1745 | model->rowArray(1)->checkClear(); |
---|
| 1746 | #endif |
---|
[328] | 1747 | // now pivot in |
---|
[336] | 1748 | unpack(model,model->rowArray(1),firstAvailable_); |
---|
| 1749 | model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(1)); |
---|
| 1750 | double alpha = model->rowArray(1)->denseVector()[newRow]; |
---|
[332] | 1751 | int updateStatus = |
---|
[336] | 1752 | model->factorization()->replaceColumn(model, |
---|
| 1753 | model->rowArray(2), |
---|
| 1754 | model->rowArray(1), |
---|
| 1755 | newRow, alpha); |
---|
| 1756 | model->rowArray(1)->clear(); |
---|
[335] | 1757 | if (updateStatus) { |
---|
| 1758 | if (updateStatus==3) { |
---|
| 1759 | // out of memory |
---|
| 1760 | // increase space if not many iterations |
---|
| 1761 | if (model->factorization()->pivots()< |
---|
| 1762 | 0.5*model->factorization()->maximumPivots()&& |
---|
| 1763 | model->factorization()->pivots()<400) |
---|
| 1764 | model->factorization()->areaFactor( |
---|
| 1765 | model->factorization()->areaFactor() * 1.1); |
---|
| 1766 | } else { |
---|
| 1767 | printf("Bad returncode %d from replaceColumn\n",updateStatus); |
---|
| 1768 | } |
---|
| 1769 | bestSequence=-1; |
---|
| 1770 | return; |
---|
| 1771 | } |
---|
[328] | 1772 | // firstAvailable_ only finally updated if good pivot (in updatePivot) |
---|
| 1773 | // otherwise it reverts to firstAvailableBefore_ |
---|
| 1774 | firstAvailable_++; |
---|
| 1775 | } else { |
---|
| 1776 | // slack key |
---|
| 1777 | model->setStatus(iSequence,ClpSimplex::basic); |
---|
| 1778 | model->djRegion()[iSequence]=0.0; |
---|
| 1779 | solution[iSequence]=shift; |
---|
| 1780 | rhsOffset_[newRow] = -shift; // sign? |
---|
| 1781 | } |
---|
| 1782 | // correct slack |
---|
| 1783 | model->costRegion()[iSequence] = 0.0; |
---|
| 1784 | model->nonLinearCost()->setOne(iSequence,solution[iSequence],columnLower[iSequence], |
---|
| 1785 | columnUpper[iSequence],0.0); |
---|
| 1786 | } |
---|
| 1787 | if (savedBestSequence_>=structuralOffset) { |
---|
| 1788 | // recompute dj and create |
---|
| 1789 | double value=cost_[bestSequence2]-savedBestGubDual_; |
---|
[604] | 1790 | for (CoinBigIndex jBigIndex=startColumn_[bestSequence2]; |
---|
| 1791 | jBigIndex<startColumn_[bestSequence2+1];jBigIndex++) { |
---|
| 1792 | int jRow=row_[jBigIndex]; |
---|
| 1793 | value -= duals[jRow]*element_[jBigIndex]; |
---|
[328] | 1794 | } |
---|
| 1795 | int gubRow = toIndex_[savedBestSet_] + numberStaticRows_; |
---|
| 1796 | double * element = matrix_->getMutableElements(); |
---|
| 1797 | int * row = matrix_->getMutableIndices(); |
---|
| 1798 | CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1799 | int * length = matrix_->getMutableVectorLengths(); |
---|
| 1800 | CoinBigIndex numberElements = startColumn[firstAvailable_]; |
---|
| 1801 | int numberThis = startColumn_[bestSequence2+1]-startColumn_[bestSequence2]+1; |
---|
| 1802 | if (numberElements+numberThis>numberElements_) { |
---|
| 1803 | // need to redo |
---|
[399] | 1804 | numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis); |
---|
[328] | 1805 | matrix_->reserve(lastDynamic_,numberElements_); |
---|
| 1806 | element = matrix_->getMutableElements(); |
---|
| 1807 | row = matrix_->getMutableIndices(); |
---|
| 1808 | // these probably okay but be safe |
---|
| 1809 | startColumn = matrix_->getMutableVectorStarts(); |
---|
| 1810 | length = matrix_->getMutableVectorLengths(); |
---|
| 1811 | } |
---|
| 1812 | // already set startColumn[firstAvailable_]=numberElements; |
---|
| 1813 | length[firstAvailable_]=numberThis; |
---|
| 1814 | model->costRegion()[firstAvailable_]=cost_[bestSequence2]; |
---|
| 1815 | CoinBigIndex base = startColumn_[bestSequence2]; |
---|
| 1816 | for (int j=0;j<numberThis-1;j++) { |
---|
| 1817 | row[numberElements]=row_[base+j]; |
---|
| 1818 | element[numberElements++]=element_[base+j]; |
---|
| 1819 | } |
---|
| 1820 | row[numberElements]=gubRow; |
---|
| 1821 | element[numberElements++]=1.0; |
---|
| 1822 | id_[firstAvailable_-firstDynamic_]=bestSequence2; |
---|
| 1823 | //printf("best %d\n",bestSequence2); |
---|
| 1824 | model->solutionRegion()[firstAvailable_]=0.0; |
---|
| 1825 | if (!columnLower_&&!columnUpper_) { |
---|
| 1826 | model->setStatus(firstAvailable_,ClpSimplex::atLowerBound); |
---|
| 1827 | columnLower[firstAvailable_] = 0.0; |
---|
| 1828 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 1829 | } else { |
---|
| 1830 | DynamicStatus status = getDynamicStatus(bestSequence2); |
---|
| 1831 | if (columnLower_) |
---|
| 1832 | columnLower[firstAvailable_] = columnLower_[bestSequence2]; |
---|
| 1833 | else |
---|
| 1834 | columnLower[firstAvailable_] = 0.0; |
---|
| 1835 | if (columnUpper_) |
---|
| 1836 | columnUpper[firstAvailable_] = columnUpper_[bestSequence2]; |
---|
| 1837 | else |
---|
| 1838 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 1839 | if (status==atLowerBound) { |
---|
| 1840 | solution[firstAvailable_]=columnLower[firstAvailable_]; |
---|
| 1841 | model->setStatus(firstAvailable_,ClpSimplex::atLowerBound); |
---|
| 1842 | } else { |
---|
| 1843 | solution[firstAvailable_]=columnUpper[firstAvailable_]; |
---|
| 1844 | model->setStatus(firstAvailable_,ClpSimplex::atUpperBound); |
---|
| 1845 | } |
---|
| 1846 | } |
---|
| 1847 | model->setObjectiveOffset(model->objectiveOffset()+cost_[bestSequence2]* |
---|
| 1848 | solution[firstAvailable_]); |
---|
| 1849 | model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_], |
---|
| 1850 | columnLower[firstAvailable_], |
---|
| 1851 | columnUpper[firstAvailable_],cost_[bestSequence2]); |
---|
| 1852 | bestSequence=firstAvailable_; |
---|
| 1853 | // firstAvailable_ only updated if good pivot (in updatePivot) |
---|
| 1854 | startColumn[firstAvailable_+1]=numberElements; |
---|
| 1855 | //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_); |
---|
| 1856 | reducedCost[bestSequence] = value; |
---|
| 1857 | } else { |
---|
| 1858 | // slack - row must just have been created |
---|
| 1859 | assert (toIndex_[savedBestSet_]==numberActiveSets_-1); |
---|
| 1860 | int newRow = numberStaticRows_+numberActiveSets_-1; |
---|
| 1861 | bestSequence = lastDynamic_ + newRow; |
---|
| 1862 | reducedCost[bestSequence] = savedBestGubDual_; |
---|
| 1863 | } |
---|
| 1864 | } |
---|
| 1865 | // clear for next iteration |
---|
| 1866 | savedBestSequence_=-1; |
---|
| 1867 | } |
---|
| 1868 | // Returns reduced cost of a variable |
---|
| 1869 | double |
---|
| 1870 | ClpDynamicMatrix::reducedCost(ClpSimplex * model,int sequence) const |
---|
| 1871 | { |
---|
| 1872 | int numberRows = model->numberRows(); |
---|
| 1873 | int slackOffset = lastDynamic_+numberRows; |
---|
| 1874 | if (sequence<slackOffset) |
---|
| 1875 | return model->djRegion()[sequence]; |
---|
| 1876 | else |
---|
| 1877 | return savedBestDj_; |
---|
| 1878 | } |
---|
| 1879 | // Does gub crash |
---|
| 1880 | void |
---|
| 1881 | ClpDynamicMatrix::gubCrash() |
---|
| 1882 | { |
---|
| 1883 | // Do basis - cheapest or slack if feasible |
---|
| 1884 | int longestSet=0; |
---|
| 1885 | int iSet; |
---|
[329] | 1886 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 1887 | int n=0; |
---|
| 1888 | int j= startSet_[iSet]; |
---|
| 1889 | while (j>=0) { |
---|
| 1890 | n++; |
---|
| 1891 | j=next_[j]; |
---|
| 1892 | } |
---|
[399] | 1893 | longestSet = CoinMax(longestSet,n); |
---|
[329] | 1894 | } |
---|
[328] | 1895 | double * upper = new double[longestSet+1]; |
---|
| 1896 | double * cost = new double[longestSet+1]; |
---|
| 1897 | double * lower = new double[longestSet+1]; |
---|
| 1898 | double * solution = new double[longestSet+1]; |
---|
[329] | 1899 | int * back = new int[longestSet+1]; |
---|
[328] | 1900 | double tolerance = model_->primalTolerance(); |
---|
| 1901 | double objectiveOffset = 0.0; |
---|
| 1902 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 1903 | int iBasic=-1; |
---|
| 1904 | double value=0.0; |
---|
| 1905 | // find cheapest |
---|
[329] | 1906 | int numberInSet = 0; |
---|
| 1907 | int j=startSet_[iSet]; |
---|
| 1908 | while (j>=0) { |
---|
| 1909 | if (!columnLower_) |
---|
| 1910 | lower[numberInSet]=0.0; |
---|
| 1911 | else |
---|
| 1912 | lower[numberInSet]= columnLower_[j]; |
---|
| 1913 | if (!columnUpper_) |
---|
| 1914 | upper[numberInSet]=COIN_DBL_MAX; |
---|
| 1915 | else |
---|
[344] | 1916 | upper[numberInSet]= columnUpper_[j]; |
---|
[329] | 1917 | back[numberInSet++]=j; |
---|
| 1918 | j = next_[j]; |
---|
[328] | 1919 | } |
---|
| 1920 | CoinFillN(solution,numberInSet,0.0); |
---|
| 1921 | // and slack |
---|
| 1922 | iBasic=numberInSet; |
---|
| 1923 | solution[iBasic]=-value; |
---|
| 1924 | lower[iBasic]=-upperSet_[iSet]; |
---|
| 1925 | upper[iBasic]=-lowerSet_[iSet]; |
---|
| 1926 | int kphase; |
---|
| 1927 | if (value<lowerSet_[iSet]-tolerance||value>upperSet_[iSet]+tolerance) { |
---|
| 1928 | // infeasible |
---|
| 1929 | kphase=0; |
---|
| 1930 | // remember bounds are flipped so opposite to natural |
---|
| 1931 | if (value<lowerSet_[iSet]-tolerance) |
---|
| 1932 | cost[iBasic]=1.0; |
---|
| 1933 | else |
---|
| 1934 | cost[iBasic]=-1.0; |
---|
| 1935 | CoinZeroN(cost,numberInSet); |
---|
| 1936 | double dualTolerance =model_->dualTolerance(); |
---|
| 1937 | for (int iphase =kphase;iphase<2;iphase++) { |
---|
| 1938 | if (iphase) { |
---|
| 1939 | cost[numberInSet]=0.0; |
---|
| 1940 | for (int j=0;j<numberInSet;j++) |
---|
[329] | 1941 | cost[j]= cost_[back[j]]; |
---|
[328] | 1942 | } |
---|
| 1943 | // now do one row lp |
---|
| 1944 | bool improve=true; |
---|
| 1945 | while (improve) { |
---|
| 1946 | improve=false; |
---|
| 1947 | double dual = cost[iBasic]; |
---|
| 1948 | int chosen =-1; |
---|
| 1949 | double best=dualTolerance; |
---|
| 1950 | int way=0; |
---|
| 1951 | for (int i=0;i<=numberInSet;i++) { |
---|
| 1952 | double dj = cost[i]-dual; |
---|
| 1953 | double improvement =0.0; |
---|
| 1954 | double distance=0.0; |
---|
| 1955 | if (iphase||i<numberInSet) |
---|
| 1956 | assert (solution[i]>=lower[i]&&solution[i]<=upper[i]); |
---|
| 1957 | if (dj>dualTolerance) |
---|
| 1958 | improvement = dj*(solution[i]-lower[i]); |
---|
| 1959 | else if (dj<-dualTolerance) |
---|
| 1960 | improvement = dj*(solution[i]-upper[i]); |
---|
| 1961 | if (improvement>best) { |
---|
| 1962 | best=improvement; |
---|
| 1963 | chosen=i; |
---|
| 1964 | if (dj<0.0) { |
---|
| 1965 | way = 1; |
---|
| 1966 | distance = upper[i]-solution[i]; |
---|
| 1967 | } else { |
---|
| 1968 | way = -1; |
---|
| 1969 | distance = solution[i]-lower[i]; |
---|
| 1970 | } |
---|
| 1971 | } |
---|
| 1972 | } |
---|
| 1973 | if (chosen>=0) { |
---|
| 1974 | improve=true; |
---|
| 1975 | // now see how far |
---|
| 1976 | if (way>0) { |
---|
| 1977 | // incoming increasing so basic decreasing |
---|
| 1978 | // if phase 0 then go to nearest bound |
---|
| 1979 | double distance=upper[chosen]-solution[chosen]; |
---|
| 1980 | double basicDistance; |
---|
| 1981 | if (!iphase) { |
---|
| 1982 | assert (iBasic==numberInSet); |
---|
| 1983 | assert (solution[iBasic]>upper[iBasic]); |
---|
| 1984 | basicDistance = solution[iBasic]-upper[iBasic]; |
---|
| 1985 | } else { |
---|
| 1986 | basicDistance = solution[iBasic]-lower[iBasic]; |
---|
| 1987 | } |
---|
| 1988 | // need extra coding for unbounded |
---|
[399] | 1989 | assert (CoinMin(distance,basicDistance)<1.0e20); |
---|
[328] | 1990 | if (distance>basicDistance) { |
---|
| 1991 | // incoming becomes basic |
---|
| 1992 | solution[chosen] += basicDistance; |
---|
| 1993 | if (!iphase) |
---|
| 1994 | solution[iBasic]=upper[iBasic]; |
---|
| 1995 | else |
---|
| 1996 | solution[iBasic]=lower[iBasic]; |
---|
| 1997 | iBasic = chosen; |
---|
| 1998 | } else { |
---|
| 1999 | // flip |
---|
| 2000 | solution[chosen]=upper[chosen]; |
---|
| 2001 | solution[iBasic] -= distance; |
---|
| 2002 | } |
---|
| 2003 | } else { |
---|
| 2004 | // incoming decreasing so basic increasing |
---|
| 2005 | // if phase 0 then go to nearest bound |
---|
| 2006 | double distance=solution[chosen]-lower[chosen]; |
---|
| 2007 | double basicDistance; |
---|
| 2008 | if (!iphase) { |
---|
| 2009 | assert (iBasic==numberInSet); |
---|
| 2010 | assert (solution[iBasic]<lower[iBasic]); |
---|
| 2011 | basicDistance = lower[iBasic]-solution[iBasic]; |
---|
| 2012 | } else { |
---|
| 2013 | basicDistance = upper[iBasic]-solution[iBasic]; |
---|
| 2014 | } |
---|
| 2015 | // need extra coding for unbounded - for now just exit |
---|
[399] | 2016 | if (CoinMin(distance,basicDistance)>1.0e20) { |
---|
[328] | 2017 | printf("unbounded on set %d\n",iSet); |
---|
| 2018 | iphase=1; |
---|
| 2019 | iBasic=numberInSet; |
---|
| 2020 | break; |
---|
| 2021 | } |
---|
| 2022 | if (distance>basicDistance) { |
---|
| 2023 | // incoming becomes basic |
---|
| 2024 | solution[chosen] -= basicDistance; |
---|
| 2025 | if (!iphase) |
---|
| 2026 | solution[iBasic]=lower[iBasic]; |
---|
| 2027 | else |
---|
| 2028 | solution[iBasic]=upper[iBasic]; |
---|
| 2029 | iBasic = chosen; |
---|
| 2030 | } else { |
---|
| 2031 | // flip |
---|
| 2032 | solution[chosen]=lower[chosen]; |
---|
| 2033 | solution[iBasic] += distance; |
---|
| 2034 | } |
---|
| 2035 | } |
---|
| 2036 | if (!iphase) { |
---|
| 2037 | if(iBasic<numberInSet) |
---|
| 2038 | break; // feasible |
---|
| 2039 | else if (solution[iBasic]>=lower[iBasic]&& |
---|
| 2040 | solution[iBasic]<=upper[iBasic]) |
---|
| 2041 | break; // feasible (on flip) |
---|
| 2042 | } |
---|
| 2043 | } |
---|
| 2044 | } |
---|
| 2045 | } |
---|
| 2046 | } |
---|
| 2047 | // do solution i.e. bounds |
---|
| 2048 | if (columnLower_||columnUpper_) { |
---|
| 2049 | for (int j=0;j<numberInSet;j++) { |
---|
| 2050 | if (j!=iBasic) { |
---|
| 2051 | objectiveOffset += solution[j]*cost[j]; |
---|
| 2052 | if (columnLower_&&columnUpper_) { |
---|
[329] | 2053 | if (fabs(solution[j]-columnLower_[back[j]])> |
---|
| 2054 | fabs(solution[j]-columnUpper_[back[j]])) |
---|
| 2055 | setDynamicStatus(back[j],atUpperBound); |
---|
[328] | 2056 | } else if (columnUpper_&&solution[j]>0.0) { |
---|
[329] | 2057 | setDynamicStatus(back[j],atUpperBound); |
---|
[328] | 2058 | } else { |
---|
[329] | 2059 | setDynamicStatus(back[j],atLowerBound); |
---|
[330] | 2060 | assert(!solution[j]); |
---|
[328] | 2061 | } |
---|
| 2062 | } |
---|
| 2063 | } |
---|
| 2064 | } |
---|
| 2065 | // convert iBasic back and do bounds |
---|
| 2066 | if (iBasic==numberInSet) { |
---|
| 2067 | // slack basic |
---|
| 2068 | setStatus(iSet,ClpSimplex::basic); |
---|
[330] | 2069 | iBasic=iSet+maximumGubColumns_; |
---|
[328] | 2070 | } else { |
---|
[329] | 2071 | iBasic = back[iBasic]; |
---|
[328] | 2072 | setDynamicStatus(iBasic,soloKey); |
---|
| 2073 | // remember bounds flipped |
---|
| 2074 | if (upper[numberInSet]==lower[numberInSet]) |
---|
| 2075 | setStatus(iSet,ClpSimplex::isFixed); |
---|
| 2076 | else if (solution[numberInSet]==upper[numberInSet]) |
---|
| 2077 | setStatus(iSet,ClpSimplex::atLowerBound); |
---|
| 2078 | else if (solution[numberInSet]==lower[numberInSet]) |
---|
| 2079 | setStatus(iSet,ClpSimplex::atUpperBound); |
---|
| 2080 | else |
---|
| 2081 | abort(); |
---|
| 2082 | } |
---|
| 2083 | keyVariable_[iSet]=iBasic; |
---|
| 2084 | } |
---|
| 2085 | model_->setObjectiveOffset(objectiveOffset_-objectiveOffset); |
---|
| 2086 | delete [] lower; |
---|
| 2087 | delete [] solution; |
---|
| 2088 | delete [] upper; |
---|
| 2089 | delete [] cost; |
---|
[329] | 2090 | delete [] back; |
---|
[328] | 2091 | // make sure matrix is in good shape |
---|
| 2092 | matrix_->orderMatrix(); |
---|
| 2093 | } |
---|
| 2094 | // Populates initial matrix from dynamic status |
---|
| 2095 | void |
---|
| 2096 | ClpDynamicMatrix::initialProblem() |
---|
| 2097 | { |
---|
| 2098 | int iSet; |
---|
| 2099 | double * element = matrix_->getMutableElements(); |
---|
| 2100 | int * row = matrix_->getMutableIndices(); |
---|
| 2101 | CoinBigIndex * startColumn = matrix_->getMutableVectorStarts(); |
---|
| 2102 | int * length = matrix_->getMutableVectorLengths(); |
---|
| 2103 | double * cost = model_->objective(); |
---|
| 2104 | double * solution = model_->primalColumnSolution(); |
---|
| 2105 | double * columnLower = model_->columnLower(); |
---|
| 2106 | double * columnUpper = model_->columnUpper(); |
---|
| 2107 | double * rowSolution = model_->primalRowSolution(); |
---|
| 2108 | double * rowLower = model_->rowLower(); |
---|
| 2109 | double * rowUpper = model_->rowUpper(); |
---|
| 2110 | CoinBigIndex numberElements=startColumn[firstDynamic_]; |
---|
| 2111 | |
---|
| 2112 | firstAvailable_ = firstDynamic_; |
---|
| 2113 | numberActiveSets_=0; |
---|
| 2114 | for (iSet=0;iSet<numberSets_;iSet++) { |
---|
| 2115 | toIndex_[iSet]=-1; |
---|
| 2116 | int numberActive=0; |
---|
[329] | 2117 | int whichKey=-1; |
---|
| 2118 | if (getStatus(iSet)==ClpSimplex::basic) |
---|
[330] | 2119 | whichKey = maximumGubColumns_+iSet; |
---|
[329] | 2120 | else |
---|
| 2121 | whichKey = -1; |
---|
| 2122 | int j= startSet_[iSet]; |
---|
| 2123 | while (j>=0) { |
---|
| 2124 | assert (getDynamicStatus(j)!=soloKey||whichKey==-1); |
---|
[328] | 2125 | if (getDynamicStatus(j)==inSmall) |
---|
| 2126 | numberActive++; |
---|
[329] | 2127 | else if (getDynamicStatus(j)==soloKey) |
---|
| 2128 | whichKey=j; |
---|
| 2129 | j = next_[j]; //onto next in set |
---|
[328] | 2130 | } |
---|
| 2131 | if (getStatus(iSet)==ClpSimplex::basic&&numberActive) |
---|
| 2132 | numberActive++; |
---|
| 2133 | if (numberActive) { |
---|
| 2134 | assert (numberActive>1); |
---|
| 2135 | int iRow = numberActiveSets_+numberStaticRows_; |
---|
| 2136 | rowSolution[iRow]=0.0; |
---|
| 2137 | double lowerValue; |
---|
| 2138 | if (lowerSet_[iSet]>-1.0e20) |
---|
| 2139 | lowerValue = lowerSet_[iSet]; |
---|
| 2140 | else |
---|
| 2141 | lowerValue=-COIN_DBL_MAX; |
---|
| 2142 | double upperValue; |
---|
| 2143 | if (upperSet_[iSet]<1.0e20) |
---|
| 2144 | upperValue = upperSet_[iSet]; |
---|
| 2145 | else |
---|
| 2146 | upperValue=COIN_DBL_MAX; |
---|
| 2147 | rowLower[iRow]=lowerValue; |
---|
| 2148 | rowUpper[iRow]=upperValue; |
---|
| 2149 | if (getStatus(iSet)==ClpSimplex::basic) { |
---|
| 2150 | model_->setRowStatus(iRow,ClpSimplex::basic); |
---|
| 2151 | rowSolution[iRow]=0.0; |
---|
| 2152 | } else if (getStatus(iSet)==ClpSimplex::atLowerBound) { |
---|
| 2153 | model_->setRowStatus(iRow,ClpSimplex::atLowerBound); |
---|
| 2154 | rowSolution[iRow]=lowerValue; |
---|
| 2155 | } else { |
---|
| 2156 | model_->setRowStatus(iRow,ClpSimplex::atUpperBound); |
---|
| 2157 | rowSolution[iRow]=upperValue; |
---|
| 2158 | } |
---|
| 2159 | j= startSet_[iSet]; |
---|
| 2160 | while (j>=0) { |
---|
| 2161 | DynamicStatus status = getDynamicStatus(j); |
---|
| 2162 | if (status==inSmall) { |
---|
| 2163 | int numberThis = startColumn_[j+1]-startColumn_[j]+1; |
---|
| 2164 | if (numberElements+numberThis>numberElements_) { |
---|
| 2165 | // need to redo |
---|
[399] | 2166 | numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis); |
---|
[328] | 2167 | matrix_->reserve(lastDynamic_,numberElements_); |
---|
| 2168 | element = matrix_->getMutableElements(); |
---|
| 2169 | row = matrix_->getMutableIndices(); |
---|
| 2170 | // these probably okay but be safe |
---|
| 2171 | startColumn = matrix_->getMutableVectorStarts(); |
---|
| 2172 | length = matrix_->getMutableVectorLengths(); |
---|
| 2173 | } |
---|
| 2174 | length[firstAvailable_]=numberThis; |
---|
| 2175 | cost[firstAvailable_]=cost_[j]; |
---|
| 2176 | CoinBigIndex base = startColumn_[j]; |
---|
| 2177 | for (int k=0;k<numberThis-1;k++) { |
---|
| 2178 | row[numberElements]=row_[base+k]; |
---|
| 2179 | element[numberElements++]=element_[base+k]; |
---|
| 2180 | } |
---|
| 2181 | row[numberElements]=iRow; |
---|
| 2182 | element[numberElements++]=1.0; |
---|
| 2183 | id_[firstAvailable_-firstDynamic_]=j; |
---|
| 2184 | solution[firstAvailable_]=0.0; |
---|
| 2185 | model_->setStatus(firstAvailable_,ClpSimplex::basic); |
---|
| 2186 | if (!columnLower_&&!columnUpper_) { |
---|
| 2187 | columnLower[firstAvailable_] = 0.0; |
---|
| 2188 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 2189 | } else { |
---|
| 2190 | if (columnLower_) |
---|
| 2191 | columnLower[firstAvailable_] = columnLower_[j]; |
---|
| 2192 | else |
---|
| 2193 | columnLower[firstAvailable_] = 0.0; |
---|
| 2194 | if (columnUpper_) |
---|
| 2195 | columnUpper[firstAvailable_] = columnUpper_[j]; |
---|
| 2196 | else |
---|
| 2197 | columnUpper[firstAvailable_] = COIN_DBL_MAX; |
---|
| 2198 | if (status==atLowerBound) { |
---|
| 2199 | solution[firstAvailable_]=columnLower[firstAvailable_]; |
---|
| 2200 | } else { |
---|
| 2201 | solution[firstAvailable_]=columnUpper[firstAvailable_]; |
---|
| 2202 | } |
---|
| 2203 | } |
---|
| 2204 | firstAvailable_++; |
---|
| 2205 | startColumn[firstAvailable_]=numberElements; |
---|
| 2206 | } |
---|
| 2207 | j = next_[j]; //onto next in set |
---|
| 2208 | } |
---|
[329] | 2209 | model_->setRowStatus(numberActiveSets_+numberStaticRows_,getStatus(iSet)); |
---|
[328] | 2210 | toIndex_[iSet]=numberActiveSets_; |
---|
| 2211 | fromIndex_[numberActiveSets_++]=iSet; |
---|
| 2212 | } |
---|
[329] | 2213 | assert (toIndex_[iSet]>=0||whichKey>=0); |
---|
| 2214 | keyVariable_[iSet]=whichKey; |
---|
[328] | 2215 | } |
---|
| 2216 | return; |
---|
| 2217 | } |
---|
[330] | 2218 | // Adds in a column to gub structure (called from descendant) |
---|
| 2219 | int |
---|
| 2220 | ClpDynamicMatrix::addColumn(int numberEntries,const int * row, const float * element, |
---|
| 2221 | float cost, float lower, float upper,int iSet, |
---|
| 2222 | DynamicStatus status) |
---|
| 2223 | { |
---|
| 2224 | // check if already in |
---|
| 2225 | int j=startSet_[iSet]; |
---|
| 2226 | while (j>=0) { |
---|
| 2227 | if (startColumn_[j+1]-startColumn_[j]==numberEntries) { |
---|
| 2228 | const int * row2 = row_+startColumn_[j]; |
---|
| 2229 | const float * element2 = element_ + startColumn_[j]; |
---|
| 2230 | bool same=true; |
---|
| 2231 | for (int k=0;k<numberEntries;k++) { |
---|
| 2232 | if (row[k]!=row2[k]||element[k]!=element2[k]) { |
---|
| 2233 | same=false; |
---|
| 2234 | break; |
---|
| 2235 | } |
---|
| 2236 | } |
---|
| 2237 | if (same) { |
---|
| 2238 | bool odd=false; |
---|
| 2239 | if (cost!=cost_[j]) |
---|
| 2240 | odd =true; |
---|
| 2241 | if (columnLower_&&lower!=columnLower_[j]) |
---|
| 2242 | odd=true; |
---|
| 2243 | if (columnUpper_&&upper!=columnUpper_[j]) |
---|
| 2244 | odd=true; |
---|
| 2245 | if (odd) |
---|
| 2246 | printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n", |
---|
| 2247 | cost,lower,upper,cost_[j], |
---|
| 2248 | columnLower_ ? columnLower_[j] : 0.0, |
---|
| 2249 | columnUpper_ ? columnUpper_[j] : 1.0e100); |
---|
| 2250 | else |
---|
| 2251 | return j; |
---|
| 2252 | } |
---|
| 2253 | } |
---|
| 2254 | j = next_[j]; |
---|
| 2255 | } |
---|
| 2256 | |
---|
| 2257 | if (numberGubColumns_==maximumGubColumns_|| |
---|
| 2258 | startColumn_[numberGubColumns_]+numberEntries>maximumElements_) { |
---|
| 2259 | CoinBigIndex j; |
---|
| 2260 | int i; |
---|
| 2261 | int put=0; |
---|
| 2262 | int numberElements=0; |
---|
| 2263 | CoinBigIndex start=0; |
---|
| 2264 | // compress - leave ones at ub and basic |
---|
| 2265 | int * which = new int [numberGubColumns_]; |
---|
| 2266 | for (i=0;i<numberGubColumns_;i++) { |
---|
| 2267 | CoinBigIndex end=startColumn_[i+1]; |
---|
| 2268 | if (getDynamicStatus(i)!=atLowerBound) { |
---|
| 2269 | // keep in |
---|
| 2270 | for (j=start;j<end;j++) { |
---|
| 2271 | row_[numberElements]=row_[j]; |
---|
| 2272 | element_[numberElements++]=element_[j]; |
---|
| 2273 | } |
---|
| 2274 | startColumn_[put+1]=numberElements; |
---|
| 2275 | cost_[put]=cost_[i]; |
---|
| 2276 | if (columnLower_) |
---|
| 2277 | columnLower_[put]=columnLower_[i]; |
---|
| 2278 | if (columnUpper_) |
---|
| 2279 | columnUpper_[put]=columnUpper_[i]; |
---|
| 2280 | dynamicStatus_[put]=dynamicStatus_[i]; |
---|
| 2281 | id_[put]=id_[i]; |
---|
| 2282 | which[i]=put; |
---|
| 2283 | put++; |
---|
| 2284 | } else { |
---|
| 2285 | // out |
---|
| 2286 | which[i]=-1; |
---|
| 2287 | } |
---|
| 2288 | start=end; |
---|
| 2289 | } |
---|
| 2290 | // now redo startSet_ and next_ |
---|
| 2291 | int * newNext = new int [maximumGubColumns_]; |
---|
| 2292 | for (int jSet=0;jSet<numberSets_;jSet++) { |
---|
| 2293 | int sequence = startSet_[jSet]; |
---|
| 2294 | while (which[sequence]<0) { |
---|
| 2295 | // out |
---|
| 2296 | assert (next_[sequence]>=0); |
---|
| 2297 | sequence=next_[sequence]; |
---|
| 2298 | } |
---|
| 2299 | startSet_[jSet]=which[sequence]; |
---|
| 2300 | int last = which[sequence]; |
---|
| 2301 | while (next_[sequence]>=0) { |
---|
| 2302 | sequence = next_[sequence]; |
---|
| 2303 | if(which[sequence]>=0) { |
---|
| 2304 | // keep |
---|
| 2305 | newNext[last]=which[sequence]; |
---|
| 2306 | last=which[sequence]; |
---|
| 2307 | } |
---|
| 2308 | } |
---|
| 2309 | newNext[last]=-jSet-1; |
---|
| 2310 | } |
---|
| 2311 | delete [] next_; |
---|
| 2312 | next_=newNext; |
---|
| 2313 | delete [] which; |
---|
| 2314 | abort(); |
---|
| 2315 | } |
---|
| 2316 | CoinBigIndex start = startColumn_[numberGubColumns_]; |
---|
| 2317 | memcpy(row_+start,row,numberEntries*sizeof(int)); |
---|
| 2318 | memcpy(element_+start,element,numberEntries*sizeof(float)); |
---|
| 2319 | startColumn_[numberGubColumns_+1]=start+numberEntries; |
---|
| 2320 | cost_[numberGubColumns_]=cost; |
---|
| 2321 | if (columnLower_) |
---|
| 2322 | columnLower_[numberGubColumns_]=lower; |
---|
| 2323 | else |
---|
| 2324 | assert (!lower); |
---|
| 2325 | if (columnUpper_) |
---|
| 2326 | columnUpper_[numberGubColumns_]=upper; |
---|
| 2327 | else |
---|
| 2328 | assert (upper>1.0e20); |
---|
| 2329 | setDynamicStatus(numberGubColumns_,status); |
---|
| 2330 | // Do next_ |
---|
| 2331 | j=startSet_[iSet]; |
---|
| 2332 | startSet_[iSet]=numberGubColumns_; |
---|
| 2333 | next_[numberGubColumns_]=j; |
---|
| 2334 | numberGubColumns_++; |
---|
| 2335 | return numberGubColumns_-1; |
---|
| 2336 | } |
---|
| 2337 | // Returns which set a variable is in |
---|
| 2338 | int |
---|
| 2339 | ClpDynamicMatrix::whichSet (int sequence) const |
---|
| 2340 | { |
---|
| 2341 | while (next_[sequence]>=0) |
---|
| 2342 | sequence = next_[sequence]; |
---|
| 2343 | int iSet = - next_[sequence]-1; |
---|
| 2344 | return iSet; |
---|
| 2345 | } |
---|